{"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":"code","source":"import gc\nimport os\nimport cv2\nimport zipfile\nimport rasterio\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport tifffile\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nfrom rasterio.windows import Window\nfrom torch.utils.data import Dataset\n\n!pip install staintools\n!pip install spams\n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:06:15.738297Z","iopub.execute_input":"2022-09-15T04:06:15.739031Z","iopub.status.idle":"2022-09-15T04:08:40.960603Z","shell.execute_reply.started":"2022-09-15T04:06:15.738914Z","shell.execute_reply":"2022-09-15T04:08:40.959052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imaging_measurements = {\n    'hpa': {\n        'pixel_size': {\n            'kidney': 0.4,\n            'prostate': 0.4,\n            'largeintestine': 0.4,\n            'spleen': 0.4,\n            'lung': 0.4\n        },\n        'tissue_thickness': {\n            'kidney': 4,\n            'prostate': 4,\n            'largeintestine': 4,\n            'spleen': 4,\n            'lung': 4\n        }\n    },\n    'hubmap': {\n        'pixel_size': {\n            'kidney': 0.5,\n            'prostate': 6.263,\n            'largeintestine': 0.229,\n            'spleen': 0.4945,\n            'lung': 0.7562\n        },\n        'tissue_thickness': {\n        'kidney': 10,\n            'prostate': 5,\n            'largeintestine': 8,\n            'spleen': 4,\n            'lung': 5\n        }\n    }\n}","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:08:40.963821Z","iopub.execute_input":"2022-09-15T04:08:40.965142Z","iopub.status.idle":"2022-09-15T04:08:40.975287Z","shell.execute_reply.started":"2022-09-15T04:08:40.965075Z","shell.execute_reply":"2022-09-15T04:08:40.974013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import staintools\n\ntarget = cv2.imread(\"/kaggle/input/hubmap-organ-segmentation/test_images/10078.tiff\")\ntarget = staintools.LuminosityStandardizer.standardize(target)\nnormalizer = staintools.StainNormalizer(method='vahadane')\nnormalizer.fit(target)\n\n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:08:40.976681Z","iopub.execute_input":"2022-09-15T04:08:40.977037Z","iopub.status.idle":"2022-09-15T04:08:52.821357Z","shell.execute_reply.started":"2022-09-15T04:08:40.976998Z","shell.execute_reply":"2022-09-15T04:08:52.820326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"OUT_TRAIN = [\"./hubmap_images\",\"./fda_images\"]\nos.makedirs(OUT_TRAIN[0], exist_ok=False)\nos.makedirs(OUT_TRAIN[1], exist_ok=False)\nMASKS = '../input/hubmap-organ-segmentation/train.csv'\nDATA = '../input/hubmap-organ-segmentation/train_images'","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:08:57.462598Z","iopub.execute_input":"2022-09-15T04:08:57.463049Z","iopub.status.idle":"2022-09-15T04:08:57.469755Z","shell.execute_reply.started":"2022-09-15T04:08:57.463009Z","shell.execute_reply":"2022-09-15T04:08:57.468459Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def augment_image(image, domain_pixel_size, target_pixel_size, domain_tissue_thickness, target_tissue_thickness, alpha=0.15):\n    \n    \"\"\"\n    Visualize raw and augmented images \n    \n    Parameters\n    ----------\n    image (numpy.ndarray of shape (height, width, 3)): Image array\n    domain_pixel_size (float): Pixel size of the domain images in micrometers\n    target_pixel_size (float): Pixel size of the target images in micrometers\n    domain_tissue_thickness (float): Tissue thickness of the domain images in micrometers\n    target_tissue_thickness (float): Tissue thickness of the target images in micrometers\n    alpha (float): Multiplier to control saturation and value scale\n    \"\"\"\n    \n    # Augment tissue thickness\n    tissue_thickness_scale_factor = target_tissue_thickness - domain_tissue_thickness\n    image_hsv = cv2.cvtColor(image, cv2.COLOR_RGB2HSV).astype(np.float32)\n    image_hsv[:, :, 1] *= (1 + (alpha * tissue_thickness_scale_factor))\n    image_hsv[:, :, 2] *= (1 - (alpha * tissue_thickness_scale_factor))\n    image_hsv = image_hsv.astype(np.uint8)\n    image_scaled = cv2.cvtColor(image_hsv, cv2.COLOR_HSV2RGB)\n    \n    # Standardize luminosity\n    image_scaled = staintools.LuminosityStandardizer.standardize(image_scaled)\n\n    # Augment pixel size\n    pixel_size_scale_factor = domain_pixel_size / target_pixel_size\n    image_resized = cv2.resize(\n        image_scaled,\n        dsize=None,\n        fx=pixel_size_scale_factor,\n        fy=pixel_size_scale_factor,\n        interpolation=cv2.INTER_CUBIC\n    )\n    image_resized = cv2.resize(\n        image_resized,\n        dsize=(\n            image.shape[1],\n            image.shape[0]\n        ),\n        interpolation=cv2.INTER_CUBIC\n    )\n    \n    # Standardize luminosity\n    image = staintools.LuminosityStandardizer.standardize(image)\n    image_augmented = staintools.LuminosityStandardizer.standardize(image_resized)\n    \n    \n    \n    \n    return image, image_augmented\n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:09:00.885282Z","iopub.execute_input":"2022-09-15T04:09:00.885896Z","iopub.status.idle":"2022-09-15T04:09:00.895655Z","shell.execute_reply.started":"2022-09-15T04:09:00.885863Z","shell.execute_reply":"2022-09-15T04:09:00.894460Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nfrom PIL import Image\nimport torch\nimport cv2\nimport torch\nimport numpy as np\n\ndef extract_ampl_phase(fft_im):\n    # fft_im: size should be bx3xhxwx2\n    fft_amp = fft_im[:,:,:,:,0]**2 + fft_im[:,:,:,:,1]**2\n    fft_amp = torch.sqrt(fft_amp)\n    fft_pha = torch.atan2( fft_im[:,:,:,:,1], fft_im[:,:,:,:,0] )\n    return fft_amp, fft_pha\n\ndef low_freq_mutate( amp_src, amp_trg, L=0.1 ):\n    _, _, h, w = amp_src.size()\n    b = (  np.floor(np.amin((h,w))*L)  ).astype(int)     # get b\n    amp_src[:,:,0:b,0:b]     = amp_trg[:,:,0:b,0:b]      # top left\n    amp_src[:,:,0:b,w-b:w]   = amp_trg[:,:,0:b,w-b:w]    # top right\n    amp_src[:,:,h-b:h,0:b]   = amp_trg[:,:,h-b:h,0:b]    # bottom left\n    amp_src[:,:,h-b:h,w-b:w] = amp_trg[:,:,h-b:h,w-b:w]  # bottom right\n    return amp_src\n\ndef low_freq_mutate_np( amp_src, amp_trg, L=0.1 ):\n    a_src = np.fft.fftshift( amp_src, axes=(-2, -1) )\n    a_trg = np.fft.fftshift( amp_trg, axes=(-2, -1) )\n\n    _, h, w = a_src.shape\n    b = (  np.floor(np.amin((h,w))*L)  ).astype(int)\n    c_h = np.floor(h/2.0).astype(int)\n    c_w = np.floor(w/2.0).astype(int)\n\n    h1 = c_h-b\n    h2 = c_h+b+1\n    w1 = c_w-b\n    w2 = c_w+b+1\n\n    a_src[:,h1:h2,w1:w2] = a_trg[:,h1:h2,w1:w2]\n    a_src = np.fft.ifftshift( a_src, axes=(-2, -1) )\n    return a_src\n\ndef FDA_source_to_target(src_img, trg_img, L=0.1):\n    # exchange magnitude\n    # input: src_img, trg_img\n\n    # get fft of both source and target\n    fft_src = torch.rfft( src_img.clone(), signal_ndim=2, onesided=False )\n    fft_trg = torch.rfft( trg_img.clone(), signal_ndim=2, onesided=False )\n\n    # extract amplitude and phase of both ffts\n    amp_src, pha_src = extract_ampl_phase( fft_src.clone())\n    amp_trg, pha_trg = extract_ampl_phase( fft_trg.clone())\n\n    # replace the low frequency amplitude part of source with that from target\n    amp_src_ = low_freq_mutate( amp_src.clone(), amp_trg.clone(), L=L )\n\n    # recompose fft of source\n    fft_src_ = torch.zeros( fft_src.size(), dtype=torch.float )\n    fft_src_[:,:,:,:,0] = torch.cos(pha_src.clone()) * amp_src_.clone()\n    fft_src_[:,:,:,:,1] = torch.sin(pha_src.clone()) * amp_src_.clone()\n\n    # get the recomposed image: source content, target style\n    _, _, imgH, imgW = src_img.size()\n    src_in_trg = torch.irfft( fft_src_, signal_ndim=2, onesided=False, signal_sizes=[imgH,imgW] )\n\n    return src_in_trg\n\ndef FDA_source_to_target_np( src_img, trg_img, L=0.1 ):\n    # exchange magnitude\n    # input: src_img, trg_img\n\n    src_img_np = src_img #.cpu().numpy()\n    trg_img_np = trg_img #.cpu().numpy()\n\n    # get fft of both source and target\n    fft_src_np = np.fft.fft2( src_img_np, axes=(-2, -1) )\n    fft_trg_np = np.fft.fft2( trg_img_np, axes=(-2, -1) )\n\n    # extract amplitude and phase of both ffts\n    amp_src, pha_src = np.abs(fft_src_np), np.angle(fft_src_np)\n    amp_trg, pha_trg = np.abs(fft_trg_np), np.angle(fft_trg_np)\n\n    # mutate the amplitude part of source with target\n    amp_src_ = low_freq_mutate_np( amp_src, amp_trg, L=L )\n\n    # mutated fft of source\n    fft_src_ = amp_src_ * np.exp( 1j * pha_src )\n\n    # get the mutated image\n    src_in_trg = np.fft.ifft2( fft_src_, axes=(-2, -1) )\n    src_in_trg = np.real(src_in_trg)\n\n    return src_in_trg\n\ndef transfer(im_src,target):\n    im_trg = target\n\n    im_src = np.asarray(im_src, np.float32)\n    im_trg = np.asarray(im_trg, np.float32)\n\n    im_src = im_src.transpose((2, 0, 1))\n    im_trg = im_trg.transpose((2, 0, 1))\n\n    src_in_trg = FDA_source_to_target_np( im_src, im_trg, L=0.001 )\n\n    src_in_trg = src_in_trg.transpose((1,2,0))\n\n    # scipy.misc.toimage(src_in_trg, cmin=0.0, cmax=255.0) # .save('demo_images/src_in_tar.png')\n\n    return src_in_trg # scipy.misc.toimage(src_in_trg, cmin=0.0, cmax=255.0)\n\n# if __name__ == \"__main__\":\n#     source=Image.open(\"demo_images/source.png\")\n#     img = transfer(source)\n#     img.save('demo_images/src_in_tar.png')\n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:09:07.037167Z","iopub.execute_input":"2022-09-15T04:09:07.037591Z","iopub.status.idle":"2022-09-15T04:09:07.066728Z","shell.execute_reply.started":"2022-09-15T04:09:07.037546Z","shell.execute_reply":"2022-09-15T04:09:07.065210Z"},"jupyter":{"source_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train = pd.read_csv(MASKS)\ndf_masks = df_train.reset_index(drop=True)\nimage_size = 1024\nfda_target = cv2.resize(target,(image_size,image_size),\n                            interpolation = cv2.INTER_LINEAR)\nfor idx, row in df_masks.iterrows():\n            #image+mask dataset\n            image = cv2.imread(os.path.join(DATA,str(row[\"id\"])+'.tiff'))\n            \n            image, image_augmented = augment_image(\n            image=image,\n            domain_pixel_size=imaging_measurements['hpa']['pixel_size'][row['organ']],\n            target_pixel_size=imaging_measurements['hubmap']['pixel_size'][row['organ']],\n            domain_tissue_thickness=imaging_measurements['hpa']['tissue_thickness'][row['organ']],\n            target_tissue_thickness=imaging_measurements['hubmap']['tissue_thickness'][row['organ']],\n            alpha=0.15\n        )\n            hub_image = normalizer.transform(image_augmented)\n            hub_image = cv2.resize(hub_image,(image_size,image_size),\n                            interpolation = cv2.INTER_LINEAR)\n            \n            fda_image = transfer(cv2.resize(image_augmented,(image_size,image_size),\n                            interpolation = cv2.INTER_LINEAR),fda_target)\n            \n            cv2.imwrite(os.path.join (OUT_TRAIN[0],str(row[\"id\"])+\".png\"),hub_image)\n            cv2.imwrite(os.path.join (OUT_TRAIN[1],str(row[\"id\"])+\".png\"),fda_image)\n\n           \n","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:09:07.541043Z","iopub.execute_input":"2022-09-15T04:09:07.541870Z","iopub.status.idle":"2022-09-15T04:11:07.123652Z","shell.execute_reply.started":"2022-09-15T04:09:07.541827Z","shell.execute_reply":"2022-09-15T04:11:07.121290Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_augmentation(image, image_augmented, metadata, path=None):\n\n    \"\"\"\n    Visualize raw and augmented images \n    \n    Parameters\n    ----------\n    image (numpy.ndarray of shape (height, width, 3)): Raw image array\n    image_augmented (numpy.ndarray of shape (height, width, 3)): Augmented image array\n    metadata (dict): Dictionary of image metadata\n    path (path-like str or None): Path of the output file or None (if path is None, plot is displayed with selected backend)\n    \"\"\"\n\n    fig, axes = plt.subplots(figsize=(36, 20), ncols=2)\n    \n    axes[0].imshow(image)\n    axes[1].imshow(image_augmented)\n\n    for i in range(2):\n        axes[i].set_xlabel('')\n        axes[i].set_ylabel('')\n        axes[i].tick_params(axis='x', labelsize=15, pad=10)\n        axes[i].tick_params(axis='y', labelsize=15, pad=10)\n\n    axes[0].set_title(f'Raw Image\\nMean: {np.mean(image):.4f} - Std: {np.std(image):.4f} - Min: {np.min(image)} - Max: {np.max(image)}', size=25, pad=15)\n    axes[1].set_title(f'Augmented Image\\nMean: {np.mean(image_augmented):.4f} - Std: {np.std(image_augmented):.4f} - Min: {np.min(image_augmented)} - Max: {np.max(image_augmented)}', size=25, pad=15)\n\n    fig.suptitle(\n        f'''\n        Image ID {metadata[\"id\"]} - {metadata[\"organ\"]} - {metadata[\"data_source\"]} - {metadata[\"age\"]} - {metadata[\"sex\"]}\n        Raw Image Shape: {metadata[\"img_height\"]}x{metadata[\"img_width\"]} - Pixel Size: {metadata[\"pixel_size\"]}µm - Tissue Thickness: {metadata[\"tissue_thickness\"]}µm\n        Augmented Image Shape: {metadata[\"img_height\"]}x{metadata[\"img_width\"]} - Pixel Size: {imaging_measurements[\"hubmap\"][\"pixel_size\"][metadata[\"organ\"]]}µm - Tissue Thickness: {imaging_measurements[\"hubmap\"][\"tissue_thickness\"][metadata[\"organ\"]]}µm\n        ''',\n        fontsize=50,\n        y=1.05\n    )\n\n    if path is None:\n        plt.show()\n    else:\n        plt.savefig(path)\n        plt.close(fig)","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:11:14.459646Z","iopub.execute_input":"2022-09-15T04:11:14.460171Z","iopub.status.idle":"2022-09-15T04:11:14.474330Z","shell.execute_reply.started":"2022-09-15T04:11:14.460116Z","shell.execute_reply":"2022-09-15T04:11:14.472885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_augmentation(\n        image=image,\n        image_augmented=hub_image,\n        metadata=row.to_dict()\n    )","metadata":{"execution":{"iopub.status.busy":"2022-09-15T04:12:27.393570Z","iopub.execute_input":"2022-09-15T04:12:27.394035Z","iopub.status.idle":"2022-09-15T04:12:30.361663Z","shell.execute_reply.started":"2022-09-15T04:12:27.394001Z","shell.execute_reply":"2022-09-15T04:12:30.360058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}