{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":84969,"databundleVersionId":10033515,"sourceType":"competition"},{"sourceId":9862305,"sourceType":"datasetVersion","datasetId":6052780},{"sourceId":9867543,"sourceType":"datasetVersion","datasetId":6040935},{"sourceId":10225818,"sourceType":"datasetVersion","datasetId":6321984},{"sourceId":10445850,"sourceType":"datasetVersion","datasetId":6465904},{"sourceId":206640467,"sourceType":"kernelVersion"},{"sourceId":211097053,"sourceType":"kernelVersion"}],"dockerImageVersionId":30823,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"from IPython.display import clear_output\n!tar xfvz /kaggle/input/ultralytics-for-offline-install/archive.tar.gz\n!pip install --no-index --find-links=./packages ultralytics\n!rm -rf ./packages\ntry:\n    import zarr\nexcept: \n    !cp -r '/kaggle/input/hengck-czii-cryo-et-01/wheel_file' '/kaggle/working/'\n    !pip install /kaggle/working/wheel_file/asciitree-0.3.3/asciitree-0.3.3\n    !pip install --no-index --find-links=/kaggle/working/wheel_file zarr\n    !pip install --no-index --find-links=/kaggle/working/wheel_file connected-components-3d\nfrom typing import List, Tuple, Union\ndeps_path = '/kaggle/input/czii-cryoet-dependencies'\n! pip install -q --no-index --find-links {deps_path} --requirement {deps_path}/requirements.txt\nimport lightning.pytorch as pl\nfrom datetime import datetime\nimport pytz\nimport sys\nsys.path.append('/kaggle/input/hengck-czii-cryo-et-01')\nfrom czii_helper import *\nfrom dataset import *\nfrom model2 import *\nclear_output()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:33:38.197229Z","iopub.execute_input":"2025-01-12T12:33:38.197548Z","iopub.status.idle":"2025-01-12T12:34:09.411610Z","shell.execute_reply.started":"2025-01-12T12:33:38.197517Z","shell.execute_reply":"2025-01-12T12:34:09.410525Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport glob\nimport time\nimport sys\nimport warnings\nimport math\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport cv2\nimport torch\nfrom tqdm import tqdm\nfrom ultralytics import YOLO\nimport zarr\nfrom scipy.spatial import cKDTree\nfrom collections import defaultdict","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:33.761003Z","iopub.execute_input":"2025-01-12T12:35:33.761362Z","iopub.status.idle":"2025-01-12T12:35:36.131934Z","shell.execute_reply.started":"2025-01-12T12:35:33.761338Z","shell.execute_reply":"2025-01-12T12:35:36.131262Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"model_path = '/kaggle/input/czii-yolo-l-trained-with-synthetic-data/best_synthetic.pt'\nmodel = YOLO(model_path)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:38.203139Z","iopub.execute_input":"2025-01-12T12:35:38.203483Z","iopub.status.idle":"2025-01-12T12:35:38.857292Z","shell.execute_reply.started":"2025-01-12T12:35:38.203458Z","shell.execute_reply":"2025-01-12T12:35:38.856270Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"runs_path = '/kaggle/input/czii-cryo-et-object-identification/test/static/ExperimentRuns/*'\nruns = sorted(glob.glob(runs_path))\nruns = [os.path.basename(run) for run in runs]\nsp = len(runs)//2\nruns1 = runs[:sp]\nruns1[:5]\n\n#add by @minfuka\nruns2 = runs[sp:]\nruns2[:5]\n\n#add by @minfuka\nassert torch.cuda.device_count() == 2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:40.388532Z","iopub.execute_input":"2025-01-12T12:35:40.388857Z","iopub.status.idle":"2025-01-12T12:35:40.426070Z","shell.execute_reply.started":"2025-01-12T12:35:40.388830Z","shell.execute_reply":"2025-01-12T12:35:40.425342Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"particle_names = [\n    'apo-ferritin',\n    'beta-amylase',\n    'beta-galactosidase',\n    'ribosome',\n    'thyroglobulin',\n    'virus-like-particle'\n]\n\nparticle_to_index = {\n    'apo-ferritin': 0,\n    'beta-amylase': 1,\n    'beta-galactosidase': 2,\n    'ribosome': 3,\n    'thyroglobulin': 4,\n    'virus-like-particle': 5\n}\n\nindex_to_particle = {index: name for name, index in particle_to_index.items()}\n\nparticle_radius = {\n    'apo-ferritin': 60,\n    'beta-amylase': 65,\n    'beta-galactosidase': 90,\n    'ribosome': 150,\n    'thyroglobulin': 130,\n    'virus-like-particle': 135,\n}\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:42.453542Z","iopub.execute_input":"2025-01-12T12:35:42.453868Z","iopub.status.idle":"2025-01-12T12:35:42.458664Z","shell.execute_reply.started":"2025-01-12T12:35:42.453840Z","shell.execute_reply":"2025-01-12T12:35:42.457863Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# add by @sesasj\nclass UnionFind:\n    def __init__(self, size):\n        self.parent = np.arange(size)\n        self.rank = np.zeros(size, dtype=int)\n\n    def find(self, u):\n        if self.parent[u] != u:\n            self.parent[u] = self.find(self.parent[u])  \n        return self.parent[u]\n\n    def union(self, u, v):\n        u_root = self.find(u)\n        v_root = self.find(v)\n        if u_root == v_root:\n            return\n            \n        if self.rank[u_root] < self.rank[v_root]:\n            self.parent[u_root] = v_root\n        else:\n            self.parent[v_root] = u_root\n            if self.rank[u_root] == self.rank[v_root]:\n                self.rank[u_root] += 1\n\nclass PredictionAggregator:\n    def __init__(self, first_conf=0.2, conf_coef=0.75):\n        self.first_conf = first_conf\n        self.conf_coef = conf_coef\n        self.particle_confs = np.array([0.5, 0.0, 0.2, 0.5, 0.2, 0.5])\n        \n    def convert_to_8bit(self, volume):\n        lower, upper = np.percentile(volume, (0.5, 99.5))\n        clipped = np.clip(volume, lower, upper)\n        scaled = ((clipped - lower) / (upper - lower + 1e-12) * 255).astype(np.uint8)\n        return scaled\n\n    def make_predictions(self, run_id, model, device_no):\n        volume_path = f'/kaggle/input/czii-cryo-et-object-identification/test/static/ExperimentRuns/{run_id}/VoxelSpacing10.000/denoised.zarr'\n        volume = zarr.open(volume_path, mode='r')[0]\n        volume_8bit = self.convert_to_8bit(volume)\n        num_slices = volume_8bit.shape[0]\n\n        detections = {\n            'particle_type': [],\n            'confidence': [],\n            'x': [],\n            'y': [],\n            'z': []\n        }\n\n        for slice_idx in range(num_slices):\n            \n            img = volume_8bit[slice_idx]\n            input_image = cv2.resize(np.stack([img]*3, axis=-1), (640, 640))\n\n            results = model.predict(\n                input_image,\n                save=False,\n                imgsz=640,\n                conf=self.first_conf,\n                device=device_no,\n                batch=1,\n                verbose=False,\n            )\n\n            for result in results:\n                boxes = result.boxes\n                if boxes is None:\n                    continue\n                cls = boxes.cls.cpu().numpy().astype(int)\n                conf = boxes.conf.cpu().numpy()\n                xyxy = boxes.xyxy.cpu().numpy()\n\n                xc = ((xyxy[:, 0] + xyxy[:, 2]) / 2.0) * 10 * (63/64) # 63/64 because of the resize\n                yc = ((xyxy[:, 1] + xyxy[:, 3]) / 2.0) * 10 * (63/64)\n                zc = np.full(xc.shape, slice_idx * 10 + 5)\n\n                particle_types = [index_to_particle[c] for c in cls]\n\n                detections['particle_type'].extend(particle_types)\n                detections['confidence'].extend(conf)\n                detections['x'].extend(xc)\n                detections['y'].extend(yc)\n                detections['z'].extend(zc)\n\n        if not detections['particle_type']:\n            return pd.DataFrame()  \n\n        particle_types = np.array(detections['particle_type'])\n        confidences = np.array(detections['confidence'])\n        xs = np.array(detections['x'])\n        ys = np.array(detections['y'])\n        zs = np.array(detections['z'])\n\n        aggregated_data = []\n\n        for idx, particle in enumerate(particle_names):\n            if particle == 'beta-amylase':\n                continue \n\n            mask = (particle_types == particle)\n            if not np.any(mask):\n                continue  \n                \n            particle_confidences = confidences[mask]\n            particle_xs = xs[mask]\n            particle_ys = ys[mask]\n            particle_zs = zs[mask]\n            # -------------modified by @sersasj ------------------------\n            coords = np.vstack((particle_xs, particle_ys, particle_zs)).T\n\n           \n            z_distance = 30 # How many slices can you \"jump\" to aggregate predictions 10 = 1, 20 = 2...\n            xy_distance = 20 # xy_tol_p2 in original code by ITK8191\n            \n            max_distance = math.sqrt(z_distance**2 + xy_distance**2)\n            tree = cKDTree(coords)            \n            pairs = tree.query_pairs(r=max_distance, p=2)\n\n            \n            uf = UnionFind(len(coords))\n            \n            coords_xy = coords[:, :2]\n            coords_z = coords[:, 2]\n            for u, v in pairs:\n                z_diff = abs(coords_z[u] - coords_z[v])\n                if z_diff > z_distance:\n                    continue  \n\n                xy_diff = np.linalg.norm(coords_xy[u] - coords_xy[v])\n                if xy_diff > xy_distance:\n                    continue  \n\n                uf.union(u, v)\n\n            roots = np.array([uf.find(i) for i in range(len(coords))])\n            unique_roots, inverse_indices, counts = np.unique(roots, return_inverse=True, return_counts=True)\n            conf_sums = np.bincount(inverse_indices, weights=particle_confidences)\n            \n            aggregated_confidences = conf_sums / (counts ** self.conf_coef)\n            cluster_per_particle = [4,1,2,9,4,8] # Update\n            valid_clusters = (counts >= cluster_per_particle[idx]) & (aggregated_confidences > self.particle_confs[idx])\n\n            if not np.any(valid_clusters):\n                continue  \n\n            cluster_ids = unique_roots[valid_clusters]\n\n            centers_x = np.bincount(inverse_indices, weights=particle_xs) / counts\n            centers_y = np.bincount(inverse_indices, weights=particle_ys) / counts\n            centers_z = np.bincount(inverse_indices, weights=particle_zs) / counts\n\n            centers_x = centers_x[valid_clusters]\n            centers_y = centers_y[valid_clusters]\n            centers_z = centers_z[valid_clusters]\n\n            aggregated_df = pd.DataFrame({\n                'experiment': [run_id] * len(centers_x),\n                'particle_type': [particle] * len(centers_x),\n                'x': centers_x,\n                'y': centers_y,\n                'z': centers_z\n            })\n\n            aggregated_data.append(aggregated_df)\n\n        if aggregated_data:\n            return pd.concat(aggregated_data, axis=0)\n        else:\n            return pd.DataFrame()  ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:44.940865Z","iopub.execute_input":"2025-01-12T12:35:44.941242Z","iopub.status.idle":"2025-01-12T12:35:44.961211Z","shell.execute_reply.started":"2025-01-12T12:35:44.941209Z","shell.execute_reply":"2025-01-12T12:35:44.960123Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# instance main class\naggregator = PredictionAggregator(first_conf=0.19,  conf_coef=0.34) #Update\naggregated_results = []\n#add by @minfuka\nfrom concurrent.futures import ProcessPoolExecutor #add by @minfuka\n\n#add by @minfuka\ndef inference(runs, model, device_no):\n    subs = []\n    for r in tqdm(runs, total=len(runs)):\n        df = aggregator.make_predictions(r, model, device_no)\n        subs.append(df)\n    \n    return subs\nstart_time = time.time()\n\nwith ProcessPoolExecutor(max_workers=2) as executor:\n    results = list(executor.map(inference, (runs1, runs2), (model, model), (\"0\", \"1\")))\n\n\nend_time = time.time()\n\nestimated_total_time = (end_time - start_time) / len(runs) * 500  \nprint(f'estimated total prediction time for 500 runs: {estimated_total_time:.4f} seconds')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:35:48.650704Z","iopub.execute_input":"2025-01-12T12:35:48.651067Z","iopub.status.idle":"2025-01-12T12:36:14.629800Z","shell.execute_reply.started":"2025-01-12T12:35:48.651027Z","shell.execute_reply":"2025-01-12T12:36:14.628571Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#change by @minfuka\nsubmission0 = pd.concat(results[0])\nsubmission1 = pd.concat(results[1])\nsubmission_ = pd.concat([submission0, submission1]).reset_index(drop=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:36:17.544432Z","iopub.execute_input":"2025-01-12T12:36:17.544790Z","iopub.status.idle":"2025-01-12T12:36:17.552826Z","shell.execute_reply.started":"2025-01-12T12:36:17.544758Z","shell.execute_reply":"2025-01-12T12:36:17.551746Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_.insert(0, 'id', range(len(submission_)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:36:20.522960Z","iopub.execute_input":"2025-01-12T12:36:20.523380Z","iopub.status.idle":"2025-01-12T12:36:20.532130Z","shell.execute_reply.started":"2025-01-12T12:36:20.523349Z","shell.execute_reply":"2025-01-12T12:36:20.531013Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Model(pl.LightningModule):\n    def __init__(\n        self, \n        spatial_dims: int = 3,\n        in_channels: int = 1,\n        out_channels: int = 7,\n        channels: Union[Tuple[int, ...], List[int]] = (48, 64, 80, 80),\n        strides: Union[Tuple[int, ...], List[int]] = (2, 2, 1),\n        num_res_units: int = 1,\n    ):\n        super().__init__()\n        self.save_hyperparameters()\n        self.model = UNet(\n            spatial_dims=self.hparams.spatial_dims,\n            in_channels=self.hparams.in_channels,\n            out_channels=self.hparams.out_channels,\n            channels=self.hparams.channels,\n            strides=self.hparams.strides,\n            num_res_units=self.hparams.num_res_units,\n        )\n    def forward(self, x):\n        return self.model(x)\n# 设置训练超参数\nchannels = (48, 64, 80, 80)\nstrides_pattern = (2, 2, 1)\nnum_res_units = 1\ndef extract_3d_patches_minimal_overlap(arrays: List[np.ndarray], patch_size: int) -> Tuple[List[np.ndarray], List[Tuple[int, int, int]]]:\n    if not arrays or not isinstance(arrays, list):\n        raise ValueError(\"Input must be a non-empty list of arrays\")\n    \n    # Verify all arrays have the same shape\n    shape = arrays[0].shape\n    if not all(arr.shape == shape for arr in arrays):\n        raise ValueError(\"All input arrays must have the same shape\")\n    \n    if patch_size > min(shape):\n        raise ValueError(f\"patch_size ({patch_size}) must be smaller than smallest dimension {min(shape)}\")\n    \n    m, n, l = shape\n    patches = []\n    coordinates = []\n    \n    # Calculate starting positions for each dimension\n    x_starts = calculate_patch_starts(m, patch_size)\n    y_starts = calculate_patch_starts(n, patch_size)\n    z_starts = calculate_patch_starts(l, patch_size)\n    \n    # Extract patches from each array\n    for arr in arrays:\n        for x in x_starts:\n            for y in y_starts:\n                for z in z_starts:\n                    patch = arr[\n                        x:x + patch_size,\n                        y:y + patch_size,\n                        z:z + patch_size\n                    ]\n                    patches.append(patch)\n                    coordinates.append((x, y, z))\n    \n    return patches, coordinates\ndef reconstruct_array(patches: List[np.ndarray], \n                     coordinates: List[Tuple[int, int, int]], \n                     original_shape: Tuple[int, int, int]) -> np.ndarray:\n    reconstructed = np.zeros(original_shape, dtype=np.int64)  # To track overlapping regions\n    \n    patch_size = patches[0].shape[0]\n    \n    for patch, (x, y, z) in zip(patches, coordinates):\n        reconstructed[\n            x:x + patch_size,\n            y:y + patch_size,\n            z:z + patch_size\n        ] = patch\n        \n    \n    return reconstructed\ndef calculate_patch_starts(dimension_size: int, patch_size: int) -> List[int]:\n    if dimension_size <= patch_size:\n        return [0]\n        \n    # Calculate number of patches needed\n    n_patches = np.ceil(dimension_size / patch_size)\n    \n    if n_patches == 1:\n        return [0]\n    \n    # Calculate overlap\n    total_overlap = (n_patches * patch_size - dimension_size) / (n_patches - 1)\n    \n    # Generate starting positions\n    positions = []\n    for i in range(int(n_patches)):\n        pos = int(i * (patch_size - total_overlap))\n        if pos + patch_size > dimension_size:\n            pos = dimension_size - patch_size\n        if pos not in positions:  # Avoid duplicates\n            positions.append(pos)\n    \n    return positions\nimport pandas as pd\n\ndef dict_to_df(coord_dict, experiment_name):\n    # Create lists to store data\n    all_coords = []\n    all_labels = []\n    \n    # Process each label and its coordinates\n    for label, coords in coord_dict.items():\n        all_coords.append(coords)\n        all_labels.extend([label] * len(coords))\n    \n    # Concatenate all coordinates\n    all_coords = np.vstack(all_coords)\n    \n    df = pd.DataFrame({\n        'experiment': experiment_name,\n        'particle_type': all_labels,\n        'x': all_coords[:, 0],\n        'y': all_coords[:, 1],\n        'z': all_coords[:, 2]\n    })\n\n    \n    return df\nfrom typing import List, Tuple, Union\nimport numpy as np\nimport torch\nfrom monai.data import DataLoader, Dataset, CacheDataset, decollate_batch\nfrom monai.transforms import (\n    Compose, \n    EnsureChannelFirstd, \n    Orientationd,  \n    AsDiscrete,  \n    RandFlipd, \n    RandRotate90d, \n    NormalizeIntensityd,\n    RandCropByLabelClassesd,\n)\nTRAIN_DATA_DIR = \"/kaggle/input/create-numpy-dataset-exp-name\"\nimport json\ncopick_config_path = TRAIN_DATA_DIR + \"/copick.config\"\n\nwith open(copick_config_path) as f:\n    copick_config = json.load(f)\n\ncopick_config['static_root'] = '/kaggle/input/czii-cryo-et-object-identification/test/static'\n\ncopick_test_config_path = 'copick_test.config'\n\nwith open(copick_test_config_path, 'w') as outfile:\n    json.dump(copick_config, outfile)\nimport copick\n\nroot = copick.from_file(copick_test_config_path)\n\ncopick_user_name = \"copickUtils\"\ncopick_segmentation_name = \"paintedPicks\"\nvoxel_size = 10\ntomo_type = \"denoised\"\ninference_transforms = Compose([\n    EnsureChannelFirstd(keys=[\"image\"], channel_dim=\"no_channel\"),\n    NormalizeIntensityd(keys=\"image\"),\n    Orientationd(keys=[\"image\"], axcodes=\"RAS\")\n])\nimport cc3d\n\nid_to_name = {1: \"apo-ferritin\", \n              2: \"beta-amylase\",\n              3: \"beta-galactosidase\", \n              4: \"ribosome\", \n              5: \"thyroglobulin\", \n              6: \"virus-like-particle\"}\nBLOB_THRESHOLD = 200\nCERTAINTY_THRESHOLD = 0.05\n\nclasses = [1, 2, 3, 4, 5, 6]\nimport torch\nimport numpy as np\nimport pandas as pd\nimport cc3d\nfrom monai.data import CacheDataset\nfrom monai.transforms import Compose, EnsureType\nfrom torch import nn\nfrom tqdm import tqdm\n\n# 假设你有一个函数来加载模型权重\ndef load_models(model_paths):\n    models = []\n    for model_path in model_paths:\n        model = torch.load(model_path)\n        model.to('cuda')\n        model.eval()\n        models.append(model)\n    return models\n\n# 假设你有3个模型的路径\nmodel_paths = [\n    '/kaggle/input/weight-all/weight/val_on_TS_6_6_best_model.pt'\n]\n\n# 加载模型\nmodels = load_models(model_paths)\ndef ensemble_prediction_tta(models, input_tensor, threshold=0.5):\n    probs_list = []\n    data_copy0 = input_tensor.clone()\n    data_copy0=torch.flip(data_copy0, dims=[2])\n    data_copy1 = input_tensor.clone()\n    data_copy1=torch.flip(data_copy1, dims=[3])\n    data_copy2 = input_tensor.clone()\n    data_copy2=torch.flip(data_copy2, dims=[4])\n    data_copy3 = input_tensor.clone()\n    data_copy3 = data_copy3.rot90(1, dims=[3, 4])\n    with torch.no_grad():\n        model_output0 = model(input_tensor)\n        model_output1 = model(data_copy0)\n        model_output1=torch.flip(model_output1, dims=[2])\n        model_output2 = model(data_copy1)\n        model_output2=torch.flip(model_output2, dims=[3])\n        model_output3 = model(data_copy2)\n        model_output3=torch.flip(model_output3, dims=[4])\n        probs0 = torch.softmax(model_output0[0], dim=0)\n        probs1 = torch.softmax(model_output1[0], dim=0)\n        probs2 = torch.softmax(model_output2[0], dim=0)\n        probs3 = torch.softmax(model_output3[0], dim=0)\n        probs_list.append(probs0)\n        probs_list.append(probs1)\n        probs_list.append(probs2)\n        probs_list.append(probs3)\n    avg_probs = torch.mean(torch.stack(probs_list), dim=0)\n    thresh_probs = avg_probs > threshold\n    _, max_classes = thresh_probs.max(dim=0)\n    return max_classes\nsub=[]\nfor model in models:\n    with torch.no_grad():\n        location_df = []\n        for run in root.runs:\n            tomo = run.get_voxel_spacing(10)\n            tomo = tomo.get_tomogram(tomo_type).numpy()\n            tomo_patches, coordinates = extract_3d_patches_minimal_overlap([tomo], 96)\n            tomo_patched_data = [{\"image\": img} for img in tomo_patches]\n            tomo_ds = CacheDataset(data=tomo_patched_data, transform=inference_transforms, cache_rate=1.0)\n            pred_masks = []\n            for i in tqdm(range(len(tomo_ds))):\n                input_tensor = tomo_ds[i]['image'].unsqueeze(0).to(\"cuda\")\n                max_classes = ensemble_prediction_tta(models, input_tensor, threshold=CERTAINTY_THRESHOLD)\n                pred_masks.append(max_classes.cpu().numpy())\n            reconstructed_mask = reconstruct_array(pred_masks, coordinates, tomo.shape)\n            location = {}\n            for c in classes:\n                cc = cc3d.connected_components(reconstructed_mask == c)\n                stats = cc3d.statistics(cc)\n                zyx = stats['centroids'][1:] * 10.012444  # 转换单位\n                zyx_large = zyx[stats['voxel_counts'][1:] > BLOB_THRESHOLD]\n                xyz = np.ascontiguousarray(zyx_large[:, ::-1])\n                location[id_to_name[c]] = xyz\n            df = dict_to_df(location, run.name)\n            location_df.append(df)\n        location_df = pd.concat(location_df)\n        location_df.insert(loc=0, column='id', value=np.arange(len(location_df)))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:36:22.391044Z","iopub.execute_input":"2025-01-12T12:36:22.391409Z","iopub.status.idle":"2025-01-12T12:37:34.399043Z","shell.execute_reply.started":"2025-01-12T12:36:22.391378Z","shell.execute_reply":"2025-01-12T12:37:34.398011Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nfrom sklearn.cluster import DBSCAN\n\n# 假设sub已经给定，拼接DataFrame\ndf = pd.concat([submission_,location_df], ignore_index=True)\n# 粒子半径映射\nparticle_names = ['apo-ferritin', 'beta-amylase', 'beta-galactosidase', 'ribosome', 'thyroglobulin', 'virus-like-particle']\nparticle_radius = {\n    'apo-ferritin': 60,\n    'beta-amylase': 65,\n    'beta-galactosidase': 90,\n    'ribosome': 150,\n    'thyroglobulin': 130,\n    'virus-like-particle': 135,\n}\n\nfinal = []  # 用于存储最终的结果\nfor pidx, p in enumerate(particle_names):\n    # 筛选出该粒子类型的所有点\n    pdf = df[df['particle_type'] == p].reset_index(drop=True)\n    p_rad = particle_radius[p]\n    \n    # 根据 experiment 分组\n    grouped = pdf.groupby(['experiment'])\n    \n    for exp, group in grouped:\n        group = group.reset_index(drop=True)\n        \n        # 使用DBSCAN进行聚类\n        coords = group[['x', 'y', 'z']].values\n        db = DBSCAN(eps=p_rad, min_samples=2, metric='euclidean').fit(coords)\n        labels = db.labels_\n        \n        # 将聚类结果添加到DataFrame中\n        group['cluster'] = labels\n        \n        # 对每个簇进行处理\n        for cluster_id in np.unique(labels):\n            if cluster_id == -1:\n                continue  # 跳过噪声点\n            \n            cluster_points = group[group['cluster'] == cluster_id]\n            \n            # 计算簇的中心（平均位置）\n            avg_x = cluster_points['x'].mean()\n            avg_y = cluster_points['y'].mean()\n            avg_z = cluster_points['z'].mean()\n            \n            # 更新簇内点的位置\n            group.loc[group['cluster'] == cluster_id, ['x', 'y', 'z']] = avg_x, avg_y, avg_z\n            group = group.drop_duplicates(subset=['x', 'y', 'z'])\n        # 将处理后的数据添加到 final 列表\n        final.append(group)\n\n# 合并处理后的数据\ndf_save = pd.concat(final, ignore_index=True)\ndf_save = df_save.drop(columns=['cluster'])\n\n# 排序按 'experiment' 和 'particle_type' 两列\ndf_save = df_save.sort_values(by=['experiment', 'particle_type']).reset_index(drop=True)\n\n# 重新生成 'id' 列，从 1 开始\ndf_save['id'] = np.arange(0, len(df_save))\n\n# 输出结果到 CSV\ndf_save.to_csv('submission.csv', index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-01-12T12:37:45.927944Z","iopub.execute_input":"2025-01-12T12:37:45.929058Z","iopub.status.idle":"2025-01-12T12:37:47.976347Z","shell.execute_reply.started":"2025-01-12T12:37:45.929018Z","shell.execute_reply":"2025-01-12T12:37:47.975584Z"}},"outputs":[],"execution_count":null}]}