{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"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"},{"sourceId":182898236,"sourceType":"kernelVersion"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Locate Lumbar Spine","metadata":{}},{"cell_type":"markdown","source":"Thanks to the organizer and the community, you made this awesome learning possible, appreciated very much. <br>This is to share few findings so far during the journey, and share the experiments of finding lumbar spine in MRI images.","metadata":{}},{"cell_type":"markdown","source":"## setup environment","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport sqlite3\nimport cv2\nimport pydicom as dicom\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\ndata_path = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification'","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:05.154655Z","iopub.execute_input":"2024-06-11T22:53:05.155034Z","iopub.status.idle":"2024-06-11T22:53:06.039541Z","shell.execute_reply.started":"2024-06-11T22:53:05.154993Z","shell.execute_reply":"2024-06-11T22:53:06.038346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# In order to run this Notebook on Kaggle, add this ( https://www.kaggle.com/code/jeff271/lumbar-spine-2024 ) as an input\n# run this cell to make \"util\" folder available in working directory\n\n!ln -s /kaggle/input/lumbar-spine-2024/util /kaggle/working/util 2> /dev/null","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:06.041466Z","iopub.execute_input":"2024-06-11T22:53:06.041925Z","iopub.status.idle":"2024-06-11T22:53:07.058891Z","shell.execute_reply.started":"2024-06-11T22:53:06.041894Z","shell.execute_reply":"2024-06-11T22:53:07.057330Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from util.lumbar_spine import gray_norm,get_hist,resize","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:07.060553Z","iopub.execute_input":"2024-06-11T22:53:07.060903Z","iopub.status.idle":"2024-06-11T22:53:09.059272Z","shell.execute_reply.started":"2024-06-11T22:53:07.060871Z","shell.execute_reply":"2024-06-11T22:53:09.058251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Normalize Sagittal T1 \n - deal with noice, low contrast. transform local representations","metadata":{}},{"cell_type":"code","source":"# select 5 Sagittal T1 images, choose instance number in the middle (1-18) say 9.\ndef get_samples(num:int=5):\n    DB  = os.path.join(os.getcwd(),'util','lumbar.db')\n    con = sqlite3.connect(DB)\n    sql = f'SELECT * FROM csv_series where series_description = \"Sagittal T1\" order by random() limit {num}'\n    ds  = pd.read_sql(sql,con)\n    con.close()\n\n    dcms ={}\n    for ix in range(len(ds)):\n        # read image\n        pID,sID = ds.iloc[ix,0],ds.iloc[ix,1]\n        fn = os.path.join(data_path,'train_images',str(pID),str(sID),'9')\n        fn+='.dcm'\n        dcms[ix] = dicom.dcmread(fn)\n    \n    return dcms","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:09.061438Z","iopub.execute_input":"2024-06-11T22:53:09.061757Z","iopub.status.idle":"2024-06-11T22:53:09.069983Z","shell.execute_reply.started":"2024-06-11T22:53:09.061732Z","shell.execute_reply":"2024-06-11T22:53:09.068774Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# show downsized images, observe histograms, compare several approaches.\n\ndcms = get_samples(5)\nimgs = {} # keep processed images\n\nfig,ax = plt.subplots(8,len(dcms),figsize=(13,13),gridspec_kw={'height_ratios':[2,1,2,1,2,1,2,1]})\nfor ix in range(len(dcms)):\n\n    # norm to [0,1], resize to 128\n    img = dcms[ix].pixel_array\n    img = img.astype(np.float32)/img.max()\n    img = resize(img)\n\n    # method 1: equalization\n    im1 = cv2.equalizeHist((img*2**8).astype(np.uint8))\n    \n    # method 2: personal recipe\n    im2 = gray_norm((img*2**16).astype(np.int32),sz:=128//8)\n    \n    # method 3: CLAHE\n    clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(sz, sz))\n    im3 = clahe.apply((img*2**16).astype(np.uint16))\n    \n    # maping to [-1,1]\n    img = img*2-1; im1 = im1/2**8 *2-1; im3 = im3/2**16 *2-1; \n    \n    # visualize\n    for i,mx in enumerate([img,im1,im2,im3]):\n        ax[i*2,ix].imshow(mx,cmap='gray')            \n        idx,pdf,cdf=get_hist(mx)\n        ax[i*2+1,ix].plot(idx,pdf); aw=ax[i*2+1,ix].twinx(); aw.plot(idx,cdf,'--');aw.get_yaxis().set_visible(False)\n        for x in ['top','right','left']:\n            ax[i*2+1,ix].spines[x].set_visible(False)\n            aw.spines[x].set_visible(False)\n\n    # keep for later use\n    imgs[ix] = [img,im1,im2,im3]\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:09.071999Z","iopub.execute_input":"2024-06-11T22:53:09.072781Z","iopub.status.idle":"2024-06-11T22:53:19.265380Z","shell.execute_reply.started":"2024-06-11T22:53:09.072737Z","shell.execute_reply":"2024-06-11T22:53:19.264002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Convolution to Extract Features","metadata":{}},{"cell_type":"code","source":"# get points of interest, POI.\nfrom scipy.signal import convolve2d\ndef conv(im):\n    h,w = im.shape\n\n    # narrow POI area\n    wL,wR = w//3,w*2//3  # left right\n    hU,hD = h//12,h*4//5 # up down\n    img = im[hU:hD,wL:wR]\n\n    # simple kernel for T1..5 : upper 3/4 white, lower 1/4 black\n    kernel = np.ones((sz:=w//12,sz),np.int32)\n    kernel[-sz//4:] = -2 ; \n\n    # like CNN but in a determinism way\n    conv = convolve2d(img,kernel,'same')\n    \n    # sort conv val, get x,y position. \n    pos = np.dstack(np.unravel_index(np.argsort(conv.ravel()), img.shape))\n    pos = pos[0]                   # dim3 [1,:,2] to dim2 [:,2]\n\n    # POI shall be in this portion. \n    num = int(np.sqrt(img.shape[0]*img.shape[0])) \n    pos = pos[-num:]\n\n    # restore to input (im) coordinates\n    pos[:,0] += hU;  pos[:,1] += wL\n    \n    return pos","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:19.266745Z","iopub.execute_input":"2024-06-11T22:53:19.267169Z","iopub.status.idle":"2024-06-11T22:53:19.275617Z","shell.execute_reply.started":"2024-06-11T22:53:19.267139Z","shell.execute_reply":"2024-06-11T22:53:19.274273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# show POI\n\nfig,ax = plt.subplots(4,len(dcms),figsize=(13,7))\nfor ix in range(len(dcms)):\n    # rows: original img + processed images with method 1 2 3; cols: samples\n    for i,mx in enumerate(imgs[ix]):\n        # plot poi on original image\n        ax[i,ix].imshow(imgs[ix][0],cmap='gray')            \n        \n        poi = conv(mx)\n        for x,y in poi:\n            ax[i,ix].scatter(y,x)\n\n        # keep POI for clusterring later on\n        imgs[f'poi_{ix}_{i}'] = poi\n        \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:19.276813Z","iopub.execute_input":"2024-06-11T22:53:19.277671Z","iopub.status.idle":"2024-06-11T22:53:38.297686Z","shell.execute_reply.started":"2024-06-11T22:53:19.277625Z","shell.execute_reply":"2024-06-11T22:53:38.296222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Cluster POI","metadata":{}},{"cell_type":"code","source":"!pip install hdbscan","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:53:38.299297Z","iopub.execute_input":"2024-06-11T22:53:38.299761Z","iopub.status.idle":"2024-06-11T22:54:08.033758Z","shell.execute_reply.started":"2024-06-11T22:53:38.299720Z","shell.execute_reply":"2024-06-11T22:54:08.032473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# choose DBSCAN since it doesn't require pre-defined num of clusters.\n# from sklearn.cluster import HDBSCAN\nfrom hdbscan import HDBSCAN\n# cluster_selection_method, not sure eom or leaf which fits better? choose default for now.\n# choose medoid over centroid since it guarantees the output is a POI data point.\nhdb = HDBSCAN(min_cluster_size=2,cluster_selection_epsilon = 128//16,cluster_selection_method = 'eom')\n\nfig,ax = plt.subplots(4,len(dcms),figsize=(13,7))\nfor ix in range(len(dcms)):\n    # rows: original img + processed images with 3 methods; cols: samples\n    for i,mx in enumerate(imgs[ix]):\n        # plot clusters on original image\n        ax[i,ix].imshow(imgs[ix][0],cmap='gray')            \n        \n        poi = imgs[f'poi_{ix}_{i}']\n        model = hdb.fit(poi)\n        medoids = [model.weighted_cluster_medoid(i) for i in np.unique(model.labels_) if i>=0]\n        \n        for x,y in medoids:\n            ax[i,ix].scatter(y,x)\n        \nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:55:01.306320Z","iopub.execute_input":"2024-06-11T22:55:01.306721Z","iopub.status.idle":"2024-06-11T22:55:04.833637Z","shell.execute_reply.started":"2024-06-11T22:55:01.306690Z","shell.execute_reply":"2024-06-11T22:55:04.832294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Locate Lumbar Spine\nwith these medoids, need more work to map to T1..5,S1. seems possible ...","metadata":{}},{"cell_type":"markdown","source":"## Appendix - DB init","metadata":{}},{"cell_type":"code","source":"# one time job. could access at /kaggle/working/util/lumbar.db\ndef db_init():\n    # dcm files and meta to db\n    con = sqlite3.connect(DB)\n    cur = con.cursor()\n    \n    cur.execute('drop table if exists imgs_meta')\n    cur.execute('create table imgs_meta (study int,series int,dcm int,PatientID int,SeriesID int,InstanceID int,PosX real,PosY real,PoxZ real);')\n    cur.execute('BEGIN TRANSACTION;')\n    for rootdir, dirs, files in os.walk(os.path.join(data_path,'train_images')):\n        for file in files:\n            if file.endswith('.dcm'):\n                study,series=rootdir.split(os.sep)[-2:]\n                dcm = dicom.dcmread(rootdir+os.sep+file)\n                IDp,IDs,IDi,pos = [dcm[c].value for c in 'PatientID,SeriesInstanceUID,InstanceNumber,ImagePositionPatient'.split(',')]\n                _ = ','.join(['?']*9)\n                cur.execute(f'insert into imgs_meta (study,series,dcm,PatientID,SeriesID,InstanceID,PosX,PosY,PoxZ) values({_})',\\\n                            (study,series,file[:-4],IDp,IDs.split('.')[1],IDi,*pos))\n    cur.execute('COMMIT;')        \n\n    # transform csv to db\n    ds = pd.read_csv(os.path.join(data_path,'train.csv'))\n    df = ds.melt(id_vars=['study_id'])\n    \n    levels = df['value'].unique()\n    levels = dict(zip(range(len(levels)),levels))\n    vk = dict([(v,k) for k,v in levels.items()])\n    df['value'] = df['value'].apply(lambda x:vk[x])\n    \n    locations = df['variable'].unique()\n    locations = dict(zip(range(len(locations)),locations))\n    vk = dict([(v,k) for k,v in locations.items()])\n    df['variable'] = df['variable'].apply(lambda x:vk[x])    \n    df.to_sql('csv_train',con,if_exists='replace',index=False)\n    \n    cur.execute('drop table if exists map_levels')\n    cur.execute('drop table if exists map_locations')\n    cur.execute('create table map_levels (k int, v text)')\n    cur.execute('create table map_locations (k int, v text)')\n    for k,v in levels.items():\n        cur.execute('insert into map_levels (k,v) values (?,?)',(k,v))\n    for k,v in locations.items():\n        cur.execute('insert into map_locations (k,v) values (?,?)',(k,v))\n    con.commit()\n    \n    ds = pd.read_csv(os.path.join(data_path,'train_series_descriptions.csv'))\n    ds.to_sql('csv_series',con,if_exists='replace',index=False)\n    \n    ds = pd.read_csv(os.path.join(data_path,'train_label_coordinates.csv'))\n    vk = dict([(v,k) for k,v in locations.items()])\n    def f(x):\n        a,b = x\n        a = a.lower().replace(' ','_')\n        b = b.lower().replace(r'/','_')\n        return vk[a+'_'+b]\n    ds['location'] = ds[['condition','level']].apply(f,axis=1)\n    ds.to_sql('csv_label',con,if_exists='replace',index=False)\n    \n    con.close()    ","metadata":{"execution":{"iopub.status.busy":"2024-06-11T22:54:12.375120Z","iopub.status.idle":"2024-06-11T22:54:12.375584Z","shell.execute_reply.started":"2024-06-11T22:54:12.375386Z","shell.execute_reply":"2024-06-11T22:54:12.375406Z"},"trusted":true},"execution_count":null,"outputs":[]}]}