{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":45867,"databundleVersionId":6924515,"sourceType":"competition"},{"sourceId":6912790,"sourceType":"datasetVersion","datasetId":3969953,"isSourceIdPinned":true},{"sourceId":6984590,"sourceType":"datasetVersion","datasetId":4014175}],"dockerImageVersionId":30626,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# TMA generator from WSI\nThis notebook allows generating multiple TMA from WSI slides. It picks some random ellipses inside WSI according to provided tumor mask or generated OTSU mask or even randomly when no mask is available.\n\nAbout magnification: UBC WSI are at x20 and TMA (only 25 instances) are at x40. We're generating TMA at x20 based on WSI with no down scaling.\n\nTMA simulator is based on Albumentations. It crops a random ellipse in a given tile, add noise on contours and apply a background color with some random.\n\nTile selection has some constrains such as:\n- Drop tile with high unicolor level (including black)\n- Including some tissue or tumor based on a minimum ratio\n- Squared tile with side among (1482, 1568, 1694). See SplitConfig.\n\nIf some black color (border) is captured then it's replaced by white.\n\nThen the following stain/color augmentations are applied randomly on TMA:\n- Vahadane (with some TMA references)\n- Macenko (with generic reference)\n- Reinhard (with some TMA references)\n\nTMA generation can run in parallel across  multiple CPUs. Limit is driven by memory available and image size.\n\n<b>Generated TMA examples displayed at the end of this notebook.<b>","metadata":{}},{"cell_type":"markdown","source":"### Packages\nInstall PyVips, StainTools/Spams and TorchStain packages","metadata":{}},{"cell_type":"code","source":"# intall the deb packages\n!yes | dpkg -i --force-depends /kaggle/input/pyvips-gpu-2-2-2/linux_packages/archives/*.deb\n# install the python wrapper\n!pip install pyvips -f /kaggle/input/pyvips-gpu-2-2-2/python_packages/ --no-index","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-01-01T14:41:10.424518Z","iopub.execute_input":"2024-01-01T14:41:10.424978Z","iopub.status.idle":"2024-01-01T14:42:26.367993Z","shell.execute_reply.started":"2024-01-01T14:41:10.424944Z","shell.execute_reply":"2024-01-01T14:42:26.366569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# https://github.com/getspams/spams-python\n!git clone https://github.com/getspams/spams-python","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:42:26.370151Z","iopub.execute_input":"2024-01-01T14:42:26.370528Z","iopub.status.idle":"2024-01-01T14:42:28.866858Z","shell.execute_reply.started":"2024-01-01T14:42:26.370495Z","shell.execute_reply":"2024-01-01T14:42:28.865416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cd spams-python","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:42:28.868814Z","iopub.execute_input":"2024-01-01T14:42:28.869412Z","iopub.status.idle":"2024-01-01T14:42:28.879201Z","shell.execute_reply.started":"2024-01-01T14:42:28.869359Z","shell.execute_reply":"2024-01-01T14:42:28.877986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!sed -i 's/np.bool/np.bool_/g' spams/spams.py","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:42:28.882755Z","iopub.execute_input":"2024-01-01T14:42:28.883242Z","iopub.status.idle":"2024-01-01T14:42:29.996211Z","shell.execute_reply.started":"2024-01-01T14:42:28.883197Z","shell.execute_reply":"2024-01-01T14:42:29.994401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -e .","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:42:29.998234Z","iopub.execute_input":"2024-01-01T14:42:29.99917Z","iopub.status.idle":"2024-01-01T14:44:15.81495Z","shell.execute_reply.started":"2024-01-01T14:42:29.99912Z","shell.execute_reply":"2024-01-01T14:44:15.813814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install staintools\n!pip install torchstain","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:15.81757Z","iopub.execute_input":"2024-01-01T14:44:15.817977Z","iopub.status.idle":"2024-01-01T14:44:48.757186Z","shell.execute_reply.started":"2024-01-01T14:44:15.817937Z","shell.execute_reply":"2024-01-01T14:44:48.755562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\npd.set_option('display.max_colwidth', 180)\nimport glob, os, time, random, gc\nos.environ['VIPS_CONCURRENCY'] = '4'\nos.environ['VIPS_DISC_THRESHOLD'] = '15gb'\n\nfrom tqdm.notebook import tqdm\nimport PIL\nfrom PIL import Image\nImage.MAX_IMAGE_PIXELS = None\nimport cv2\nimport pyvips\nimport joblib\nfrom concurrent.futures import ProcessPoolExecutor\nimport seaborn as sns\nfrom matplotlib import pyplot as plt\nfrom mpl_toolkits.axes_grid1 import ImageGrid\nimport shutil\nfrom torchvision import transforms\nimport torch\n\nimport albumentations as A\nimport staintools\nimport torchstain\n\nprint('CV2:', cv2.__version__)\nprint('PIL:', PIL.__version__)\nprint('pyvips:', pyvips.__version__)\nprint('joblib:', joblib.__version__)\nprint(\"Pytorch\", torch.__version__)\nprint(\"Albumentations\", A.__version__)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:48.759936Z","iopub.execute_input":"2024-01-01T14:44:48.760342Z","iopub.status.idle":"2024-01-01T14:44:55.617524Z","shell.execute_reply.started":"2024-01-01T14:44:48.760305Z","shell.execute_reply":"2024-01-01T14:44:55.615783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def seed_everything(seed):\n    \"\"\"\n    Seeds basic parameters for reproducibility of results.\n    Args:\n        seed (int): Number of the seed.\n    \"\"\"\n    random.seed(seed)\n    os.environ[\"PYTHONHASHSEED\"] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    # torch.backends.cudnn.deterministic = True\n    # torch.backends.cudnn.benchmark = False\n    \nseed_everything(42)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.620307Z","iopub.execute_input":"2024-01-01T14:44:55.622728Z","iopub.status.idle":"2024-01-01T14:44:55.651065Z","shell.execute_reply.started":"2024-01-01T14:44:55.622673Z","shell.execute_reply":"2024-01-01T14:44:55.648757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HOME = \"/kaggle/input/UBC-OCEAN\"\nDATA_HOME = os.path.join(HOME, \"train_images\")\nDATA_WSI_THUMBNAILS_HOME = os.path.join(HOME, \"train_thumbnails\")\nMASK_HOME = \"/kaggle/input/ubc-ovarian-cancer-competition-supplemental-masks\"\nTMA_FOLDER = \"/kaggle/working/TMA\"\nMIN_CPU, MAX_CPU = 1, 4\nMAX_MEM = 5000","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.65383Z","iopub.execute_input":"2024-01-01T14:44:55.655922Z","iopub.status.idle":"2024-01-01T14:44:55.666981Z","shell.execute_reply.started":"2024-01-01T14:44:55.655861Z","shell.execute_reply":"2024-01-01T14:44:55.66441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### TMA simulation\n- Pick a random ellipse\n- Add noise on contours\n- add background","metadata":{}},{"cell_type":"code","source":"class SimulateTMA(A.DualTransform):\n\n    def __init__(self, std, radius_ratio=(0.9, 1.0), ellipse_ratio=(0.9, 1.0), angle=(-90., 90.), background_color=(-1, -1, -1), background_color_ratio=1.0, noise_level=(0.0, 0.0), black_replacement_color=None, always_apply=False, p=1.0):\n        super(SimulateTMA, self).__init__(always_apply, p)\n        self.std = std\n        self.radius_ratio = radius_ratio\n        self.ellipse_ratio = ellipse_ratio\n        self.background_color = background_color\n        self.background_color_ratio = background_color_ratio\n        self.angle = angle\n        self.noise_level = noise_level\n        self.black_replacement_color = black_replacement_color\n\n    def apply(self, img, **params):\n        height, width = img.shape[:2]\n        # Replace the black regions with the replacement color\n        if self.black_replacement_color is not None:\n            black_mask = np.all(img == [0, 0, 0], axis=-1)\n            img[black_mask] = self.black_replacement_color\n        img_std = np.std(img) if self.std[0] != -1 else 0  # (20, 50)\n        if (self.std[0] == -1) or ((img_std <= self.std[1]) and (img_std >= self.std[0])):\n            # Draw circle\n            x_center = width // 2\n            y_center = height // 2\n            radius_w = int((width//2)*random.uniform(self.radius_ratio[0], self.radius_ratio[1]))  # Random radius\n            radius_h = int(radius_w*random.uniform(self.ellipse_ratio[0], self.ellipse_ratio[1]))  # int((height//2)*random.uniform(self.radius_ratio[0], self.radius_ratio[1]))  # Random radius\n            angle = int(random.uniform(self.angle[0], self.angle[1]))\n            mask = cv2.ellipse(np.zeros_like(img), (x_center, y_center), (radius_w, radius_h), angle, 0, 360, color=(255, 255, 255), thickness=-1)\n            # Add noise to the contour to mimic TMA\n            if self.noise_level[1] > 0:\n                contour = cv2.findContours(mask[:, :, 0], cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[0]\n                contour_with_noise = contour[0] + np.random.randint(-self.noise_level[0], self.noise_level[1], contour[0].shape)\n                int(random.uniform(self.angle[0], self.angle[1]))\n                mask = cv2.drawContours(np.zeros_like(img), [contour_with_noise], -1, (255, 255, 255), -1)\n            # Apply masks\n            inverse_mask = cv2.bitwise_not(mask)\n            # Background color\n            bg_color = self.background_color\n            if self.background_color == (-1, -1, -1):\n                bg_ratio = random.uniform(self.background_color_ratio[0], self.background_color_ratio[1])\n                bg_color = tuple((np.max(img, axis=(0,1))*bg_ratio).astype(np.uint8)) # Auto color\n            color_outside_circle = np.zeros_like(img) # Black image\n            color_outside_circle[:] = bg_color\n            color_outside_circle = cv2.bitwise_and(color_outside_circle, inverse_mask)\n            img = cv2.bitwise_and(img, mask)\n            img = cv2.add(img, color_outside_circle)\n        return img\n\n    def get_transform_init_args_names(self):\n        return (\"std\", \"radius_ratio\", \"ellipse_ratio\", \"angle\", \"background_color\", \"background_color_ratio\", \"noise_level\", \"black_replacement_color\")","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.682634Z","iopub.execute_input":"2024-01-01T14:44:55.687414Z","iopub.status.idle":"2024-01-01T14:44:55.730559Z","shell.execute_reply.started":"2024-01-01T14:44:55.687328Z","shell.execute_reply":"2024-01-01T14:44:55.729242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Stain/Color augmentations\n- vahadane (with some TMA references)\n- macenko (with generic reference)\n- reinhard (with some TMA references)","metadata":{}},{"cell_type":"code","source":"class Stainer(A.DualTransform):\n\n    def __init__(self, ref_images, method, luminosity=True, always_apply=False, p=1.0):\n        super(Stainer, self).__init__(always_apply, p)\n        self.luminosity = luminosity\n        self.method = method\n        self.stain_normalizer = []\n        self.torchstain_T = transforms.Compose([\n            transforms.ToTensor(),\n            transforms.Lambda(lambda x: x * 255)\n        ])\n        if method == 'macenko':\n            stain_normalizer = torchstain.normalizers.MacenkoNormalizer(backend='torch')\n            self.stain_normalizer.append(stain_normalizer)\n        else:\n            for ref_image in ref_images:\n                ref_image = np.array(Image.open(ref_image))\n                if method == 'reinhard':\n                    stain_normalizer = torchstain.normalizers.ReinhardNormalizer(backend='torch')\n                    stain_normalizer.fit(self.torchstain_T(ref_image))\n                    self.stain_normalizer.append(stain_normalizer)\n                else:\n                    ref_image = staintools.LuminosityStandardizer.standardize(ref_image) if self.luminosity == True else ref_image\n                    stain_normalizer = staintools.StainNormalizer(method=method)\n                    stain_normalizer.fit(ref_image)\n                    self.stain_normalizer.append(stain_normalizer)\n\n    def apply(self, img, **params):\n        # Standardize brightness (optional, can improve the tissue mask calculation)\n        if self.luminosity == True:\n            img = staintools.LuminosityStandardizer.standardize(img)\n        if self.method == 'macenko':\n            stain_normalizer = self.stain_normalizer[0]\n            img, _, _ = stain_normalizer.normalize(I=self.torchstain_T(img), stains=False)\n            img = img.contiguous().cpu().numpy().astype(np.uint8)\n        elif self.method == 'reinhard':\n            stain_normalizer = np.random.choice(self.stain_normalizer, 1)[0]\n            img = stain_normalizer.normalize(I=self.torchstain_T(img))\n            img = img.contiguous().cpu().numpy().astype(np.uint8)\n        else:\n            stain_normalizer = np.random.choice(self.stain_normalizer, 1)[0]\n            img = stain_normalizer.transform(img)\n        return img\n\n    def get_transform_init_args_names(self):\n        return (\"ref_images\", \"method\", \"luminosity\")","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.733612Z","iopub.execute_input":"2024-01-01T14:44:55.734739Z","iopub.status.idle":"2024-01-01T14:44:55.75397Z","shell.execute_reply.started":"2024-01-01T14:44:55.734686Z","shell.execute_reply":"2024-01-01T14:44:55.752754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### OTSU mask to capture tissue","metadata":{}},{"cell_type":"code","source":"def get_otsu_mask(img, wsi_scale, mthresh=7, sthresh=20, sthresh_up = 255, use_otsu = True):\n    img_hsv = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)  # Convert to HSV space\n    img_med = cv2.medianBlur(img_hsv[:,:,1], mthresh)  # Apply median blurring\n    # Thresholding\n    if use_otsu:\n        _, img_otsu = cv2.threshold(img_med, sthresh, sthresh_up, cv2.THRESH_OTSU+cv2.THRESH_BINARY)\n    else:\n        _, img_otsu = cv2.threshold(img_med, sthresh, sthresh_up, cv2.THRESH_BINARY)\n    # Morphological closing\n    close = int(32*wsi_scale)\n    if close > 0:\n        kernel = np.ones((close, close), np.uint8)\n        img_otsu = cv2.morphologyEx(img_otsu, cv2.MORPH_CLOSE, kernel)\n    return img_otsu","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.756038Z","iopub.execute_input":"2024-01-01T14:44:55.7569Z","iopub.status.idle":"2024-01-01T14:44:55.771245Z","shell.execute_reply.started":"2024-01-01T14:44:55.756854Z","shell.execute_reply":"2024-01-01T14:44:55.770087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generate thumbnail to visualize where are the TMAs","metadata":{}},{"cell_type":"code","source":"def generate_thumbnail(image, tiles, conf):\n    scale = conf.tma_thumbnail_scale\n    if scale is not None:\n        image = image.resize(scale, kernel=conf.kernel)\n        image = image.numpy()[..., :3]    \n        for t in tiles:\n            x1, y1, x2, y2 = t[4], t[5], t[6], t[7]\n            uid = str(t[-1])\n            image_id = t[1]\n            x1 = int(x1*scale)\n            y1 = int(y1*scale)\n            x2 = int(x2*scale)\n            y2 = int(y2*scale)\n            x_center = int((x1 + x2)/2)\n            y_center = int((y1 + y2)/2)    \n            font = cv2.FONT_HERSHEY_DUPLEX\n            font_scale = 16.0*scale\n            font_thickness = int(30*scale)        \n            uid_textsize = cv2.getTextSize(uid, font, font_scale, font_thickness)[0]\n            \n            # image = cv2.rectangle(image, (x1, y1), (x2, y2), color=(0,255,255), thickness=int(2))\n            image = cv2.circle(image, (x_center, y_center), (y2-y1)//2, color=(0,255,255), thickness=int(2))\n            \n            image = cv2.putText(image, text = uid, org = (x1+2, y1+2 + int(uid_textsize[1]*1.2)), fontFace = font, fontScale = font_scale, color = (0,0,0), thickness = font_thickness)\n            image = cv2.putText(image, text = uid, org = (x1, y1 + int(uid_textsize[1]*1.2)), fontFace = font, fontScale = font_scale, color = (0,255,255), thickness = font_thickness)\n        if conf.tma_folder is not None:\n            img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n            os.makedirs(img_dir, exist_ok=True)\n            map_file = os.path.join(img_dir, \"thumbnail_map.png\")\n            th = Image.fromarray(image)\n            th.save(map_file)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.773305Z","iopub.execute_input":"2024-01-01T14:44:55.774269Z","iopub.status.idle":"2024-01-01T14:44:55.791689Z","shell.execute_reply.started":"2024-01-01T14:44:55.774209Z","shell.execute_reply":"2024-01-01T14:44:55.790463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generate TMAs from WSI\n- Use mask (if available) with tumor ratio\n- Use OTSU mask (if enough memory to manage it) with ratio\n- Use standard deviation criteria if none of mask nor OTSU mask available","metadata":{}},{"cell_type":"code","source":"def tile_single_image(file, conf):\n    seed_everything(42)\n    \n    max_tiles = conf.max_tiles\n    drop_tile_color_ratio = conf.drop_tile_color_ratio\n        \n    tiles = []\n    filepath, is_tma, has_mask = file\n    image_id = int(filepath.split(\"/\")[-1].replace(\".png\", \"\"))\n    \n    # TMA image x40\n    if is_tma:\n        image = pyvips.Image.new_from_file(filepath)\n        if conf.tma_generation:\n            if isinstance(conf.tma_scale, tuple):\n                image = image.thumbnail_image(conf.tma_scale[1], height=conf.tma_scale[0])\n            else:\n                image = image.resize(conf.tma_scale, kernel=conf.kernel) if conf.tma_scale != 1. else image # x40 => x20, x10 ...         \n        else:\n            if isinstance(conf.tma_scale, tuple):\n                image = image.thumbnail_image(conf.tma_scale[1], height=conf.tma_scale[0])\n            else:\n                image = image.resize(conf.tma_scale, kernel=conf.kernel) if conf.tma_scale != 1. else image # x40 => x20, x10 ...\n\n    else:\n        # Open image\n        image = pyvips.Image.new_from_file(filepath) # default access is random and needs more memory/time # access='sequential'            \n        image = image.resize(conf.wsi_scale, kernel=conf.kernel) if conf.wsi_scale != 1. else image # x20 => x10\n\n    image_width = image.width\n    image_height = image.height\n\n    # Resume\n    if conf.tma_folder is not None:\n        img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n        if os.path.exists(img_dir):\n            return tiles\n\n    if is_tma == False:\n        # Pick box\n        crop_side_w = conf.tma_crop[-1] # Biggest crop\n        crop_side_h = crop_side_w\n        if has_mask:\n            # Open mask\n            mask = pyvips.Image.new_from_file(os.path.join(MASK_HOME, \"%s.png\" % image_id))\n            mask = mask.resize(conf.wsi_scale, kernel=pyvips.enums.Kernel.NEAREST) if conf.wsi_scale != 1. else mask # x20 => x10\n            mask_width = mask.width\n            mask_height = mask.height\n            assert(mask_width == image_width)\n            assert(mask_height == image_height)\n\n            # Find ROI with tumor based on mask\n            idxs = [(y, y + crop_side_h, x, x + crop_side_w, 0, 0, 0) for y in range(0, mask_height, crop_side_h) for x in range(0, mask_width, crop_side_w)]\n            random.shuffle(idxs)\n            for uid, (y, y_, x, x_, _, _, _) in enumerate(idxs):                                                \n                # Update crop size randomly\n                crop_side_w_ = np.random.choice(conf.tma_crop, 1)[0]\n                crop_side_h_ = crop_side_w_                \n                x1, y1, w1, h1 = x, y, min(crop_side_w_, mask_width - x), min(crop_side_h_, mask_height - y)\n                if isinstance(mask, np.ndarray):\n                    x1, y1, x2, y2 = x1, y1, x1 + w1, y1 + h1\n                    mask_tumor_tile = mask[y1:y2, x1:x2]\n                else:\n                    mask_tumor_tile = mask.crop(x1, y1, w1, h1).numpy()[..., :3]\n                    mask_tumor_tile = mask_tumor_tile[:, :, 0]\n                tumoralr = (np.sum((mask_tumor_tile != 0).astype(int)))/(crop_side_w_*crop_side_h_)\n                if tumoralr > conf.tma_tumoral_ratio:\n                    tile_image = image.crop(x1, y1, w1, h1).numpy()[..., :3]\n                    mask_black = np.all(tile_image == [0,0,0], axis=2) # Black to white\n                    tile_image[mask_black] = [255, 255, 255]                                                    \n                     # Make sure the tile is square\n                    if (tile_image.shape[0] != crop_side_h_) or (tile_image.shape[1] != crop_side_w_):\n                        continue\n                    # Random augmentation\n                    if np.random.random() > conf.tma_simulation_prob:\n                        tile_image = conf.tma_simulation(image=tile_image)[\"image\"] if conf.tma_simulation is not None else tile_image\n                    tumor_tile = Image.fromarray(tile_image)\n                    if conf.tma_folder is not None:\n                        img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n                        os.makedirs(img_dir, exist_ok=True)\n                        tile_file = os.path.join(img_dir, os.path.basename(filepath).replace(\".svs\", \".png\").replace(\".png\", \"_%.2f_%d.png\" % (tumoralr, uid)))\n                        tumor_tile.save(tile_file)\n                    tile_height, tile_width = tumor_tile.width, tumor_tile.height                            \n                    # Track TMA selected\n                    x1, y1, x2, y2 = x1, y1, x1 + w1, y1 + h1\n                    tiles.append((filepath, image_id, image_width, image_height, x1, y1, x2, y2, tile_width, tile_height, 1, is_tma, tile_file, \n                                  has_mask, tumoralr, 0, uid))            \n                    if len(tiles) >= conf.tma_max_tiles:\n                        break\n            generate_thumbnail(image, tiles, conf)\n            del mask\n        else:\n            # No mask, pick tile randomly with otsu mask threshold\n            if (conf.otsu_mask_zero_ratio is not None) and (image_width*image_height <= conf.otsu_mask_size_limit):\n                tmp_img = image.numpy()[..., :3]\n                mask = get_otsu_mask(tmp_img, conf.wsi_scale)                        \n                del tmp_img                        \n                gc.collect()\n                mask_width = mask.shape[1]\n                mask_height = mask.shape[0]\n                assert(mask_width == image_width)\n                assert(mask_height == image_height)\n\n                img_dir = os.path.join(conf.tma_folder, str(image_id))\n                os.makedirs(img_dir, exist_ok=True)             \n                (Image.fromarray((mask).astype(np.uint8)).resize((mask_width//4, mask_height//4))).save(os.path.join(img_dir, \"otsu_mask.png\"))\n\n                # Find ROI with tumor based on mask\n                idxs = [(y, y + crop_side_h, x, x + crop_side_w, 0, 0, 0) for y in range(0, mask_height, crop_side_h) for x in range(0, mask_width, crop_side_w)]\n                random.shuffle(idxs)\n                for uid, (y, y_, x, x_, _, _, _) in enumerate(idxs):                                                \n                    # Update crop size randomly\n                    crop_side_w_ = np.random.choice(conf.tma_crop, 1)[0]\n                    crop_side_h_ = crop_side_w_                \n                    x1, y1, w1, h1 = x, y, min(crop_side_w_, mask_width - x), min(crop_side_h_, mask_height - y)\n                    x1, y1, x2, y2 = x1, y1, x1 + w1, y1 + h1                            \n                    tile_mask = mask[y1:y2, x1:x2]\n                    tumoralr = (tile_mask == 0).sum()/(tile_mask.shape[1]*tile_mask.shape[0])\n                    if tumoralr >= conf.otsu_mask_zero_ratio:\n                        continue\n                    tile_image = image.crop(x1, y1, w1, h1).numpy()[..., :3]\n                    mask_black = np.all(tile_image == [0,0,0], axis=2) # Black to white\n                    tile_image[mask_black] = [255, 255, 255]                                \n                     # Make sure the tile is square\n                    if (tile_image.shape[0] != crop_side_h_) or (tile_image.shape[1] != crop_side_w_):\n                        continue\n                    # Ignore uniform title\n                    pix, pixcnt = np.unique(np.array(Image.fromarray(tile_image).convert(\"L\")), return_counts=True)\n                    pixcnt = (pixcnt/(tile_image.shape[0]*tile_image.shape[1]))\n                    if (pixcnt >= conf.drop_tile_color_ratio).any():\n                        continue                 \n                    # Random augmentation\n                    if np.random.random() > conf.tma_simulation_prob:\n                        tile_image = conf.tma_simulation(image=tile_image)[\"image\"] if conf.tma_simulation is not None else tile_image\n                    tumor_tile = Image.fromarray(tile_image)\n                    if conf.tma_folder is not None:\n                        img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n                        os.makedirs(img_dir, exist_ok=True)\n                        tile_file = os.path.join(img_dir, os.path.basename(filepath).replace(\".png\", \"_%d.png\" % uid))\n                        tumor_tile.save(tile_file)\n                    tile_height, tile_width = tumor_tile.width, tumor_tile.height                            \n                    tiles.append((filepath, image_id, image_width, image_height, x1, y1, x2, y2, tile_width, tile_height, 1, is_tma, tile_file, \n                                  has_mask, tumoralr, 0, uid))            \n                    if len(tiles) >= conf.tma_max_tiles:\n                        break\n                generate_thumbnail(image, tiles, conf) \n                del mask\n            # No mask, pick tile randomly with std in a range and not single color\n            else:\n                idxs = [(y, y + crop_side_h, x, x + crop_side_w, 0, 0, 0) for y in range(0, image_height, crop_side_h) for x in range(0, image_width, crop_side_w)]\n                random.shuffle(idxs)\n                for uid, (y, y_, x, x_, _, _, _) in enumerate(idxs):                        \n                    # Update crop size randomly\n                    crop_side_w_ = np.random.choice(conf.tma_crop, 1)[0]\n                    crop_side_h_ = crop_side_w_                        \n                    x1, y1, w1, h1 = x, y, min(crop_side_w_, image_width - x), min(crop_side_h_, image_height - y)\n                    image_tile = image.crop(x1, y1, w1, h1).numpy()[..., :3]\n                    mask_black = np.all(image_tile == [0,0,0], axis=2) # Black to white\n                    image_tile[mask_black] = [255, 255, 255]                           \n                    # Make sure the tile is square\n                    if (image_tile.shape[0] != crop_side_h_) or (image_tile.shape[1] != crop_side_w_):\n                        continue                        \n                    # Ignore uniform title\n                    pix, pixcnt = np.unique(np.array(Image.fromarray(image_tile).convert(\"L\")), return_counts=True)\n                    pixcnt = (pixcnt/(image_tile.shape[0]*image_tile.shape[1]))\n                    if (pixcnt >= conf.drop_tile_color_ratio).any():\n                        continue\n                    image_std = np.std(image_tile) if conf.std is not None else None\n                    if (image_std == None) or ((image_std > conf.std[0]) and (image_std <= conf.std[1])):\n                        # Random augmentation\n                        if np.random.random() > conf.tma_simulation_prob:\n                            image_tile = conf.tma_simulation(image=image_tile)[\"image\"] if conf.tma_simulation is not None else image_tile                      \n                        if conf.tma_folder is not None:\n                            img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n                            os.makedirs(img_dir, exist_ok=True)\n                            tile_file = os.path.join(img_dir, os.path.basename(filepath).replace(\".png\", \"_%d.png\" % uid))\n                            tile = Image.fromarray(image_tile)\n                            tile.save(tile_file)\n                        tile_height, tile_width = tile.width, tile.height                            \n                        # Track TMA selected\n                        x1, y1, x2, y2 = x1, y1, x1 + w1, y1 + h1\n                        tiles.append((filepath, image_id, image_width, image_height, x1, y1, x2, y2, tile_width, tile_height, 1, is_tma, tile_file, \n                                      has_mask, 0, 0, uid))            \n                        if len(tiles) >= conf.tma_max_tiles:\n                            break\n                generate_thumbnail(image, tiles, conf)      \n    else:\n        # TMA\n        if conf.tma_folder is not None:\n            img_dir = os.path.join(conf.tma_folder, \"%s\"%image_id)\n            os.makedirs(img_dir, exist_ok=True)\n            tile_file = os.path.join(img_dir, os.path.basename(filepath).replace(\".png\", \"_%d.png\" % 0))\n            (Image.fromarray(image.numpy()[..., :3])).save(tile_file)       \n        x1, y1, x2, y2 = 0, 0, image_width-1, image_height-1\n        tile_width, tile_height = image_width, image_height\n        tumoralr = 1.0\n        tiles.append((filepath, image_id, image_width, image_height, x1, y1, x2, y2, tile_width, tile_height, 1, is_tma, tile_file, \n                      has_mask, tumoralr, 0, 0))                \n    del image\n    gc.collect()                        \n    return tiles","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.795241Z","iopub.execute_input":"2024-01-01T14:44:55.79606Z","iopub.status.idle":"2024-01-01T14:44:55.865434Z","shell.execute_reply.started":"2024-01-01T14:44:55.796015Z","shell.execute_reply":"2024-01-01T14:44:55.864156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Files and configuration\nThe more TMA references the longer initialization.","metadata":{}},{"cell_type":"code","source":"# List WSI files, thumbnails and masks\nfiles = glob.glob(os.path.join(DATA_HOME, \"*.png\"))\nthumbnails_files = glob.glob(os.path.join(DATA_WSI_THUMBNAILS_HOME, \"*.png\"))\nmasks_files = glob.glob(os.path.join(MASK_HOME, \"*.png\"))\nfiles_pd = pd.DataFrame(files, columns=[\"file\"])\nfiles_pd[\"image_id\"] = files_pd[\"file\"].apply(lambda x: int(x.split(\"/\")[-1].replace(\".png\", \"\")))\n# Thumbnails are only for WSI so we know which images are TMA\nthumbnails_files_pd = pd.DataFrame(thumbnails_files, columns=[\"file\"])\nthumbnails_files_pd[\"image_id\"] = thumbnails_files_pd[\"file\"].apply(lambda x: int(x.split(\"/\")[-1].replace(\"_thumbnail.png\", \"\")))\ntma_list = thumbnails_files_pd[\"image_id\"].unique()\nmasks_files_pd = pd.DataFrame(masks_files, columns=[\"file\"])\nmasks_files_pd[\"image_id\"] = masks_files_pd[\"file\"].apply(lambda x: int(x.split(\"/\")[-1].replace(\".png\", \"\")))\nmasks_list = masks_files_pd[\"image_id\"].unique()\nfiles_pd[\"is_tma\"] = False\nfiles_pd.loc[(~files_pd[\"image_id\"].isin(tma_list)), \"is_tma\"] = True\nfiles_pd[\"has_mask\"] = False\nfiles_pd.loc[(files_pd[\"image_id\"].isin(masks_list)), \"has_mask\"] = True\nprint(\"WSI:\", files_pd[files_pd[\"is_tma\"] == False].shape)\nprint(\"TMA:\", files_pd[files_pd[\"is_tma\"] == True].shape)\nprint(\"MASK:\", files_pd[files_pd[\"has_mask\"] == True].shape)\ntma_pd = files_pd[(files_pd[\"is_tma\"] == True)].reset_index(drop=True)\nfiles = files_pd[[\"file\", \"is_tma\", \"has_mask\"]].values\nprint(files.shape, files[0:10])","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:55.866993Z","iopub.execute_input":"2024-01-01T14:44:55.867803Z","iopub.status.idle":"2024-01-01T14:44:56.128783Z","shell.execute_reply.started":"2024-01-01T14:44:55.867762Z","shell.execute_reply":"2024-01-01T14:44:56.127631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# TMA references for stain augmentation\ntma_images = tma_pd[\"file\"].unique()\ntma_pd.head()","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:56.130234Z","iopub.execute_input":"2024-01-01T14:44:56.130778Z","iopub.status.idle":"2024-01-01T14:44:56.152027Z","shell.execute_reply.started":"2024-01-01T14:44:56.130744Z","shell.execute_reply":"2024-01-01T14:44:56.151099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def tma_augmentation(p=1.0):\n    return A.Compose([\n        SimulateTMA((-1, -1), radius_ratio=(0.6, 1.0), ellipse_ratio=(0.85, 1.15), angle=(-90., 90.), background_color=(-1, -1, -1), background_color_ratio=(0.80, 1.0), noise_level=(20./5, 100./5), black_replacement_color=None, p=1.0, always_apply=True),\n        A.OneOf([\n            Stainer(ref_images=tma_images, method='vahadane', luminosity=True, p=0.34),\n            Stainer(ref_images=None, method='macenko', luminosity=False, p=0.33),\n            Stainer(ref_images=tma_images, method='reinhard', luminosity=False, p=0.33),\n        ], p=0.60),        \n    ], p=p)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:56.153602Z","iopub.execute_input":"2024-01-01T14:44:56.154254Z","iopub.status.idle":"2024-01-01T14:44:56.162457Z","shell.execute_reply.started":"2024-01-01T14:44:56.154219Z","shell.execute_reply":"2024-01-01T14:44:56.161195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SplitConfig:\n\n    # Tiles\n    max_tiles = 500\n    std = (20, 80)\n\n    # Resize interpolation\n    resize_interpolation=Image.LANCZOS\n    kernel = pyvips.enums.Kernel.LANCZOS3\n    \n    wsi_scale = 1.0 # Keep x20\n    tma_scale = 0.5 # x40 => x20\n    tma_crop = [1482, 1568, 1694]\n    tma_folder = TMA_FOLDER\n    tma_tumoral_ratio = 0.75\n    tma_max_tiles = 500\n    tma_simulation = tma_augmentation(p=1.0)\n    tma_simulation_prob = 0.0 # 0.25\n    drop_tile_color_ratio = 0.30 # 0.50\n    otsu_mask_size_limit = 40000*40000 # To prevent OOM for large WSI\n    otsu_mask_zero_ratio = 0.80\n    tma_thumbnail_scale = 0.10\n\nsplit_config = SplitConfig","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:44:56.164164Z","iopub.execute_input":"2024-01-01T14:44:56.164837Z","iopub.status.idle":"2024-01-01T14:55:36.166174Z","shell.execute_reply.started":"2024-01-01T14:44:56.164803Z","shell.execute_reply":"2024-01-01T14:55:36.164793Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sanity check\n- Image size vs mask size\n- Maximum TMAs per WSI","metadata":{}},{"cell_type":"code","source":"def probe_single_image(file, conf):\n    tiles = []\n    h, w = conf.tma_crop[0], conf.tma_crop[0]\n    Image.MAX_IMAGE_PIXELS = None        \n    # Open image\n    filepath, is_tma, has_mask = file\n    image_id = int(filepath.split(\"/\")[-1].replace(\".png\", \"\"))    \n    image = pyvips.Image.new_from_file(filepath)    \n    image_width = image.width\n    image_height = image.height\n    if has_mask:\n        mask_filepath = os.path.join(MASK_HOME, \"%s.png\"%image_id)\n        if os.path.exists(mask_filepath):\n            mask = Image.open(mask_filepath)            \n            mask_width = mask.width\n            mask_height = mask.height\n            assert(image_width == mask_width)\n            assert(image_height == mask_height)\n        else:\n            raise Exception(\"Expected mask not found %s\" % mask_filepath)\n    # Tile image\n    idxs = [(y, y + h, x, x + w) for y in range(0, image_height, h) for x in range(0, image_width, w)]    \n    tiles.append((filepath, image_id, image_width, image_height, len(idxs), is_tma, has_mask))\n    del image\n    return tiles\n\nimages_info = joblib.Parallel(n_jobs=1)(joblib.delayed(probe_single_image)(file, split_config) for file in files)\nimages_probe = []\nfor c in images_info:\n    images_probe.extend(c)\nimages_probe_pd = pd.DataFrame(images_probe, columns=[\"file\", \"image_id\", \"width\", \"height\", \"tiles\", \"is_tma\", \"has_mask\"])\nimages_probe_pd[\"surface\"] = images_probe_pd[\"width\"] * images_probe_pd[\"height\"] / 1000000\n# Create groups for adaptive jobs\nimages_probe_pd[\"jobs\"] = images_probe_pd[\"surface\"].apply(lambda x: max(MIN_CPU, min(MAX_CPU, np.floor(MAX_MEM/x)))).astype(np.int16)\n# Keep non-TMA images\nimages_probe_pd = images_probe_pd.sort_values([\"surface\"], ascending=[True])\nimages_probe_pd = images_probe_pd[images_probe_pd[\"is_tma\"] == False].reset_index(drop=True)\nimages_probe_pd","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:55:36.16844Z","iopub.execute_input":"2024-01-01T14:55:36.168987Z","iopub.status.idle":"2024-01-01T14:55:42.530249Z","shell.execute_reply.started":"2024-01-01T14:55:36.168941Z","shell.execute_reply":"2024-01-01T14:55:42.529019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Limit to a few WSI\nMAX_WSI = 16\nimages_groups = pd.concat([images_probe_pd[images_probe_pd[\"has_mask\"] == True].head(MAX_WSI//2), images_probe_pd[images_probe_pd[\"has_mask\"] == False].head(MAX_WSI//2)], ignore_index=True)\nimages_groups","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:55:42.531825Z","iopub.execute_input":"2024-01-01T14:55:42.532206Z","iopub.status.idle":"2024-01-01T14:55:42.555401Z","shell.execute_reply.started":"2024-01-01T14:55:42.532172Z","shell.execute_reply":"2024-01-01T14:55:42.554355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"images_groups = images_groups.groupby(\"jobs\")[[\"file\", \"is_tma\", \"has_mask\"]].apply(lambda x: list(map(tuple,x.values))).reset_index().sort_values([\"jobs\"], ascending=[False])\nimages_groups.rename(columns={0: \"file\"}, inplace=True)\n\ntiles = []\n# Run by batch, each in parallel with CPUs depending on acceptable memory limit\nfor idx, row in images_groups.iterrows():\n    jobs = int(row[\"jobs\"])\n    images_list = row[\"file\"]\n    images_tiled = joblib.Parallel(n_jobs=jobs)(joblib.delayed(tile_single_image)(file, split_config) for file in tqdm(images_list, total=len(images_list)))\n    for c in images_tiled:\n        tiles.extend(c)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T14:55:42.557151Z","iopub.execute_input":"2024-01-01T14:55:42.557888Z","iopub.status.idle":"2024-01-01T15:08:18.405051Z","shell.execute_reply.started":"2024-01-01T14:55:42.557852Z","shell.execute_reply":"2024-01-01T15:08:18.402592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### TMA visualization","metadata":{}},{"cell_type":"code","source":"for filepath, is_tma, has_mask in images_list:\n    try:\n        image_id = int(filepath.split(\"/\")[-1].replace(\".png\", \"\")) \n        files = glob.glob(os.path.join(TMA_FOLDER, str(image_id), \"*.png\"))\n        fig, ax = plt.subplots(1, 1, figsize=(32, 20))\n        if has_mask:\n            d = ax.imshow(Image.open(os.path.join(MASK_HOME, \"%d.png\" % image_id)), interpolation='none')\n            d = ax.set_title(\"Tumor mask - We keep red only\")\n        elif os.path.join(TMA_FOLDER, str(image_id),'otsu_mask.png') in files:\n            d = ax.imshow(Image.open(os.path.join(TMA_FOLDER, str(image_id),\"otsu_mask.png\")), interpolation='none')\n            d = ax.set_title(\"OTSU mask\")\n        else:\n            d = ax.imshow(Image.open(os.path.join(DATA_WSI_THUMBNAILS_HOME,\"%d_thumbnail.png\"%image_id)))\n        d = plt.show()\n\n        fig, ax = plt.subplots(1, 1, figsize=(32, 32))\n        if has_mask:\n            d = ax.imshow(Image.open(os.path.join(TMA_FOLDER, str(image_id),\"thumbnail_map.png\")))\n            d = ax.set_title(\"Tiles selected for TMA - %d\"%image_id)\n        elif os.path.join(TMA_FOLDER, str(image_id),'otsu_mask.png') in files:\n            d = ax.imshow(Image.open(os.path.join(TMA_FOLDER, str(image_id),\"thumbnail_map.png\")))\n            d = ax.set_title(\"Tiles selected for TMA - %d\"%image_id)\n        else:\n            d = ax.imshow(Image.open(os.path.join(TMA_FOLDER, str(image_id),\"thumbnail_map.png\")))\n            d = ax.set_title(\"Tiles selected for TMA - %d\"%image_id)\n        d = plt.show()    \n\n        files = [f for f in files  if \"otsu_mask.png\" not in f]\n        files = [f for f in files  if \"thumbnail_map.png\" not in f]\n\n        T = 4 # 8\n        chunks = len(files)//T\n        chunks = 1 if chunks ==0 else chunks\n        files = files[0:T*chunks]\n        for j, tmas in enumerate(np.array_split(files, chunks)):\n            fig, ax = plt.subplots(1, T, figsize=(32, 20))\n            for i, tma in enumerate(tmas):\n                d = ax[i].imshow(Image.open(tma))\n                d = ax[i].set_title(tma.split(\"/\")[-1])\n            plt.show()\n            if j > 6:\n                break\n    except Exception as ex:\n        print(ex)","metadata":{"execution":{"iopub.status.busy":"2024-01-01T15:08:18.409125Z","iopub.execute_input":"2024-01-01T15:08:18.409842Z","iopub.status.idle":"2024-01-01T15:22:13.852853Z","shell.execute_reply.started":"2024-01-01T15:08:18.409778Z","shell.execute_reply":"2024-01-01T15:22:13.851549Z"},"trusted":true},"execution_count":null,"outputs":[]}]}