{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"please see: \nhttps://blog.redbrickai.com/blog-posts/introduction-to-dicom-coordinate  \n\n- it explains the 2d and 3d coordinate systems used by dicom\n- the math for 2d to 3d affine transformation","metadata":{}},{"cell_type":"code","source":"if 1: \n    #natural sort\n    !pip install natsort","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:45:57.972975Z","iopub.execute_input":"2024-07-30T03:45:57.973389Z","iopub.status.idle":"2024-07-30T03:46:16.043154Z","shell.execute_reply.started":"2024-07-30T03:45:57.973354Z","shell.execute_reply":"2024-07-30T03:46:16.041705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kaggle_dir = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification'\n\nimport matplotlib\nimport matplotlib.pyplot as plt\n\nimport glob\nimport pydicom\nfrom natsort import natsorted, ns\nimport numpy as np\nimport pandas as pd\n\npd.options.mode.chained_assignment = None  # default='warn'\n# disable SettingWithCopyWarning\n\n\n#helper\nclass dotdict(dict):\n    __setattr__ = dict.__setitem__\n    __delattr__ = dict.__delitem__\n\n    def __getattr__(self, name):\n        try:\n            return self[name]\n        except KeyError:\n            raise AttributeError(name)\n            \nprint('import ok')","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:16.045611Z","iopub.execute_input":"2024-07-30T03:46:16.046090Z","iopub.status.idle":"2024-07-30T03:46:16.802172Z","shell.execute_reply.started":"2024-07-30T03:46:16.046051Z","shell.execute_reply":"2024-07-30T03:46:16.801057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#dicom reader\n\ndef normalise_to_8bit(x, lower=0.1, upper=99.9): # 1, 99 #0.05, 99.5 #0, 100\n    lower, upper = np.percentile(x, (lower, upper))\n    x = np.clip(x, lower, upper)\n    x = x - np.min(x)\n    x = x / np.max(x)\n    return (x * 255).astype(np.uint8)\n\n# discontinous volume will be splitted to continous chunks\ndef load_and_split_mri_from_dicom_dir(\n        study_id,\n        series_id,\n        series_description,\n    ):\n    dicom_dir= f'{kaggle_dir}/train_images/{study_id}/{series_id}'\n\n    #---\n    dicom_file = natsorted(glob.glob( f'{dicom_dir}/*.dcm'))\n    instance_number = [int(f.split('/')[-1].split('.')[0]) for f in dicom_file]\n    dicom = [pydicom.dcmread(f) for f in dicom_file]\n\n    dicom_df = []\n    for i,d in zip(instance_number,dicom): #d__.dict__\n        dicom_df.append(\n            dotdict(\n                study_id=study_id,\n                series_id=series_id,\n                series_description=series_description,\n                instance_number=i,\n                #InstanceNumber = d.InstanceNumber,\n                ImagePositionPatient=tuple([float(v) for v in d.ImagePositionPatient]),  # not continous\n                ImageOrientationPatient=tuple([float(v) for v in d.ImageOrientationPatient]),\n                PixelSpacing=tuple([float(v) for v in d.PixelSpacing]),\n                SpacingBetweenSlices=float(d.SpacingBetweenSlices),\n                SliceThickness=float(d.SliceThickness),\n            )\n        )\n    dicom_df = pd.DataFrame(dicom_df)\n    dicom_df = [d for _,d in dicom_df.groupby('ImageOrientationPatient')]\n\n    #sort slice\n    mri=[]\n    for df in dicom_df:\n        position = np.array(df['ImagePositionPatient'].values.tolist())\n        orientation = np.array(df['ImageOrientationPatient'].values.tolist())\n        normal = np.cross(orientation[:,:3], orientation[:,3:])\n        projection = np.sum(normal*position,1) #np.dot(normal, position)\n        df.loc[:,'projection'] = projection\n        df = df.sort_values('projection')\n\n        #todo: assert all slices are continous\n        assert len(df.SliceThickness.unique())==1\n        assert len(df.ImageOrientationPatient.unique())==1\n        assert len(df.ImageOrientationPatient.unique())==1\n        assert len(df.SpacingBetweenSlices.unique())==1\n\n        volume =[\n            dicom[instance_number.index(i)].pixel_array for i in df.instance_number\n        ]\n        volume = np.stack(volume)\n        volume = normalise_to_8bit(volume)\n        mri.append(dotdict(\n            df=df,\n            volume=volume,\n        ))\n    return mri\n\n\n#convert 2d x,y to 3d X,Y,Z for point in label.csv\ndef add_XYZ_to_label_df(study_id_df):\n\n    # add shape W,H; world coords xx,yy,zz\n    for col in ['W','H']:\n        study_id_df.loc[:,col]=0\n    for col in ['xx','yy','zz']:\n        study_id_df.loc[:,col]=0.0\n\n    for t,d in study_id_df.iterrows():\n        #print(d)\n        #print('')\n        dicom_file = f'{kaggle_dir}/train_images/{d.study_id}/{d.series_id}/{d.instance_number}.dcm'\n        dicom = pydicom.dcmread(dicom_file)\n        H,W = dicom.pixel_array.shape\n        sx,sy,sz = [float(v) for v in dicom.ImagePositionPatient]\n        o0, o1, o2, o3, o4, o5, = [float(v) for v in dicom.ImageOrientationPatient]\n        delx,dely = dicom.PixelSpacing\n\n        xx =  o0*delx*d.x + o3*dely*d.y + sx\n        yy =  o1*delx*d.x + o4*dely*d.y + sy\n        zz =  o2*delx*d.x + o5*dely*d.y + sz\n\n        study_id_df.loc[t,'W'] = W\n        study_id_df.loc[t,'H'] = H\n        study_id_df.loc[t,'xx'] = xx\n        study_id_df.loc[t,'yy'] = yy\n        study_id_df.loc[t,'zz'] = zz\n    return study_id_df\n\n#read all mri for one patient\ndef load_for_one(study_id_df):\n    gb = study_id_df.groupby(['series_description','series_id']).agg('first').index #[['series_id', 'series_description']]\n    mri=[]\n    for series_description,series_id in gb:\n        mri += load_and_split_mri_from_dicom_dir(\n            study_id=study_id,\n            series_description=series_description,\n            series_id= series_id\n        )\n    return mri\n\n\n#back project 3D to 2d\ndef backproject_XYZ(xx,yy,zz,mri):\n    r = mri\n    d0 = r.df.iloc[0] \n\n    sx, sy, sz = [float(v) for v in d0.ImagePositionPatient]\n    o0, o1, o2, o3, o4, o5, = [float(v) for v in d0.ImageOrientationPatient]\n    delx, dely = d0.PixelSpacing\n    delz = d0.SpacingBetweenSlices\n\n    ax = np.array([o0,o1,o2])\n    ay = np.array([o3,o4,o5])\n    az = np.cross(ax,ay)\n\n    p = np.array([xx-sx,yy-sy,zz-sz])\n    x = np.dot(ax, p)/delx\n    y = np.dot(ay, p)/dely\n    z = np.dot(az, p)/delz\n    x = int(round(x))\n    y = int(round(y))\n    z = int(round(z))\n\n    D,H,W = r.volume.shape\n    inside = \\\n        (x>=0) & (x<W) &\\\n        (y>=0) & (y<H) &\\\n        (z>=0) & (z<D)\n    if not inside:\n        #print('out-of-bound')\n        return False,0,0,0,0\n\n    n = r.df.instance_number.values[z]\n    return True,x,y,z,n\n        ","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:16.803679Z","iopub.execute_input":"2024-07-30T03:46:16.804152Z","iopub.status.idle":"2024-07-30T03:46:16.839062Z","shell.execute_reply.started":"2024-07-30T03:46:16.804122Z","shell.execute_reply":"2024-07-30T03:46:16.837862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load kaggle csv\ndesc_df = pd.read_csv(f'{kaggle_dir}/train_series_descriptions.csv')\nlabel_df = pd.read_csv(f'{kaggle_dir}/train_label_coordinates.csv')\nlabel_df = label_df.merge(desc_df, on=['study_id', 'series_id'])","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:16.841704Z","iopub.execute_input":"2024-07-30T03:46:16.842084Z","iopub.status.idle":"2024-07-30T03:46:17.046646Z","shell.execute_reply.started":"2024-07-30T03:46:16.842053Z","shell.execute_reply":"2024-07-30T03:46:17.045417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#verify code\n\n'''\nstudy_id\tseries_id\tseries_description\n0\t4003253\t702807833\tSagittal T2/STIR\n1\t4003253\t1054713880\tSagittal T1\n2\t4003253\t2448190387\tAxial T2\n'''\n\nstudy_id = 4003253\nstudy_id_df = label_df[label_df.study_id==study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\nprint(study_id_df.iloc[0])\n\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))\nprint('')\n\n# backproject within the same view as truth point:\n# - given x,y,z of view v from label.csv\n# - project to 3d xx,yy,zz\n# - for the same view v, backproject from xx,yy,zz to vx,vy,vz\n# - check (vx,vy,vz) must be same as (x,y,z)\n\nprint('START VERIFICATION !!!')\nfor t,d in study_id_df.iterrows():\n    print('=================================')\n    print(d)\n    print('')\n    xx, yy, zz = d.xx, d.yy, d.zz\n    \n    found=0\n    for r in mri:\n        d0 = r.df.iloc[0]\n        if not (\n            (d0.study_id ==d.study_id) & \n            (d0.series_id ==d.series_id) &\n            (d.instance_number in r.df.instance_number)\n        ): continue\n        inside,x,y,z,n = backproject_XYZ(xx,yy,zz,r)\n        found+=1\n        print('truth:', d.instance_number, d.x, d.y)\n        print('predict:', n, x, y, f'inside={inside}', f'array index z={z}')\n    print(f'found={found}')\n    print('')","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:17.048237Z","iopub.execute_input":"2024-07-30T03:46:17.049432Z","iopub.status.idle":"2024-07-30T03:46:19.553572Z","shell.execute_reply.started":"2024-07-30T03:46:17.049375Z","shell.execute_reply":"2024-07-30T03:46:19.552281Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#example of correct cross view 3d to 2d projection\nd = study_id_df.iloc[0] \nprint(d)\nprint('')\n\n\nfor r in mri:\n    d0 = r.df.iloc[0]\n    print(d0.series_id)\n    print(d0.series_description)\n    \n    inside,x,y,z,n = backproject_XYZ(d.xx,d.yy,d.zz,r) \n    print('predict:', n, x, y, f'inside={inside}', f'array index z={z}')\n \n    \n    if inside:\n        slice = r.volume[z].copy()\n        slice[y] = 255 #slice[y] // 2 + 127\n        slice[:,x] = 255 #slice[:,x] // 2 + 127\n    else:\n        D,H,W = r.volume.shape\n        slice = np.zeros((H,W),dtype=np.uint8)\n    plt.imshow(slice, cmap='gray')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:19.555489Z","iopub.execute_input":"2024-07-30T03:46:19.556195Z","iopub.status.idle":"2024-07-30T03:46:20.859954Z","shell.execute_reply.started":"2024-07-30T03:46:19.556150Z","shell.execute_reply":"2024-07-30T03:46:20.858736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#example of wrong cross view 3d to 2d projection ?????\nstudy_id = 4096820034\nstudy_id_df = label_df[label_df.study_id==study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\n#print(study_id_df.iloc[0])\n\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))\nprint('')\n\n\nd = study_id_df.iloc[0] \nprint(d)\nprint('')\n\n\nfor r in mri:\n    d0 = r.df.iloc[0]\n    print(d0.series_id)\n    print(d0.series_description)\n    \n    inside,x,y,z,n = backproject_XYZ(d.xx,d.yy,d.zz,r) \n    print('predict:', n, x, y, f'inside={inside}', f'array index z={z}')\n \n    \n    if inside:\n        slice = r.volume[z].copy()\n        slice[y] = 255 #slice[y] // 2 + 127\n        slice[:,x] = 255 #slice[:,x] // 2 + 127\n    else:\n        D,H,W = r.volume.shape\n        slice = np.zeros((H,W),dtype=np.uint8)\n    plt.imshow(slice, cmap='gray')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:20.861413Z","iopub.execute_input":"2024-07-30T03:46:20.861786Z","iopub.status.idle":"2024-07-30T03:46:27.676480Z","shell.execute_reply.started":"2024-07-30T03:46:20.861754Z","shell.execute_reply":"2024-07-30T03:46:27.675142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#here is the magic !!!\n#plot axial slices\n\ncolor_table = [\n    [1,0,0],\n    [0,1,0],\n    [0,0,1],\n    [1,1,0],\n    [0,1,1],\n]\n\nstudy_id = 40745534 #88465004 #40745534 #4003253 #11340341 #4096820034\nstudy_id_df = label_df[label_df.study_id == study_id]\nstudy_id_df = add_XYZ_to_label_df(study_id_df)\nstudy_id_df = study_id_df.sort_values(['series_id','instance_number'])\nmri = load_for_one(study_id_df)\nprint('len(mri): ', len(mri))\nprint('')\n\n\nfig = plt.figure()\nax = fig.add_subplot(projection='3d')\nax.set_aspect('equal')\n\ncolor_i=0\nfor j,r in enumerate(mri):\n    D, H, W = r.volume.shape\n    if not 'Axial' in r.df.iloc[0].series_description: continue\n\n    color = color_table[color_i]\n    color_i += 1\n\n    N = len(r.df)\n    for i in range(N):\n        d0 = r.df.iloc[i]\n\n        o0, o1, o2, o3, o4, o5 = d0.ImageOrientationPatient\n        ox = np.array([o0, o1, o2])\n        oy = np.array([o3, o4, o5])\n        sx,sy,sz = d0.ImagePositionPatient\n        s = np.array([sx,sy,sz])\n        delx, dely = d0.PixelSpacing\n\n        p0 = s\n        p1 = s + W*delx*ox\n        p2 = s + H*dely*oy\n        p3 = s + H*dely*oy + W*delx*ox\n\n        grid = np.stack([p0, p1, p2, p3]).reshape(2,2,3)\n        gx = grid[:,:,0]\n        gy = grid[:, :, 1]\n        gz = grid[:, :, 2]\n\n        if i==0:\n            ax.plot_surface(gx, gy, gz, alpha=0.7, color=color)\n            ax.scatter([sx], [sy], [sz], color='black')\n\n        else:\n            ax.plot_surface(gx, gy, gz, alpha=0.3, color=color)\n            ax.scatter([sx], [sy], [sz], alpha=0.1, color='black')\n\n    #---\n\n    #plt.show()\n\nxLabel = ax.set_xlabel('x')\nyLabel = ax.set_ylabel('y')\nzLabel = ax.set_zlabel('z')\nplt.show()\nzz=0","metadata":{"execution":{"iopub.status.busy":"2024-07-30T03:46:27.678292Z","iopub.execute_input":"2024-07-30T03:46:27.678814Z","iopub.status.idle":"2024-07-30T03:46:31.609003Z","shell.execute_reply.started":"2024-07-30T03:46:27.678778Z","shell.execute_reply":"2024-07-30T03:46:31.607819Z"},"trusted":true},"execution_count":null,"outputs":[]}]}