{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport matplotlib.pyplot as plt\nimport pydicom\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"cd ../input/second-annual-data-science-bowl/","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Dataset visualisation & understanding**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"validate=os.listdir(\"validate/validate\")\nvalidate=sorted(validate,key=lambda x:int(x))\nprint(validate)\nprint(\"\\n length :\",len(validate))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"validate=os.listdir(\"test/test\")\nvalidate=sorted(validate,key=lambda x:int(x))\nprint(validate)\nprint(\"\\n length :\",len(validate))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"validate=os.listdir(\"train/train\")\nvalidate=sorted(validate,key=lambda x:int(x))\nprint(validate)\nprint(\"\\n length :\",len(validate))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n# Data to plot\nlabels = 'Train', 'test', 'validate'\nsizes = [500, 440, 200]\ncolors = ['yellowgreen', 'lightcoral', 'lightskyblue']\nexplode = (0, 0, 0)  # explode 1st slice\n\n# Plot\nplt.pie(sizes, explode=explode, labels=labels, colors=colors,\nautopct='%1.1f%%', shadow=True, startangle=140)\n\nplt.axis('equal')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndf=pd.read_csv(\"train.csv\")\nprint(df)\n ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"path=\"train/train/\"\npatients=os.listdir(path)\nnum_files=[]\nsax_files=[]\nfor patient in patients:\n    count=0\n    patient_files=os.listdir(path+patient+\"/study\")\n    num_files.append(len(patient_files))\n    for file in patient_files:\n        if file[:3]==\"sax\":\n            count+=1\n    sax_files.append(count)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(20,9))\nax = fig.add_axes([0,0,1,1])\nax.bar(patients[:50],num_files[:50])\nax.set_title('number of files per patient')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_sax=[]\nfor patient in patients:\n    files=os.listdir(path+patient+\"/study\")\n    sax_files=[ss for ss in files if ss[0]==\"s\"]\n    all_sax=all_sax+sax_files\nsax=dict()\nfor s in all_sax:\n    if s not in sax:\n        sax[s]=all_sax.count(s)\n\n        \nall_sax2=list(sax.keys())\nall_sax2=sorted(all_sax2,key=lambda x:int(x[4:]))\nprint(all_sax2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pos=[int(s[4:]) for s in all_sax2]\nval=sax.values()\nfig = plt.figure()\nax = fig.add_axes([0,0,1,1])\nax.bar(pos,val)\nax.set_title('number of slices by position')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"files=os.listdir(path+patient+\"/study\")\nslices_paths=os.listdir(path+patient+\"/study/\"+files[0])\na_slice=pydicom.dcmread(path+patient+\"/study/\"+files[0]+\"/\"+slices_paths[0])\nplt.imshow(a_slice.pixel_array,cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig=plt.figure(figsize=(20, 20))\ncolumns=6\nrows=5\nfor i in range(1,columns*rows+1):\n    fig.add_subplot(rows,columns,i)\n    slices_paths=sorted(slices_paths,key=lambda x: int(x[-8:-4]))\n    a_slice=pydicom.dcmread(path+patient+\"/study/\"+files[0]+\"/\"+slices_paths[i-1])\n    print(path+patient+\"/study/\"+files[0]+\"/\"+slices_paths[i-1])\n    plt.imshow(a_slice.pixel_array,cmap=plt.cm.bone)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **ROI extraction**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import cv2\nimport pydicom\nimport matplotlib.pyplot as plt\n\n\n\npath=\"train/train\"\npatients_list=os.listdir(path) #self explanatory\n\n\ndef patient_folders(patient):\n    \"\"\"returns all the short axis slices folders, sorted by position from sax_MIN to sax_MAX\"\"\"\n    files = os.listdir(path+\"/\"+patient+\"/study\")\n    files = [f for f in files if f[0]==\"s\"]                           #remove 2 and 4 chamber views\n    files = sorted(files, key = lambda x: int(x.split(\"_\")[1]))       #sort the folders on spacial position\n    return files\n\n\n\n\ndef normalize_image(dcm_slice):\n    \"\"\"returns a normalized image where pixels are between 0-255 and each pixel is one mm (original size of the images)\"\"\"\n    \n    image = dcm_slice.pixel_array\n    \n    scale = dcm_slice.PixelSpacing\n    \n    new_size = (int(dcm_slice.Rows*scale[0]), int(dcm_slice.Columns*scale[1]))\n    \n    #normalizing the image\n    \n    img_2d = image.astype(float)\n    img_2d_scaled = (np.maximum(img_2d,0) / img_2d.max()) * 255.0\n    image = np.uint8(img_2d_scaled)\n    \n    \n    #resizing the image to real dimensions (mm)\n    \n    image = cv2.resize(image, dsize=new_size, interpolation=cv2.INTER_CUBIC)\n    \n    return image\n\n\n\ndef create_stack(patient,folder):\n    \"\"\"creates a 3D matrix of 30 slices from a single folder, that will be used to define the ROI\"\"\"\n    \n    stack=[]\n    if folder in patient_folders(patient):\n        \n        slices_path = \"train/train/\"+patient+\"/study/\"+folder\n        slices_names = os.listdir(slices_path)\n        \n        slices_names=sorted(slices_names,key=lambda x: int(x[-8:-4]))\n        for s in slices_names:\n            \n            dcm_slice = pydicom.dcmread(slices_path+\"/\"+s)\n            image = normalize_image(dcm_slice)\n            stack.append(image)\n        \n        image_stack = np.dstack(stack)\n    \n        return image_stack\n        \n    else:\n        print(\"error, folder and patient don't match, consider calling print(patient_folders('{}'))\".format(patient))\n        \n        \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def crop_ROI(patient, file, offset=35):\n    \"\"\"crops the region of interest in our images returning a smaller image focusing on the heart using standard deviation\n    to find the region of the heart exploiting the movement of muscles\"\"\"\n    \n    #creat a stack of 30 images of the same section in diffrent frames\n\n    patient_stack = create_stack(patient,file)             \n    \n    #calculates the standard deviation of the stack, returns an image containing only the moving pixels during the 30 frames\n                                          \n    std_image = np.std(patient_stack, axis=2) \n    \n    #normalizing the image again to have a range of 0-255\n    \n    img_2d         = std_image.astype(float)\n    img_2d_scaled  = (np.maximum(img_2d,0) / img_2d.max()) * 255.0\n    movement_image = np.uint8(img_2d_scaled)\n    \n    #applying gaussian filter for noise reduction\n    \n    movement_image = cv2.blur(movement_image,(6,6))\n    \n    #applying Canny edge detection to find the edges of the heart \n    \n    edge = cv2.Canny(movement_image,50,100)\n\n    #blurring again\n    \n    edge = cv2.GaussianBlur(edge,(3,3),cv2.BORDER_DEFAULT)\n\n    #using Hough transform to find an approximation circular patterns (left ventricule) position and radius\n    \n    circles = cv2.HoughCircles(edge,cv2.HOUGH_GRADIENT,1,100,param1=100,param2=34,minRadius=0,maxRadius=0)\n    circles = np.uint16(np.around(circles))\n    \n    #unpacking the coordinates of the circle\n    \n    x,y,r = circles[0][0][0], circles[0][0][1], circles[0][0][2]\n\n    #defining the ROI box\n    \n    x_1 = x-r-offset\n    y_1 = y-r-offset\n    \n    x_2 = x+r+offset\n    y_2 = y+r+offset\n    \n    return (x_1, x_2, y_1, y_2)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_cropped(patient, file):\n    \"\"\" displays the 30 slices cropped ROI \"\"\"\n    \n    slices_path = \"train/train/\"+patient+\"/study/\"\n\n    #verifies if the file is ammong the patient's files\n    \n    if file not in os.listdir(slices_path):\n        print(\"Error: {} file not found in this directory : {}\".format(file,slices_path))\n        return 1\n    \n    #this bloc is executed if the file is found\n    \n    #finds the file names\n    \n    slices_path = \"train/train/\"+patient+\"/study/\"+file\n    slices_names = os.listdir(slices_path)\n    \n    #sort the slices through time\n    \n    slices_names=sorted(slices_names,key=lambda x: int(x[-8:-4]))\n    \n    #creates a list of relative paths to the slices\n    \n    slices = [slices_path+\"/\"+s for s in slices_names]\n    \n    #check if we want to view the cropped version\n    \n    x1, x2, y1, y2 = crop_ROI(patient, file)\n    \n    \n    \n    for s in slices:\n        \n        my_slice = pydicom.dcmread(s)\n        \n        #normalize the image , pixels range 0-255\n        image = normalize_image(my_slice)\n        \n        #create the cropped image\n        cropped = image[y1:y2, x1:x2]\n        \n        #create the plot\n        \n        \n        \n        plt.imshow(image,cmap=\"bone\")\n        plt.title(\"original\")\n        plt.show()\n        \n        plt.imshow(cropped,cmap=\"bone\")\n        plt.title(\"cropped image\")\n        \n        plt.show()\n        \n        \n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patient_folders(\"85\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display_cropped(\"85\",\"sax_10\")","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}