{"cells":[{"metadata":{},"cell_type":"markdown","source":"# RSNA PE Detection Submission\nRobert Turley\nOctober 24, 2020\n\n"},{"metadata":{},"cell_type":"markdown","source":"### Libraries and Packages"},{"metadata":{"trusted":true},"cell_type":"code","source":"import time\nstart_time = time.time()\nTIME_LIMIT = 60*60*8.25  # 8 and one quarter hours","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install ../input/fastai2013py3/fastcore-1.0.12-py3-none-any.whl\n!pip install ../input/fastai2013py3/torch-1.6.0cu101-cp37-cp37m-linux_x86_64.whl\n!pip install ../input/fastai2013py3/torchvision-0.7.0cu101-cp37-cp37m-linux_x86_64.whl\n!pip install ../input/fastai2013py3/fastai-2.0.13-py3-none-any.whl","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!cp ../input/gdcm-conda-install/gdcm.tar .\n!tar -xvzf gdcm.tar\n!conda install --offline ./gdcm/gdcm-2.8.9-py37h71b2a6d_0.tar.bz2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import pydicom\nfrom os import path\nfrom fastai.vision.all import *\nfrom scipy import signal\nimport dask\nimport lightgbm as lgb\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Confirm Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"data_path = Path('/kaggle/input/rsna-str-pulmonary-embolism-detection/')\nmodel_path = Path('/kaggle/input/pe-models')\n\nMIN_LUNG_FRACTION = 0.05\nDO_PARTIAL = False","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nif path.exists('../input/rsna-str-pulmonary-embolism-detection/train'):\n    test_df = pd.read_csv(data_path/'test.csv').head(400)\nelse:\n    test_df = pd.read_csv(data_path/'test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"############################################\n#### MODELS\n############################################\n\nventricle_model = load_learner(model_path/'ventricle-model.pkl',cpu=False)\nventclass_model = load_learner(model_path/'ventclass-model.pkl',cpu=False)\ncontrastclass_model = load_learner(model_path/'contrastclass-model.pkl',cpu=False)\nmotionclass_model =  load_learner(model_path/'motionclass-model.pkl',cpu=False)\npeclass_model = load_learner(model_path/'peclass_model.pkl',cpu=False)\n\n\ngbm_outcome = []\ngbm_vent = []\ngbm_left = []\ngbm_center = []\ngbm_right = []\ngbm_poi = []\n\nfor fold in range(5):\n    gbm_outcome.append( lgb.Booster(model_file=str(model_path/('gbm_outcome_f'+str(fold)+'.txt'))))\n    gbm_vent.append( lgb.Booster(model_file=str(model_path/('gbm_vent_f'+str(fold)+'.txt'))))\n    gbm_left.append( lgb.Booster(model_file=str(model_path/('gbm_left_f'+str(fold)+'.txt'))))\n    gbm_center.append( lgb.Booster(model_file=str(model_path/('gbm_center_f'+str(fold)+'.txt'))))\n    gbm_right.append( lgb.Booster(model_file=str(model_path/('gbm_right_f'+str(fold)+'.txt'))))\n    gbm_poi.append( lgb.Booster(model_file=str(model_path/('gbm_poi_f'+str(fold)+'.txt'))))\n\n    \n\n###################################################################\n# Feature Variables\nfeature_cols = [   #DCIM Metadata\n                   'TableHeight', \n                   'PixelSpacing', \n                   'n_file_images',\n                   'n_crop_row', \n                   'n_crop_col',\n                   'body_crop_width',\n                   'body_crop_height',\n                   'lung_fraction_max',\n                   'lung_length_cm',\n                   'lung_area_max_cm2',\n                   'lung_volume_cm3',\n                   # ventricle visibility\n                   'ventricle_visibility_max',\n                   'ventricle_visibility_num_above75',\n                   # right side vs left side prediction\n                   'rvlv_gte1_classification_top',\n                   'rvlv_gte1_classification_min',\n                   'rvlv_gte1_classification_max',\n                   'rvlv_gte1_classification_mean',\n                   # motion prediction\n                   'motion_pred_top',\n                   'motion_pred_min',\n                   'motion_pred_max',\n                   'motion_pred_mean',\n                   # contrast prediction\n                   'contrast_pred_top',\n                   'contrast_pred_min',\n                   'contrast_pred_max',\n                   'contrast_pred_mean',\n                   # PE classification model predictions\n                   'pe_acute_L_max',\n                   'pe_acute_L_above20',\n                   'pe_acute_L_above50',\n                   'pe_acute_L_lungmean',\n                   'pe_acute_C_max',\n                   'pe_acute_C_above20',\n                   'pe_acute_C_above50',\n                   'pe_acute_C_lungmean',\n                   'pe_acute_R_max',\n                   'pe_acute_R_above20',\n                   'pe_acute_R_above50',\n                   'pe_acute_R_lungmean',\n                   'pe_chronic_L_max',\n                   'pe_chronic_L_above20',\n                   'pe_chronic_L_above50',\n                   'pe_chronic_L_lungmean',\n                   'pe_chronic_C_max',\n                   'pe_chronic_C_above20',\n                   'pe_chronic_C_above50',\n                   'pe_chronic_C_lungmean',\n                   'pe_chronic_R_max',\n                   'pe_chronic_R_above20',\n                   'pe_chronic_R_above50',\n                   'pe_chronic_R_lungmean',\n                   'pe_neg_L_max',\n                   'pe_neg_L_above20',\n                   'pe_neg_L_above50',\n                   'pe_neg_L_lungmean',\n                   'pe_neg_C_max',\n                   'pe_neg_C_above20',\n                   'pe_neg_C_above50',\n                   'pe_neg_C_lungmean',\n                   'pe_neg_R_max',\n                   'pe_neg_R_above20',\n                   'pe_neg_R_above50',\n                   'pe_neg_R_lungmean']\n\n\ndirect_image_features = ['image_order', 'lung_fraction', 'ventricle_pred', 'top_vent_image',\n       'top_lung_image', 'pe_acute_L', 'pe_acute_C', 'pe_acute_R',\n       'pe_chronic_L', 'pe_chronic_C', 'pe_chronic_R', 'pe_neg_L', 'pe_neg_C',\n       'pe_neg_R', 'rvlv_gte1_pred', 'motion_pred', 'contrast_pred',\n       'lung_order_pct', 'lung_area_cm2',\n       'lung_area_vs_max']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Code"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nseries_metadata_fields = ['TableHeight','KVP','SliceThickness','Rows','Columns','PixelSpacing',\n                                    'HighBit','BitsStored','BitsAllocated','WindowWidth',\n                                    'WindowCenter','RescaleSlope','RescaleIntercept']\n\nseries_derived_fields = ['n_file_images','n_all_images','n_crop_row','n_crop_col',\n                         'WindowWidthDim','WindowWidth2','WindowCenter2','PixelSpacing2']\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# RGB values to emphasize soft tissue and contrast solution\nall_xp = np.array([50,  51, 75, 100,200,300,600,601])\nr_fm = np.array(  [127,255,255,   0,  0,  0,  0,127])\ng_fm = np.array(  [127,  0,255, 255,255,  0,  0,127])\nb_fm = np.array(  [127,  0,  0,   0,255,255,255,127])\n\nall_xp2 = np.array([-876,-875,-500,-499])\nr_fm2 = np.array(  [ 0,   255, 255,   0])\ng_fm2 = np.array(  [ 0,   230, 200,   0])\nb_fm2 = np.array(  [ 0,   230, 200,   0])\n\n\ndef artery_base_image(raw_pixel3d,all_xp=all_xp,r_fm=r_fm,g_fm=g_fm,b_fm=b_fm):\n    artery_base = np.array([np.interp(raw_pixel3d,all_xp,r_fm),\n                                   np.interp(raw_pixel3d,all_xp,g_fm),\n                                   np.interp(raw_pixel3d,all_xp,b_fm)],dtype='uint8')\n    return np.moveaxis(artery_base,0,-1)\n    \n\n    \n    \n# customized CT bins that map CT image to image data\nctbins =torch.FloatTensor([bval for bval in range(-1000,-900,100)] +     # air and background\n                            [bval for bval in range(-900,-800,25)] +\n                            [bval for bval in range(-800,-750,10)] +\n                            [bval for bval in range(-750,-700,5)] +\n                            [bval for bval in range(-700,-600,5)] +       # focus on the lungs\n                            [bval for bval in range(-600,-500,10)] +\n                            [bval for bval in range(-500,-400,25)] +\n                            [bval for bval in range(-400,-100,50)] + \n                            [bval for bval in range(-100,-20,5)] + \n                            [bval for bval in range(-20,-10,1)] + \n                            [bval for bval in range(-10,90,1)] +          # focus on the water, blood, vessels, and soft tissues\n                            [bval for bval in range(90,150,2)] + \n                            [bval for bval in range(150,200,5)] + \n                            [bval for bval in range(200,400,10)] +         # soft bone tissue\n                            [bval for bval in range(400,600,50)] +\n                            [bval for bval in range(600,1100,100)] )\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\n\ndef rounded_borders(rgb_array):\n    out_array = rgb_array.copy()\n    hght,wdth,_ = rgb_array.shape\n    x,y = np.indices((hght,wdth))\n    border_indices= ((1/(((hght-2)/2)**2))*np.square(x-(hght-1)/2) + (1/(((wdth-2)/2)**2))*np.square(y-(wdth-1)/2) > 1)\n    out_array[border_indices,:] = 127\n    return out_array\n\ndef rounded_black_borders(rgb_array):\n    out_array = rgb_array.copy()\n    hght,wdth,_ = rgb_array.shape\n    x,y = np.indices((hght,wdth))\n    border_indices= ((1/(((hght-40)/2)**2))*np.square(x-(hght-1)/2) + (1/(((wdth-40)/2)**2))*np.square(y-(wdth-1)/2) > 1)\n    out_array[border_indices,:] = 0\n    return out_array\n\ndef variation_crop(raw_pixel3d,other3d, threshold=100):\n    \n    ###########################################################################   \n    #crop areas of predominant black, which helps to center area of interest\n    n_images,n_rows,n_cols = raw_pixel3d.shape\n    if n_rows>500 and n_cols>500:\n        min_variation = 100\n    else:\n        min_variation = 0\n\n    horizontal_variation = np.mean(np.std(raw_pixel3d,axis=0),axis=0)\n    no_horizontals = np.where(horizontal_variation<min_variation)[0]\n\n    if len(no_horizontals)>0:\n        if np.min(no_horizontals)<42:\n            first_a2 = np.max(no_horizontals[no_horizontals<42])\n        else:\n            first_a2 = 0\n\n        if np.max(no_horizontals)>(raw_pixel3d.shape[2]-42):\n            last_a2 = np.min(no_horizontals[no_horizontals>(raw_pixel3d.shape[2]-42)])\n        else:\n            last_a2 = raw_pixel3d.shape[2]\n    else:\n        first_a2 = 0\n        last_a2 = raw_pixel3d.shape[2]\n\n    vertical_variation = np.mean(np.std(raw_pixel3d,axis=0),axis=1)\n    no_verticals = np.where(vertical_variation<min_variation)[0]\n    if len(no_verticals)>0:\n        if np.min(no_verticals)<106:\n            first_a1 = np.max(no_verticals[no_verticals<106])\n        else:\n            first_a1 = 0\n\n        if np.max(no_verticals)>(raw_pixel3d.shape[1]-106):\n            last_a1 = np.min(no_verticals[no_verticals>(raw_pixel3d.shape[1]-106)])\n        else:\n            last_a1 = raw_pixel3d.shape[1]\n    else:\n        first_a1 = 0\n        last_a1 = raw_pixel3d.shape[1]\n\n    return raw_pixel3d[:,first_a1:last_a1,first_a2:last_a2],other3d[:,first_a1:last_a1,first_a2:last_a2]\n\n\n# function create image weight vector\ndef image_wt_vector(image_num,im_present_tf):\n    pre_vector = np.zeros(len(im_present_tf),dtype='float')\n    prewts = np.array([1/64,1/64,1/32,1/16,1/8,1/4])\n    pre_vector[max(0,image_num-5):(image_num+1)] = prewts[max(0,5-image_num):]\n    pre_vector[im_present_tf==False] = 0\n    pre_vector = pre_vector*(0.5/np.sum(pre_vector))\n\n\n\n    post_vector = np.zeros(len(im_present_tf),dtype='float')\n    postwts = np.flip(prewts)\n    post_vector[image_num:min(len(im_present_tf),image_num+6)] = postwts[0:min(7,(len(im_present_tf)-image_num))]\n    post_vector[im_present_tf==False] = 0\n    post_vector = post_vector*(0.5/np.sum(post_vector))\n    return (pre_vector+post_vector)\n\n\ndef lung_crop(rgb_pixels,lung_pixels, threshold = 0.01):\n    \n    good_crop = True\n    \n    sh = rgb_pixels.shape\n    vsize=sh[0]\n    hsize=sh[1]\n    \n    hmeans = np.mean(lung_pixels,axis=0)\n    vmeans = np.mean(lung_pixels,axis=1)\n    \n    # left crop\n    if not good_crop or (np.max(hmeans[0:hsize//4])<threshold):\n        good_crop = False\n        left_edge = hsize//4\n    else:\n        # calculate the left edge\n        if len(np.where(hmeans[0:hsize//4]<threshold)[0])>0:\n            left_edge = np.max(np.where(hmeans[0:hsize//4]<threshold)[0])\n            if left_edge>0:\n                left_edge -= 1\n        else:\n            left_edge = 0\n\n    # right crop        \n    if not good_crop or (np.max(hmeans[3*hsize//4:])<threshold):\n        good_crop = False\n        right_edge = 3*hsize//4\n    else:   \n        # calculate right edge\n        if len(np.where(hmeans[3*hsize//4:]<threshold)[0])>0:\n            right_edge = np.min(3*hsize//4+np.where(hmeans[3*hsize//4:]<threshold)[0])\n            if right_edge<(hsize-1):\n                left_edge += 1\n        else:\n            right_edge = hsize-1\n\n    # top crop\n    if good_crop and (np.max(vmeans[0:vsize//3])<threshold):\n        good_crop = False\n        top_edge = 0\n    else:\n        # calculate the left edge\n        if len(np.where(vmeans[0:vsize//3]<threshold)[0])>0:\n            top_edge = np.max(np.where(vmeans[0:vsize//3]<threshold)[0])\n            if top_edge>0:\n                top_edge -= 1\n        else: \n            top_edge = 0\n            \n    # bottom crop    \n    if not good_crop or (np.max(vmeans[(3*vsize//4):])<threshold):\n        good_crop = False\n        bottom_edge = 3*vsize//4\n    else:\n        if len(np.where(vmeans[3*vsize//4:]<threshold)[0])>0:\n            bottom_edge = np.min(3*vsize//4+np.where(vmeans[3*vsize//4:]<threshold)[0])\n            if bottom_edge<(vsize-1):\n                bottom_edge += 1\n        else:\n            bottom_edge = (vsize-1)\n            \n    return rgb_pixels[top_edge:bottom_edge,left_edge:right_edge,:], good_crop\n\n\ndef extract_study_info(ct_study,ct_series,base_path,data_df=test_df):\n    \n    dcim_list = data_df.loc[ ((data_df['StudyInstanceUID']==ct_study) & (data_df['SeriesInstanceUID']==ct_series))]['SOPInstanceUID']\n    \n    ###########################################################################################\n    # Collect general metadata for the series\n    general_dcm = pydicom.dcmread(base_path/ct_study/ct_series/(dcim_list.iloc[0]+'.dcm'))\n    \n    # important info for picture interpretation\n    n_rows = general_dcm.Rows\n    n_cols = general_dcm.Columns\n    n_images = len(dcim_list)\n    \n    \n    # add metadata from image characteristics\n    mdt_dict = {kv:getattr(general_dcm,kv) for kv in series_metadata_fields}\n    mdt_dict['n_file_images']=n_images\n    \n    \n    \n    if isinstance(general_dcm.WindowWidth, pydicom.multival.MultiValue):\n        mdt_dict['WindowWidthDim'] = len(general_dcm.WindowWidth)\n        mdt_dict['WindowWidth'] = float(general_dcm.WindowWidth[0])\n        mdt_dict['WindowWidth2'] = float(general_dcm.WindowWidth[1])\n    else:\n        mdt_dict['WindowWidthDim'] = 1\n        mdt_dict['WindowWidth'] = float(general_dcm.WindowWidth)\n        mdt_dict['WindowWidth2'] = 0\n\n    if isinstance(general_dcm.WindowCenter, pydicom.multival.MultiValue):\n        mdt_dict['WindowCenterDim'] = len(general_dcm.WindowCenter)\n        mdt_dict['WindowCenter'] = float(general_dcm.WindowCenter[0])\n        mdt_dict['WindowCenter2'] = float(general_dcm.WindowCenter[1])\n    else:\n        mdt_dict['WindowCenterDim'] = 1\n        mdt_dict['WindowCenter'] = float(general_dcm.WindowCenter)\n        mdt_dict['WindowCenter2'] = 0\n\n\n    mdt_dict['PixelSpacing'] = float(general_dcm.PixelSpacing[0])\n    mdt_dict['PixelSpacing2'] = float(general_dcm.PixelSpacing[1])\n    \n    ###########################################################################################\n    # Collect image data for series\n    \n    # initialize variables\n    image_order_dict = {}\n    raw_pixel3d = np.zeros((n_images,n_rows,n_cols),dtype='int16')\n    lung3d = np.zeros((n_images,n_rows,n_cols),dtype='bool')\n    im_present_tf = np.zeros(n_images,dtype='bool')\n\n    # collect info for each image\n    for n,sop in enumerate(dcim_list):\n        sop_dcm = pydicom.dcmread(base_path/ct_study/ct_series/(sop+'.dcm'))\n        ordervalue = int(sop_dcm.get_item('InstanceNumber').value)-1\n        image_order_dict[sop] = ordervalue\n        # check and see if the instance number is greater than number of images\n        if ordervalue>=n_images:\n            \n            n_extra = int(sop_dcm.get_item('InstanceNumber').value) - n_images\n            raw_pixel3d = np.concatenate((raw_pixel3d,np.zeros((n_extra,n_rows,n_cols),dtype='int16')),axis=0)\n            lung3d = np.concatenate((lung3d,np.zeros((n_extra,n_rows,n_cols),dtype='bool')),axis=0)\n            im_present_tf = np.concatenate((im_present_tf,np.zeros(n_extra,dtype='bool')),axis=0)\n            n_images = ordervalue\n\n        raw_pixels = sop_dcm.pixel_array.astype('int') * int(sop_dcm.RescaleSlope) + int(sop_dcm.RescaleIntercept)\n        raw_pixel3d[ordervalue] = raw_pixels.astype('int16')\n        lung3d[ordervalue] = (signal.convolve2d((raw_pixels>-900)&(raw_pixels<-500),np.ones((6,6)),mode='same')>25)\n        im_present_tf[ordervalue] = True\n\n    # Crop unused space\n    raw_pixel3d, lung3d = variation_crop(raw_pixel3d,lung3d, threshold=100)\n    new_imagenum,new_rows,new_cols = raw_pixel3d.shape\n\n    # Lung Presence\n    central_lung = lung3d[:,int(lung3d.shape[1]*0.05):-1*int(lung3d.shape[1]*0.2),int(lung3d.shape[1]*0.05):-1*int(lung3d.shape[1]*0.03)]\n    lung_fraction = np.sum(np.sum(central_lung,axis=1),axis=1)/(central_lung.shape[1]*central_lung.shape[2])\n    \n    # additional metadata\n    mdt_dict['n_all_images']=lung3d.shape[0]\n    mdt_dict['n_crop_row']=lung3d.shape[1]\n    mdt_dict['n_crop_col']=lung3d.shape[2]\n    \n    return mdt_dict, dcim_list, image_order_dict, raw_pixel3d, im_present_tf, lung3d, lung_fraction \n\ndef lrc_lung_images(sop, artery_base, image_order_dict, im_present_tf, lung3d):\n    wt_vec = image_wt_vector(image_order_dict[sop],im_present_tf)\n    nonz = np.where(wt_vec)\n    wt_vec = wt_vec[nonz]\n    im = np.sum(wt_vec.reshape(len(wt_vec),1,1,1) * artery_base[nonz],axis=0)\n    im = im.astype('uint8')\n\n    cropped_sample,good_tf = lung_crop(im,lung3d[image_order_dict[sop]])\n    left_image = rounded_borders(cropped_sample[:,0:int(cropped_sample.shape[1]*0.45),:])\n    right_image = rounded_borders(cropped_sample[:,int(cropped_sample.shape[1]*0.55):,:])\n    center_image = rounded_borders(cropped_sample[:,int(cropped_sample.shape[1]*0.3):int(cropped_sample.shape[1]*0.75),:])\n    \n    return left_image, center_image, right_image\n    \n    \ndef lrc_heart_images(sop, heart_base, image_order_dict, im_present_tf):\n    wt_vec = image_wt_vector(image_order_dict[sop],im_present_tf)\n    nonz = np.where(wt_vec)\n    wt_vec = wt_vec[nonz]\n    im = np.sum(wt_vec.reshape(len(wt_vec),1,1) * heart_base[nonz],axis=0)\n    return im.astype('uint8')\n\n\ndef create_lung_images2(dcim_list,image_order_dict,raw_pixel3d,im_present_tf,\n                       ctbins,lung_fraction,lung3d,output_mode='list',filestem='tmp/',min_lung_fraction=0.03):\n        \n    all_images = []\n    lung_tf = np.zeros(len(dcim_list),dtype='bool')\n\n    # Artery Base Images\n    artery_base = artery_base_image(raw_pixel3d)\n    \n    for n,sop in enumerate(dcim_list):\n        if lung_fraction[image_order_dict[sop]]>min_lung_fraction:\n            lung_tf[n] = True\n            all_image = dask.delayed(lrc_lung_images)(sop, artery_base, image_order_dict, im_present_tf, lung3d)\n            all_images.append(all_image)\n\n    all_images = dask.compute(*all_images)\n    left_lung_images = [all_image[0] for all_image in all_images]\n    center_lung_images = [all_image[1] for all_image in all_images]\n    right_lung_images = [all_image[2] for all_image in all_images]\n            \n    return left_lung_images,center_lung_images,right_lung_images, lung_tf\n\ndef create_heart_images2(dcim_list,image_order_dict,raw_pixel3d,im_present_tf,ctbins,output_mode='list',filestem='/tmp'):\n\n    heart_images = []\n\n    new_imagenum,new_rows,new_cols = raw_pixel3d.shape\n    heart_base = np.digitize(raw_pixel3d,ctbins).astype('uint8')\n    half_im_size = int(0.5*min(0.8*new_cols,min(0.9*new_rows,0.5*(0.8*new_rows+0.6*new_cols))))\n    new_row_center = int(0.45*new_rows)\n    new_col_center = int(0.6*new_cols)\n    heart_base = heart_base[:,\n                            (new_row_center-half_im_size):(new_row_center+half_im_size),\n                            (new_col_center-half_im_size):(new_col_center+half_im_size)]\n    \n    for n,sop in enumerate(dcim_list):\n        heart_image = dask.delayed(lrc_heart_images)(sop, heart_base, image_order_dict, im_present_tf)\n        heart_images.append(heart_image)\n\n    heart_images = dask.compute(*heart_images)\n            \n    return heart_images\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def magic(ct_study,test_df):\n    ct_series_list = (test_df.loc[(test_df['StudyInstanceUID']==ct_study)]['SeriesInstanceUID']).unique()\n    ct_series = ct_series_list[0]\n\n\n\n    mdt_dict, dcim_list, image_order_dict, raw_pixel3d, im_present_tf, lung3d, lung_fraction  = extract_study_info(ct_study,\n                                                                                                                   ct_series,\n                                                                                                                   base_path=data_path/\"test\",\n                                                                                                                   data_df=test_df)\n\n    # takes 3.5 sec\n    L_images,C_images,R_images,TF_images = create_lung_images2(dcim_list,image_order_dict,raw_pixel3d,im_present_tf,\n                                        ctbins,lung_fraction,lung3d,output_mode='list',filestem='tmp/',\n                                        min_lung_fraction=MIN_LUNG_FRACTION)\n\n    heart_images = create_heart_images2(dcim_list,image_order_dict,raw_pixel3d,im_present_tf,ctbins,output_mode='list',filestem='/tmp')\n\n    all_image_features_df = pd.DataFrame(dcim_list.values,columns=['SOPInstanceUID'])\n    \n    all_image_features_df['SeriesInstanceUID'] = ct_series\n    all_image_features_df['StudyInstanceUID'] = ct_study\n\n    \n    dl = ventricle_model.dls.test_dl(heart_images, num_workers=8,bs=128)\n    with ventricle_model.no_bar():\n        preds,_ = ventricle_model.get_preds(dl=dl)\n    all_image_features_df['ventricle_pred'] = preds[:,1]\n\n    \n    all_image_features_df['lung_fraction'] = lung_fraction[im_present_tf]\n    all_image_features_df['image_order'] = list(image_order_dict.values())\n    all_image_features_df['top_vent_image'] = 0\n    all_image_features_df['top_lung_image'] = 0\n    \n    top_vent_index = np.argmax(preds[:,1].numpy())\n    \n    all_image_features_df.loc[top_vent_index,'top_vent_image'] = 1\n    all_image_features_df.loc[np.argmax(lung_fraction[im_present_tf]),'top_lung_image'] = 1\n    \n    \n    for lung_section in ['L','C','R']:\n        acute_pred = np.zeros(len(all_image_features_df)) + 0.01\n        chronic_pred = np.zeros(len(all_image_features_df)) + 0.01\n        negative_pred = np.zeros(len(all_image_features_df)) + 0.01\n        \n        if len(L_images)>0:\n            if lung_section == 'L':\n                dl = peclass_model.dls.test_dl(L_images, num_workers=8,bs=128)\n            elif lung_section == 'C':\n                dl = peclass_model.dls.test_dl(C_images, num_workers=8,bs=128)\n            elif lung_section == 'R':\n                dl = peclass_model.dls.test_dl(R_images, num_workers=8,bs=128)\n\n            with peclass_model.no_bar():\n                preds,_ = peclass_model.get_preds(dl=dl)\n            acute_pred[TF_images] = preds[:,0]\n            chronic_pred[TF_images] = preds[:,1]\n            negative_pred[TF_images] = preds[:,2]\n            \n        all_image_features_df['pe_acute_'+lung_section] = acute_pred\n        all_image_features_df['pe_chronic_'+lung_section] = chronic_pred\n        all_image_features_df['pe_neg_'+lung_section] = negative_pred\n\n    \n    contrast_pred = np.zeros(len(all_image_features_df)) + -1.0\n    motion_pred = np.zeros(len(all_image_features_df)) + -1.0\n    ventricle_pred = np.zeros(len(all_image_features_df)) + -1.0\n    \n    TF_heart = (all_image_features_df['ventricle_pred']>0.5)|(all_image_features_df['top_vent_image']==1)\n    if sum(TF_heart)>0:\n        select_heart_images = [np.repeat(heart_images[idx][:, :, np.newaxis], 3, axis=2) for idx in range(len(heart_images)) if TF_heart[idx]]\n\n        # contrast\n        dl = contrastclass_model.dls.test_dl(select_heart_images, num_workers=8,bs=128)\n        with contrastclass_model.no_bar():\n            preds,_ = contrastclass_model.get_preds(dl=dl)\n        contrast_pred[TF_heart] = preds[:,0]\n        \n\n        # motion\n        dl = motionclass_model.dls.test_dl(select_heart_images, num_workers=8,bs=128)\n        with motionclass_model.no_bar():\n            preds,_ = motionclass_model.get_preds(dl=dl)\n        motion_pred[TF_heart] = preds[:,0]\n        \n\n        # ventricle class\n        dl = ventclass_model.dls.test_dl(select_heart_images, num_workers=8,bs=128)\n        with ventclass_model.no_bar():\n            preds,_ = ventclass_model.get_preds(dl=dl)\n\n        ventricle_pred[TF_heart] = preds[:,1]\n        \n\n    all_image_features_df['contrast_pred'] = contrast_pred\n    all_image_features_df['motion_pred'] = motion_pred\n    all_image_features_df['rvlv_gte1_pred'] = ventricle_pred\n    \n    \n    ##############################################################################################################################\n    #######  Series features\n    series_features_df = pd.DataFrame(mdt_dict,index=[0])\n    series_features_df['body_crop_width'] = series_features_df['PixelSpacing2']*series_features_df['n_crop_col']\n    series_features_df['body_crop_height'] = series_features_df['PixelSpacing']*series_features_df['n_crop_row']\n    series_features_df['image_gaps'] = np.array(series_features_df['n_file_images']!=series_features_df['n_all_images'])*1\n\n\n    # Lung Image Order\n    if len(all_image_features_df.loc[(all_image_features_df['lung_fraction']>MIN_LUNG_FRACTION),'image_order'])>0:\n        min_lung_image = min( all_image_features_df.loc[(all_image_features_df['lung_fraction']>MIN_LUNG_FRACTION),'image_order'])\n        max_lung_image = max( all_image_features_df.loc[(all_image_features_df['lung_fraction']>MIN_LUNG_FRACTION),'image_order'])\n    else:\n        min_lung_image = min(all_image_features_df['image_order'])\n        max_lung_image = max(all_image_features_df['image_order'])\n\n    all_image_features_df['lung_order_pct'] = (all_image_features_df['image_order'].values - min_lung_image)/(max_lung_image - min_lung_image)\n\n    # Lung area relative to max lung area\n    lung_max = np.max(all_image_features_df['lung_fraction'])\n    all_image_features_df['lung_area_cm2'] = (all_image_features_df['lung_fraction']\n                                                             *(series_features_df['body_crop_width'][0]\n                                                               *series_features_df['body_crop_height'][0]))\n    all_image_features_df['lung_area_vs_max'] = all_image_features_df['lung_fraction']/lung_max\n\n    # Series Calculations\n    series_features_df['lung_fraction_max'] = lung_max\n    series_features_df['lung_area_max_cm2'] = 0.01*lung_max*series_features_df['body_crop_width'][0]*series_features_df['body_crop_height'][0]\n    series_features_df['lung_length_cm'] = 0.1*(max_lung_image - min_lung_image)*series_features_df['SliceThickness'][0]\n    series_features_df['lung_volume_cm3'] = 0.001*(np.sum(all_image_features_df['lung_fraction'])\n                                                        *series_features_df['body_crop_width'][0]\n                                                        *series_features_df['body_crop_height'][0]\n                                                        *series_features_df['SliceThickness'][0])\n\n    # Series aggregate statistics for classification models\n\n    #ventricle visibility\n    series_features_df['ventricle_visibility_max'] = np.max(all_image_features_df['ventricle_pred'])\n    series_features_df['ventricle_visibility_num_above75'] = np.sum(all_image_features_df['ventricle_pred']>0.75)\n\n\n    # pe classification\n    for pe_type in ['acute','chronic','neg']:\n        for pe_loc in ['L','C','R']:\n            colname = 'pe_'+pe_type+'_'+pe_loc\n            for stat_type in ['max','above20','above50','lungmean']:\n                if stat_type=='max':       \n                    stat_calc = np.max(all_image_features_df[colname])\n                elif stat_type=='above20':\n                    stat_calc = np.sum(all_image_features_df[colname]>0.2)/len(all_image_features_df)\n                elif stat_type=='above50':\n                    stat_calc = np.sum(all_image_features_df[colname]>0.5)/len(all_image_features_df)\n                elif stat_type=='lungmean':\n                    stat_calc = (np.sum(all_image_features_df[colname]*all_image_features_df['lung_fraction'])\n                             /np.sum(all_image_features_df['lung_fraction']) )\n\n                series_features_df[colname+'_'+stat_type] = stat_calc\n\n\n    high_vent_indices = set(all_image_features_df[all_image_features_df['ventricle_pred']>0.8].index)\n    high_vent_indices.add(top_vent_index)\n    high_vent_indices = list(high_vent_indices)\n\n    # ventricle right/left classification\n    series_features_df['rvlv_gte1_classification_top'] = all_image_features_df.loc[top_vent_index,'rvlv_gte1_pred']\n    series_features_df['rvlv_gte1_classification_min'] = np.min(all_image_features_df.loc[high_vent_indices,'rvlv_gte1_pred'])\n    series_features_df['rvlv_gte1_classification_max'] = np.max(all_image_features_df.loc[high_vent_indices,'rvlv_gte1_pred'])\n    series_features_df['rvlv_gte1_classification_mean'] = np.mean(all_image_features_df.loc[high_vent_indices,'rvlv_gte1_pred'])\n\n    # motion classification\n    series_features_df['motion_pred_top'] = all_image_features_df.loc[top_vent_index,'motion_pred']\n    series_features_df['motion_pred_min'] = np.min(all_image_features_df.loc[high_vent_indices,'motion_pred'])\n    series_features_df['motion_pred_max'] = np.max(all_image_features_df.loc[high_vent_indices,'motion_pred'])\n    series_features_df['motion_pred_mean'] = np.mean(all_image_features_df.loc[high_vent_indices,'motion_pred'])\n\n\n    # contrast classification\n    series_features_df['contrast_pred_top'] = all_image_features_df.loc[top_vent_index,'contrast_pred']\n    series_features_df['contrast_pred_min'] = np.min(all_image_features_df.loc[high_vent_indices,'contrast_pred'])\n    series_features_df['contrast_pred_max'] = np.max(all_image_features_df.loc[high_vent_indices,'contrast_pred'])\n    series_features_df['contrast_pred_mean'] = np.mean(all_image_features_df.loc[high_vent_indices,'contrast_pred'])\n\n    \n    all_image_features_df.loc[(all_image_features_df['lung_order_pct']<-1),'lung_order_pct']=-1\n    all_image_features_df.loc[(all_image_features_df['lung_order_pct']>2),'lung_order_pct']=2\n    all_image_features_df.loc[np.isnan(all_image_features_df['lung_order_pct']),'lung_order_pct']=-1\n    \n\n\n    pred_outcome = np.mean(np.array([gbm_outcome[fold].predict(series_features_df[feature_cols]) for fold in range(5)]),axis=0)\n    pred_vent =  np.mean(np.array([gbm_vent[fold].predict(series_features_df[feature_cols]) for fold in range(5)]),axis=0)\n    pred_left =  np.mean(np.array([gbm_left[fold].predict(series_features_df[feature_cols]) for fold in range(5)]),axis=0)\n    pred_center = np.mean(np.array([gbm_center[fold].predict(series_features_df[feature_cols]) for fold in range(5)]),axis=0)\n    pred_right =  np.mean(np.array([gbm_right[fold].predict(series_features_df[feature_cols]) for fold in range(5)]),axis=0)\n\n\n    for study_feature in feature_cols:\n        all_image_features_df[study_feature] = series_features_df[study_feature]\n\n    pe_present_on_image = np.mean(np.array([gbm_poi[fold].predict(all_image_features_df[direct_image_features+feature_cols]) for fold in range(5)]),axis=0)\n\n\n    negative_exam_for_pe = pred_outcome[0][1]\n    chronic_pe = pred_outcome[0][2]\n    acute_and_chronic_pe = pred_outcome[0][3]\n    indeterminate  = pred_outcome[0][4]\n\n    rv_lv_ratio_gte_1 = pred_vent[0][1]\n    rv_lv_ratio_lt_1 = pred_vent[0][2]\n\n    rightsided_pe = pred_right[0]\n    leftsided_pe = pred_left[0]\n    central_pe = pred_center[0]\n\n\n    if (negative_exam_for_pe+indeterminate)<0.5:  # should be positive\n        if max([rv_lv_ratio_gte_1,rv_lv_ratio_lt_1])<0.5:\n            xtra = 0.5001- max([rv_lv_ratio_gte_1,rv_lv_ratio_lt_1])\n            rv_lv_ratio_gte_1 += xtra\n            rv_lv_ratio_lt_1 += xtra\n        if max([rightsided_pe,leftsided_pe,central_pe])<0.5:\n            xtra = 0.5001- max([rightsided_pe,leftsided_pe,central_pe])\n            rightsided_pe += xtra\n            leftsided_pe += xtra\n            central_pe += xtra\n        if np.max(pe_present_on_image)<0.5:\n            xtra = 0.5001-np.max(pe_present_on_image)\n            pe_present_on_image += xtra\n    else:\n        if max([rv_lv_ratio_gte_1,rv_lv_ratio_lt_1])>0.5:\n            xtra = max([rv_lv_ratio_gte_1,rv_lv_ratio_lt_1])-0.5\n            rv_lv_ratio_gte_1 = max([0,rv_lv_ratio_gte_1-xtra])\n            rv_lv_ratio_lt_1 = max([0,rv_lv_ratio_lt_1-xtra])\n        if max([rightsided_pe,leftsided_pe,central_pe])>0.5:\n            xtra = max([rightsided_pe,leftsided_pe,central_pe]) - 0.5\n            rightsided_pe = max([0,rightsided_pe-xtra])\n            leftsided_pe = max([0,leftsided_pe-xtra])\n            central_pe =  max([0,central_pe-xtra])\n        if np.max(pe_present_on_image)>0.5:\n            pe_present_on_image[pe_present_on_image>0.5] = 0.49\n\n\n            \n    ids = [\n    ct_study+'_negative_exam_for_pe',\n    ct_study+'_rv_lv_ratio_gte_1',\n    ct_study+'_rv_lv_ratio_lt_1',\n    ct_study+'_leftsided_pe',\n    ct_study+'_chronic_pe',\n    ct_study+'_rightsided_pe',\n    ct_study+'_acute_and_chronic_pe',\n    ct_study+'_central_pe',\n    ct_study+'_indeterminate'] + list(dcim_list)\n    \n\n            \n    out_df = pd.DataFrame(ids,columns=['id'])\n    out_df['label'] = [negative_exam_for_pe,\n        rv_lv_ratio_gte_1,\n        rv_lv_ratio_lt_1,\n        leftsided_pe,\n        chronic_pe,\n        rightsided_pe,\n        acute_and_chronic_pe,\n        central_pe,\n        indeterminate]+list(pe_present_on_image)\n            \n    return out_df\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def sadness(ct_study,test_df):\n    ct_series_list = (test_df.loc[(test_df['StudyInstanceUID']==ct_study)]['SeriesInstanceUID']).unique()\n    ct_series = ct_series_list[0]\n    dcim_list = test_df.loc[ ((test_df['StudyInstanceUID']==ct_study) & (test_df['SeriesInstanceUID']==ct_series))]['SOPInstanceUID']\n    \n    ids = [\n    ct_study+'_negative_exam_for_pe',\n    ct_study+'_rv_lv_ratio_gte_1',\n    ct_study+'_rv_lv_ratio_lt_1',\n    ct_study+'_leftsided_pe',\n    ct_study+'_chronic_pe',\n    ct_study+'_rightsided_pe',\n    ct_study+'_acute_and_chronic_pe',\n    ct_study+'_central_pe',\n    ct_study+'_indeterminate'] + list(dcim_list)\n    \n    out_df = pd.DataFrame(ids,columns=['id'])\n    \n    out_df['label'] = [0.674681,\n                        0.129139,\n                        0.174612,\n                        0.212117,\n                        0.040115,\n                        0.257590,\n                        0.019920,\n                        0.055090,\n                        0.021569] + list(np.ones(len(dcim_list))*0.053915)\n    return out_df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ct_studies = test_df['StudyInstanceUID'].unique()\n\nout_df = magic(ct_studies[0],test_df)\n\nfor ct_study in tqdm(ct_studies[1:]):\n    \n    if (time.time() - start_time)<TIME_LIMIT:\n        tmp_df = magic(ct_study,test_df)\n    else:\n        tmp_df = sadness(ct_study,test_df)\n        \n    out_df = pd.concat([out_df,tmp_df])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"out_df.to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}