{"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":"<center><h1>Learning and Applying the Fourier transform to the BraTS21' dataset</h1></center>\n\n## Source of Learning:\n * [Prof. Sreeni Video](https://www.youtube.com/watch?v=9mLeVn8xzMw&ab_channel=Apeer_micro)\n * [Sreeni github](https://github.com/bnsreenu/python_for_image_processing_APEER/blob/master/tutorial41_image_filters_using_fourier_transform_DFT.py)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"### This notebook acts as three things:\n* a follow along as to how to apply a fourier transform to an image(video in sources)\n* proof of concept that the principle works just fine on an mri image\n* base for moving to next notebook which will aim to glean more value from the data using the transform more explicitly in our context","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport os\nimport glob\nimport re\n\nimport cv2\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.826633Z","iopub.execute_input":"2021-09-10T11:46:05.827511Z","iopub.status.idle":"2021-09-10T11:46:05.832172Z","shell.execute_reply.started":"2021-09-10T11:46:05.827477Z","shell.execute_reply":"2021-09-10T11:46:05.831502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"hidden cell for initial attempts at more graceful sorting(need to come back to)","metadata":{}},{"cell_type":"code","source":"def tryint(s):\n    try:\n        return int(s)\n    except ValueError:\n        return s\n    \ndef alphanum_key(s):\n    \"\"\" Turn a string into a list of string and number chunks.\n        \"z23a\" -> [\"z\", 23, \"a\"]\n    \"\"\"\n    return [ tryint(c) for c in re.split('([0-9]+)', s) ]\n\ndef sort_nicely(l):\n    \"\"\" Sort the given list in the way that humans expect.\n    \"\"\"\n    l.sort(key=alphanum_key)\n#source: https://nedbatchelder.com/blog/200712/human_sorting.html","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-10T11:46:05.836241Z","iopub.execute_input":"2021-09-10T11:46:05.836536Z","iopub.status.idle":"2021-09-10T11:46:05.846916Z","shell.execute_reply.started":"2021-09-10T11:46:05.836498Z","shell.execute_reply":"2021-09-10T11:46:05.846231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# flair_zero = sorted(glob_values, key= (lambda x:  alphanum_key(glob_values[x])))","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.867347Z","iopub.execute_input":"2021-09-10T11:46:05.867799Z","iopub.status.idle":"2021-09-10T11:46:05.870927Z","shell.execute_reply.started":"2021-09-10T11:46:05.867761Z","shell.execute_reply":"2021-09-10T11:46:05.870372Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"glob_values = glob.glob(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/test/00001/FLAIR/*\")","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.848287Z","iopub.execute_input":"2021-09-10T11:46:05.848688Z","iopub.status.idle":"2021-09-10T11:46:05.865921Z","shell.execute_reply.started":"2021-09-10T11:46:05.848643Z","shell.execute_reply":"2021-09-10T11:46:05.865030Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"glob_values[0]","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.872529Z","iopub.execute_input":"2021-09-10T11:46:05.872997Z","iopub.status.idle":"2021-09-10T11:46:05.884746Z","shell.execute_reply.started":"2021-09-10T11:46:05.872957Z","shell.execute_reply":"2021-09-10T11:46:05.884145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## load image, resize, crop","metadata":{}},{"cell_type":"code","source":"dicom = pydicom.read_file(glob_values[0])\ndata = dicom.pixel_array\ndata = cv2.resize(data, (256,256,)) #should edit this to be correct dimensions for dft shift","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.897951Z","iopub.execute_input":"2021-09-10T11:46:05.898737Z","iopub.status.idle":"2021-09-10T11:46:05.907559Z","shell.execute_reply.started":"2021-09-10T11:46:05.898697Z","shell.execute_reply":"2021-09-10T11:46:05.906553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(data)","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:05.926310Z","iopub.execute_input":"2021-09-10T11:46:05.927027Z","iopub.status.idle":"2021-09-10T11:46:06.102066Z","shell.execute_reply.started":"2021-09-10T11:46:05.926992Z","shell.execute_reply":"2021-09-10T11:46:06.101475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y = np.expand_dims(data, axis=-1)\nplt.imshow(y)\ny.shape","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.103315Z","iopub.execute_input":"2021-09-10T11:46:06.103702Z","iopub.status.idle":"2021-09-10T11:46:06.329656Z","shell.execute_reply.started":"2021-09-10T11:46:06.103660Z","shell.execute_reply":"2021-09-10T11:46:06.328761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from PIL import Image\nim = Image.fromarray(data)\nim.save(\"image.png\")\nimg = cv2.imread('image.png',0)","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.331000Z","iopub.execute_input":"2021-09-10T11:46:06.331305Z","iopub.status.idle":"2021-09-10T11:46:06.343622Z","shell.execute_reply.started":"2021-09-10T11:46:06.331266Z","shell.execute_reply":"2021-09-10T11:46:06.342719Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"crop_img = img[40:160, 60:200]\nplt.imshow(crop_img)","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.346294Z","iopub.execute_input":"2021-09-10T11:46:06.346882Z","iopub.status.idle":"2021-09-10T11:46:06.585496Z","shell.execute_reply.started":"2021-09-10T11:46:06.346829Z","shell.execute_reply":"2021-09-10T11:46:06.584571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"hidden cell below for testing the array values with different loading method","metadata":{}},{"cell_type":"code","source":"# im = Image.fromarray(crop_img)\n# im.save(\"image.png\")\n# img = cv2.imread('image.png',0)\n# dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)\n# dft","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-10T11:46:06.587037Z","iopub.execute_input":"2021-09-10T11:46:06.587458Z","iopub.status.idle":"2021-09-10T11:46:06.592117Z","shell.execute_reply.started":"2021-09-10T11:46:06.587398Z","shell.execute_reply":"2021-09-10T11:46:06.591126Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## for filtering out low frequency","metadata":{}},{"cell_type":"code","source":"dft = cv2.dft(np.float32(crop_img), flags=cv2.DFT_COMPLEX_OUTPUT)\ndft.shape #comparing sizes","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-10T11:46:06.593499Z","iopub.execute_input":"2021-09-10T11:46:06.593731Z","iopub.status.idle":"2021-09-10T11:46:06.607199Z","shell.execute_reply.started":"2021-09-10T11:46:06.593706Z","shell.execute_reply":"2021-09-10T11:46:06.606341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dft_shift = np.fft.fftshift(dft)\nmagnitude_spectrum = 20 * np.log(cv2.magnitude(dft_shift[:, :, 0], dft_shift[:, :, 1]))\ndft_shift.shape #comparing sizes","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-10T11:46:06.608757Z","iopub.execute_input":"2021-09-10T11:46:06.609003Z","iopub.status.idle":"2021-09-10T11:46:06.619980Z","shell.execute_reply.started":"2021-09-10T11:46:06.608978Z","shell.execute_reply":"2021-09-10T11:46:06.619066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### resize image to work with script","metadata":{}},{"cell_type":"code","source":"img = cv2.resize(data, (140,120,))","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.621371Z","iopub.execute_input":"2021-09-10T11:46:06.621612Z","iopub.status.idle":"2021-09-10T11:46:06.633747Z","shell.execute_reply.started":"2021-09-10T11:46:06.621586Z","shell.execute_reply":"2021-09-10T11:46:06.632781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## FFT Processing Step (needs more customization for quality usage; acts as a proof of concept for now)","metadata":{}},{"cell_type":"code","source":"rows, cols = img.shape\ncrow, ccol = int(rows / 2), int(cols / 2)\n\nmask = np.ones((rows, cols, 2), np.uint8)\n# r = 80\n# center = [crow, ccol]\n# x, y = np.ogrid[:rows, :cols]\n# mask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r*r\n# mask[mask_area] = 0\n\n\n# Circular LPF mask, center circle is 1, remaining all zeros\n# Only allows low frequency components - smooth regions\n#Can smooth out noise but blurs edges.\n#\n\"\"\"\nrows, cols = img.shape\ncrow, ccol = int(rows / 2), int(cols / 2)\nmask = np.zeros((rows, cols, 2), np.uint8)\nr = 100\ncenter = [crow, ccol]\nx, y = np.ogrid[:rows, :cols]\nmask_area = (x - center[0]) ** 2 + (y - center[1]) ** 2 <= r*r\nmask[mask_area] = 1\n# Band Pass Filter - Concentric circle mask, only the points living in concentric circle are ones\nrows, cols = img.shape\ncrow, ccol = int(rows / 2), int(cols / 2)\nmask = np.zeros((rows, cols, 2), np.uint8)\nr_out = 80\nr_in = 10\ncenter = [crow, ccol]\nx, y = np.ogrid[:rows, :cols]\nmask_area = np.logical_and(((x - center[0]) ** 2 + (y - center[1]) ** 2 >= r_in ** 2),\n                           ((x - center[0]) ** 2 + (y - center[1]) ** 2 <= r_out ** 2))\nmask[mask_area] = 1\n\"\"\"\n\n\n# apply mask and inverse DFT: Multiply fourier transformed image (values)\n#with the mask values. \nfshift = dft_shift * mask #dft_shift\n\n#Get the magnitude spectrum (only for plotting purposes)\nfshift_mask_mag = 20 * np.log(cv2.magnitude(fshift[:, :, 0], fshift[:, :, 1]))\n\n#Inverse shift to shift origin back to top left.\nf_ishift = np.fft.ifftshift(fshift)\n\n#Inverse DFT to convert back to image domain from the frequency domain. \n#Will be complex numbers\nimg_back = cv2.idft(f_ishift)\n\n#Magnitude spectrum of the image domain\nimg_back = cv2.magnitude(img_back[:, :, 0], img_back[:, :, 1])","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.634938Z","iopub.execute_input":"2021-09-10T11:46:06.635615Z","iopub.status.idle":"2021-09-10T11:46:06.647498Z","shell.execute_reply.started":"2021-09-10T11:46:06.635580Z","shell.execute_reply":"2021-09-10T11:46:06.646464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## plotting step","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(12, 12))\nax1 = fig.add_subplot(2,2,1)\nax1.imshow(img, cmap='gray')\nax1.title.set_text('Input Image')\nax2 = fig.add_subplot(2,2,2)\nax2.imshow(magnitude_spectrum, cmap='gray')\nax2.title.set_text('FFT of image')\nax3 = fig.add_subplot(2,2,3)\nax3.imshow(fshift_mask_mag, cmap='gray')\nax3.title.set_text('FFT + Mask')\nax4 = fig.add_subplot(2,2,4)\nax4.imshow(img_back, cmap='gray')\nax4.title.set_text('After inverse FFT')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-10T11:46:06.649224Z","iopub.execute_input":"2021-09-10T11:46:06.649476Z","iopub.status.idle":"2021-09-10T11:46:07.329485Z","shell.execute_reply.started":"2021-09-10T11:46:06.649449Z","shell.execute_reply":"2021-09-10T11:46:07.328594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}