{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# SIIM COVID-19: Process to JPG 512 with bboxes\nImage processing includes:\n- reading DICOM, correcting BW and VOI_LUT\n- removing empty borders, \n- square crop at the center attempting to keep boxes uncropped, \n- resizing down to 512x512 pixels size, \n- histogram equalization, \n- transforming to 8-bits,\n- saving images to JPG using a good quality=95.\n\nAnnotation processing: \n- recalculating bboxes to match new crop and size, \n- add human-readable labels from study level,\n- add some DICOM fields.\n\nAdditional annotation file contains only samples with the unique StudyInstanceUID. I would suggest to keep only one image per study, because the rest of the images are variation of the same image with increased sharpness or other augmentations.\n\nAnnotation file 'test_hwtblr.csv' is for the test images. Column 'hwtblr' contains a list of values for the: original_height, original_width, crop_top_line, crop_bottom_line, crop_left_pixel, crop_right_pixel. This data will help to re-calculate bboxes to the original image size.\n\n### Upvote, thanks!","metadata":{}},{"cell_type":"markdown","source":"# Acknowlegdements\nI have started this notebook from this: https://www.kaggle.com/code/xhlulu/siim-covid-19-convert-to-jpg-256px, but later stepped away to another direction.\n\nInspiration for image processing came from: https://www.kaggle.com/code/davidbroberts/export-processed-jpg-512/notebook and https://www.kaggle.com/code/raddar/convert-dicom-to-np-array-the-correct-way","metadata":{}},{"cell_type":"markdown","source":"## TODO\n- File ID '06f6423be3f9' is corrupted and can be removed\n- can create a list of useless IDs and exclude those from the final annotation\n- border crop can be smarter than simply check for equal pixels, some borders contain garbage","metadata":{}},{"cell_type":"code","source":"!pip install python-gdcm -q","metadata":{"execution":{"iopub.status.busy":"2022-11-06T16:58:11.657290Z","iopub.execute_input":"2022-11-06T16:58:11.657596Z","iopub.status.idle":"2022-11-06T16:58:20.555408Z","shell.execute_reply.started":"2022-11-06T16:58:11.657568Z","shell.execute_reply":"2022-11-06T16:58:20.554445Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from fastai.vision.all import *\nfrom fastai.medical.imaging import *\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport pandas as pd\nimport skimage\nfrom skimage import transform, exposure, io","metadata":{"execution":{"iopub.status.busy":"2022-11-06T16:58:20.557075Z","iopub.execute_input":"2022-11-06T16:58:20.557397Z","iopub.status.idle":"2022-11-06T16:58:23.835212Z","shell.execute_reply.started":"2022-11-06T16:58:20.557355Z","shell.execute_reply":"2022-11-06T16:58:23.834379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_items = get_dicom_files('../input/siim-covid19-detection/train')\nprint(f'Got {len(train_items)} DICOM files from the TRAIN path')\ntest_items = get_dicom_files('../input/siim-covid19-detection/test')\nprint(f'Got {len(test_items)} DICOM files from the TEST path')","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:21.506077Z","iopub.execute_input":"2022-11-06T17:32:21.506423Z","iopub.status.idle":"2022-11-06T17:32:28.967914Z","shell.execute_reply.started":"2022-11-06T17:32:21.506390Z","shell.execute_reply":"2022-11-06T17:32:28.966845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_items = train_items[:12] # for test\n# test_items = train_items[:12] # for test","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:31.579535Z","iopub.execute_input":"2022-11-06T17:32:31.579994Z","iopub.status.idle":"2022-11-06T17:32:31.589891Z","shell.execute_reply.started":"2022-11-06T17:32:31.579955Z","shell.execute_reply":"2022-11-06T17:32:31.589297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preprocess annotations","metadata":{}},{"cell_type":"code","source":"# read image and study annotations\nsdf = pd.read_csv('../input/siim-covid19-detection/train_study_level.csv')\nidf = pd.read_csv('../input/siim-covid19-detection/train_image_level.csv')","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:34.055549Z","iopub.execute_input":"2022-11-06T17:32:34.056049Z","iopub.status.idle":"2022-11-06T17:32:34.082501Z","shell.execute_reply.started":"2022-11-06T17:32:34.055990Z","shell.execute_reply":"2022-11-06T17:32:34.081651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create human-readable classes\ndef lab_y( r ):\n    if r['Negative for Pneumonia']: return 'Neg'\n    if r['Typical Appearance']: return 'Typ'\n    if r['Atypical Appearance']: return 'Atyp'\n    if r['Indeterminate Appearance']: return 'Indet'\n\nsdf['label_y'] = sdf.apply(lab_y, axis=1)\nsdf.id = sdf.apply(lambda r: r.id.replace('_study',''), axis=1)\nsdf = sdf.set_index('id')","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:37.221568Z","iopub.execute_input":"2022-11-06T17:32:37.221899Z","iopub.status.idle":"2022-11-06T17:32:37.383293Z","shell.execute_reply.started":"2022-11-06T17:32:37.221866Z","shell.execute_reply":"2022-11-06T17:32:37.382709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idf.drop(columns=['label'],inplace=True) # irrelevant column\n# add text labels from study csv\nidf['label_y'] = idf.apply(lambda r: sdf.loc[r.StudyInstanceUID].label_y, axis=1)\nidf.id = idf.apply(lambda r: r.id.replace('_image',''), axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:43.304969Z","iopub.execute_input":"2022-11-06T17:32:43.305489Z","iopub.status.idle":"2022-11-06T17:32:44.410869Z","shell.execute_reply.started":"2022-11-06T17:32:43.305454Z","shell.execute_reply":"2022-11-06T17:32:44.410193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# calculate a box, bounding all boxes\n# it will help to avoid/minimise crop of boxes\ndef box_all(r): # input: boxes from dataframe\n    boxes = r.boxes\n    if type(boxes) == list:\n        boxl = boxes\n        w = 'w'\n        h = 'h'\n    else: \n        if pd.isnull(boxes): return boxes\n        boxl = eval(boxes)\n        w = 'width'\n        h = 'height'\n    minx = min([b['x'] for b in boxl])\n    miny = min([b['y'] for b in boxl])\n    maxx = max([b['x']+b[w] for b in boxl])\n    maxy = max([b['y']+b[h] for b in boxl])\n    return [minx, miny, maxx, maxy]\n\nidf['bb_min_xy_max_xy'] = idf.apply(box_all, axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:47.683762Z","iopub.execute_input":"2022-11-06T17:32:47.684213Z","iopub.status.idle":"2022-11-06T17:32:47.887948Z","shell.execute_reply.started":"2022-11-06T17:32:47.684181Z","shell.execute_reply":"2022-11-06T17:32:47.887113Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# keep copy with the original bboxes for visualization check below\nidf_copy = idf.copy()","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:52.925045Z","iopub.execute_input":"2022-11-06T17:32:52.925381Z","iopub.status.idle":"2022-11-06T17:32:52.930844Z","shell.execute_reply.started":"2022-11-06T17:32:52.925348Z","shell.execute_reply":"2022-11-06T17:32:52.929422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Process and convert DICOM images","metadata":{}},{"cell_type":"code","source":"# size of the output image\nDSIZE = 512 ","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:57.154813Z","iopub.execute_input":"2022-11-06T17:32:57.155160Z","iopub.status.idle":"2022-11-06T17:32:57.159313Z","shell.execute_reply.started":"2022-11-06T17:32:57.155129Z","shell.execute_reply":"2022-11-06T17:32:57.158302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Functions to remove borders and do maximum center crop, trying to keep boxes\ndef topbot_borders(pixels):\n    t = 0\n    b = 0\n    h = pixels.shape[0]\n    for i in range(h):\n        if not np.all(pixels[i] == pixels[i][0]):\n            t = i\n            break\n              \n    for i in range(h-1, 0, -1):\n        if not np.all(pixels[i] == pixels[i][0]):\n            b = i\n            break\n    return t, b\n\ndef remove_borders(pixels, bb_min_xy_max_xy):\n    h = pixels.shape[0]\n    w = pixels.shape[1]\n    t,b = topbot_borders(pixels)\n    pixels = pixels[t:b,:] \n    pixels = np.rot90(pixels)\n    r,l = topbot_borders(pixels)\n    pixels = pixels[r:l,:]\n    pixels = np.rot90(pixels, 3)\n    l = w-l\n    r = w-r\n    # crop square shape\n    nh = pixels.shape[0]\n    nw = pixels.shape[1]\n    d = abs(nh - nw)\n    st = d // 2\n    sp = d-st\n    if nh == nw: pass\n    elif nh > nw:\n        # check to include boxes\n        if bb_min_xy_max_xy is not None:\n            if bb_min_xy_max_xy[0] < st:\n                st = bb_min_xy_max_xy[0]\n                sp = d-st\n            elif bb_min_xy_max_xy[2] > nw-sp:\n                sp = nw - bb_min_xy_max_xy[2]\n                st = d - sp\n        pixels = pixels[st:-sp,:]\n        t += st\n        b -= sp\n    elif nw > nh:\n        # check to include boxes\n        if bb_min_xy_max_xy is not None:\n            if bb_min_xy_max_xy[1] < st:\n                st = bb_min_xy_max_xy[1]\n                sp = d-st\n            elif bb_min_xy_max_xy[3] > nh-sp:\n                sp = nh - bb_min_xy_max_xy[3]\n                st = d - sp\n        pixels = pixels[:,st:-sp]\n        l += st\n        r -= sp\n    return pixels, [h,w,t,b,l,r]","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:32:59.514983Z","iopub.execute_input":"2022-11-06T17:32:59.515321Z","iopub.status.idle":"2022-11-06T17:32:59.526733Z","shell.execute_reply.started":"2022-11-06T17:32:59.515288Z","shell.execute_reply":"2022-11-06T17:32:59.525502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# process DICOM files\ndef proc_dcm(items, folder, bdf, voi_lut = True):\n\n    proc_df = pd.DataFrame(columns =['id','hwtblr', 'Modality', 'ImagerPixelSpacing', 'StudyDate', 'PatientSex', 'BodyPartExamined'])\n\n    for ti in items:\n        id = ti.stem\n        img = ti.dcmread()\n        if voi_lut:\n            pixels = apply_voi_lut(img.pixel_array, img)\n        else:\n            pixels = img.pixel_array\n        #pixels = img.pixel_array\n        bb_mimax = None\n        if bdf is not None:\n            if id in bdf.index:\n                bb_mimax = bdf.loc[id].bb_min_xy_max_xy\n        pixels,hwtblr = remove_borders(pixels, bb_mimax)\n        if img.PhotometricInterpretation == \"MONOCHROME1\":\n            pixels = np.amax(pixels) - pixels\n        pixels = skimage.transform.resize(pixels, (DSIZE, DSIZE))\n        pixels = skimage.exposure.equalize_hist(pixels,nbins=256)\n        #pixels = skimage.filters.unsharp_mask(pixels, radius=5, amount=1)\n        pixels = (pixels * 255).astype(np.uint8)\n        # fig, ax = plt.subplots(figsize=(12,12))\n        # im = ax.imshow(pixels, cmap='gray')\n        # plt.show()\n\n        skimage.io.imsave(folder + '/' + id + '.jpg', pixels, quality=95)\n        proc_df = proc_df.append({'id': id, 'hwtblr': hwtblr, 'Modality': img.Modality,\n            'ImagerPixelSpacing': img.ImagerPixelSpacing, 'StudyDate': img.StudyDate, \n            'PatientSex': img.PatientSex, 'BodyPartExamined': img.BodyPartExamined}, ignore_index=True)\n    return proc_df","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:31.004317Z","iopub.execute_input":"2022-11-06T17:33:31.004799Z","iopub.status.idle":"2022-11-06T17:33:31.012456Z","shell.execute_reply.started":"2022-11-06T17:33:31.004752Z","shell.execute_reply":"2022-11-06T17:33:31.011500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Processing train images')\n!mkdir train\ntrain_proc_df = proc_dcm(train_items, 'train', idf )","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:35.536551Z","iopub.execute_input":"2022-11-06T17:33:35.536867Z","iopub.status.idle":"2022-11-06T17:33:38.465520Z","shell.execute_reply.started":"2022-11-06T17:33:35.536837Z","shell.execute_reply":"2022-11-06T17:33:38.464805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Processing test images')\n!mkdir test\ntest_proc_df = proc_dcm(test_items, 'test', None )","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:38.735016Z","iopub.execute_input":"2022-11-06T17:33:38.735515Z","iopub.status.idle":"2022-11-06T17:33:41.462974Z","shell.execute_reply.started":"2022-11-06T17:33:38.735476Z","shell.execute_reply":"2022-11-06T17:33:41.462128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Process annotation files\nAnnotation processing: recalculating bboxes to match new crop and size, add human-readable labels from study level.","metadata":{}},{"cell_type":"code","source":"# save and index train image crop data\ntrain_proc_df.to_csv('train_proc.csv',index=False) # don't need to save this file, for debug purpose\ntrain_proc_df = train_proc_df.set_index('id')\ntrain_proc_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:46.539637Z","iopub.execute_input":"2022-11-06T17:33:46.539962Z","iopub.status.idle":"2022-11-06T17:33:46.560173Z","shell.execute_reply.started":"2022-11-06T17:33:46.539925Z","shell.execute_reply":"2022-11-06T17:33:46.559268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# save dataframe with test image crop data\ntest_proc_df.to_csv('test_proc.csv',index=False) \ntest_proc_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:50.378912Z","iopub.execute_input":"2022-11-06T17:33:50.379546Z","iopub.status.idle":"2022-11-06T17:33:50.397883Z","shell.execute_reply.started":"2022-11-06T17:33:50.379500Z","shell.execute_reply":"2022-11-06T17:33:50.396948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Recalculate boxes for the new crop and new size","metadata":{}},{"cell_type":"code","source":"def calc_boxes(o, bdf):\n    # o is row of the train annotation file\n    # bdf is the dataframe of crop results\n    if pd.isnull(o.boxes): return\n    if o.id not in bdf.index: return o.boxes\n    boxl = eval(o.boxes)\n    hwtblr = bdf.loc[o.id].hwtblr\n    new_boxl = []\n    for bb in boxl:\n        cw = hwtblr[5] - hwtblr[4] # width after crop\n        ch = hwtblr[3] - hwtblr[2] # height after crop\n        \n        bbx = bb['x'] - hwtblr[4] # left crop\n        if bbx < 0: bbx = 0\n        bbx = round((bbx * DSIZE) / cw, 2)\n            \n        bby = bb['y'] - hwtblr[2] # top crop\n        if bby < 0: bby = 0\n        bby = round((bby * DSIZE) / ch, 2)\n            \n        bbw = round((bb['width'] * DSIZE) / cw, 2)\n        if bbx+bbw > DSIZE: bbw = DSIZE-bbx\n        \n        bbh = round((bb['height'] * DSIZE) / ch, 2)\n        if bby+bbh > DSIZE: bbh = DSIZE-bby\n        \n        bbd = {'x':bbx, 'y':bby, 'w':bbw, 'h':bbh}\n        new_boxl.append(bbd)\n    #print(o.id,new_boxl)\n    return new_boxl ","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:33:56.535274Z","iopub.execute_input":"2022-11-06T17:33:56.536250Z","iopub.status.idle":"2022-11-06T17:33:56.544447Z","shell.execute_reply.started":"2022-11-06T17:33:56.536206Z","shell.execute_reply":"2022-11-06T17:33:56.543389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# replace boxes with updated values\nidf.boxes = idf.apply(lambda r: calc_boxes(r, train_proc_df), axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:34:02.572202Z","iopub.execute_input":"2022-11-06T17:34:02.572516Z","iopub.status.idle":"2022-11-06T17:34:02.733561Z","shell.execute_reply.started":"2022-11-06T17:34:02.572485Z","shell.execute_reply":"2022-11-06T17:34:02.732432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idf = idf.set_index('id')\nidf = idf.join(train_proc_df)\n\n# update boxes outline with new box values\nidf['bb_min_xy_max_xy'] = idf.apply(box_all, axis=1)\n\n# save processed CSV for future use\nidf.to_csv('train_512.csv',index=True)\nidf.head()","metadata":{"execution":{"iopub.status.busy":"2022-11-06T17:36:45.565315Z","iopub.execute_input":"2022-11-06T17:36:45.565621Z","iopub.status.idle":"2022-11-06T17:36:45.590047Z","shell.execute_reply.started":"2022-11-06T17:36:45.565593Z","shell.execute_reply":"2022-11-06T17:36:45.588639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Output train annotation file","metadata":{}},{"cell_type":"code","source":"udf = idf.drop_duplicates(subset=['StudyInstanceUID'], keep='first')\nudf.to_csv('train_512_uniqueStudy.csv',index=True)\nprint(udf.shape)","metadata":{"execution":{"iopub.status.busy":"2022-10-30T14:03:19.654934Z","iopub.execute_input":"2022-10-30T14:03:19.655327Z","iopub.status.idle":"2022-10-30T14:03:19.689628Z","shell.execute_reply.started":"2022-10-30T14:03:19.655291Z","shell.execute_reply":"2022-10-30T14:03:19.688025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Have a nice day!\nPlease, upvote, thanks!","metadata":{}}]}