{"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":"GITHUB_REPO_URL = \"https://github.com/caethan/vesuvius_image.git\"\n!git clone {GITHUB_REPO_URL}\n!git -C /kaggle/working/vesuvius_image/ pull","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-11T02:39:36.810357Z","iopub.execute_input":"2023-04-11T02:39:36.811065Z","iopub.status.idle":"2023-04-11T02:39:40.013147Z","shell.execute_reply.started":"2023-04-11T02:39:36.811024Z","shell.execute_reply":"2023-04-11T02:39:40.011744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nimport os\nclass HiddenPrints:\n    def __enter__(self):\n        self._original_stdout = sys.stdout\n        sys.stdout = open(os.devnull, 'w')\n\n    def __exit__(self, exc_type, exc_val, exc_tb):\n        sys.stdout.close()\n        sys.stdout = self._original_stdout","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:39:40.019706Z","iopub.execute_input":"2023-04-11T02:39:40.020132Z","iopub.status.idle":"2023-04-11T02:39:40.033598Z","shell.execute_reply.started":"2023-04-11T02:39:40.020076Z","shell.execute_reply":"2023-04-11T02:39:40.031296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with HiddenPrints():\n    !mamba env update --name base --file vesuvius_image/conda_env.yaml","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:39:40.037622Z","iopub.execute_input":"2023-04-11T02:39:40.038814Z","iopub.status.idle":"2023-04-11T02:44:10.600467Z","shell.execute_reply.started":"2023-04-11T02:39:40.038764Z","shell.execute_reply":"2023-04-11T02:44:10.599059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.optim as optim\nimport torch.utils.data as data\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:44:10.604311Z","iopub.execute_input":"2023-04-11T02:44:10.604724Z","iopub.status.idle":"2023-04-11T02:44:11.696775Z","shell.execute_reply.started":"2023-04-11T02:44:10.604675Z","shell.execute_reply":"2023-04-11T02:44:11.695616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if \"/kaggle/working/vesuvius_image/\" not in sys.path:\n    sys.path.append(\"/kaggle/working/vesuvius_image/\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:44:11.698647Z","iopub.execute_input":"2023-04-11T02:44:11.699831Z","iopub.status.idle":"2023-04-11T02:44:11.704890Z","shell.execute_reply.started":"2023-04-11T02:44:11.699787Z","shell.execute_reply":"2023-04-11T02:44:11.703713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from vesuvius.image import PapyrusImage\nfrom vesuvius.volume import PapyrusVolume\nfrom vesuvius.utils import Timer\n\nfrom tqdm import tqdm\nimport numpy as np\nfrom typing import Union\nimport scipy\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:44:11.706853Z","iopub.execute_input":"2023-04-11T02:44:11.707776Z","iopub.status.idle":"2023-04-11T02:44:11.884847Z","shell.execute_reply.started":"2023-04-11T02:44:11.707732Z","shell.execute_reply":"2023-04-11T02:44:11.883666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_FOLDER = \"/kaggle/input/vesuvius-challenge-ink-detection/\"\nWORKING_FOLDER = \"/kaggle/working/\"","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:44:11.886200Z","iopub.execute_input":"2023-04-11T02:44:11.886576Z","iopub.status.idle":"2023-04-11T02:44:11.894237Z","shell.execute_reply.started":"2023-04-11T02:44:11.886534Z","shell.execute_reply":"2023-04-11T02:44:11.892837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This cell will generate the zarr files if they don't already exist.\n# Note we're running with \"multiscale=False\" to avoid using more space than we need for the\n# downsampled copies.\nfor subdir in [\"1\"]:\n    train_input_dir = os.path.join(INPUT_FOLDER, \"train\", subdir)\n    train_output_dir = os.path.join(WORKING_FOLDER, f\"train_{subdir}\")\n    os.makedirs(train_output_dir, exist_ok=True)\n    print(train_output_dir)\n    # First the surface volume stack\n    frag_output = os.path.join(train_output_dir, \"frag.volzarr\")\n    if not os.path.exists(frag_output):\n        # TODO: experiment with different chunk sizes here\n        PapyrusVolume.build_from_tiffdir(\n            os.path.join(train_input_dir, \"surface_volume\"), \n            frag_output, \n            multiscale=False,\n        )\n    # The fragment mask\n    mask_output = os.path.join(train_output_dir, \"mask.imagezarr\")\n    if not os.path.exists(mask_output):\n        PapyrusImage.build_from_image(\n            os.path.join(train_input_dir, \"mask.png\"),\n            mask_output,\n            multiscale=False,\n        )\n    # The truth set\n    ink_output = os.path.join(train_output_dir, \"inklabels.imagezarr\")\n    if not os.path.exists(ink_output):\n        PapyrusImage.build_from_image(\n            os.path.join(train_input_dir, \"inklabels.png\"),\n            ink_output,\n            multiscale=False\n        )","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:47:44.102507Z","iopub.execute_input":"2023-04-11T02:47:44.103402Z","iopub.status.idle":"2023-04-11T02:47:44.115400Z","shell.execute_reply.started":"2023-04-11T02:47:44.103363Z","shell.execute_reply":"2023-04-11T02:47:44.114242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BUFFER = 30  # Buffer size in x and y direction\nZ_START = 27 # First slice in the z direction to use\nZ_DIM = 10   # Number of slices in the z direction\n\nlabelzarr = PapyrusImage(f\"{WORKING_FOLDER}/train_1/inklabels.imagezarr\")\nlabel = torch.from_numpy(\n    labelzarr[:,:,0].astype(np.float32)\n).to(DEVICE)\nfragzarr = PapyrusVolume(f\"{WORKING_FOLDER}/train_1/frag.volzarr\")\nimage_stack = torch.from_numpy(\n    # N.B. the raw values are uint16, we need to cast to float and rescale to (0,1)\n    fragzarr[:,:,Z_START:Z_START+Z_DIM,0].astype(np.float32) / np.iinfo(np.uint16).max\n).to(DEVICE)\nmaskzarr = PapyrusImage(f\"{WORKING_FOLDER}/train_1/mask.imagezarr\")\nmask = maskzarr[:,:,0]","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:47:48.408922Z","iopub.execute_input":"2023-04-11T02:47:48.411795Z","iopub.status.idle":"2023-04-11T02:47:59.575066Z","shell.execute_reply.started":"2023-04-11T02:47:48.411735Z","shell.execute_reply":"2023-04-11T02:47:59.573885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with Timer():\n    print(\"Generating pixel lists...\")\n    pixels_list = set()\n    for pixel in zip(*np.where(mask)):\n        if pixel[1] < BUFFER or pixel[1] >= mask.shape[1]-BUFFER or pixel[0] < BUFFER or pixel[0] >= mask.shape[0]-BUFFER:\n            continue # Too close to the edge of the image\n        pixels_list.add(pixel)\n    pixels_list = list(pixels_list)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:47:59.577243Z","iopub.execute_input":"2023-04-11T02:47:59.577635Z","iopub.status.idle":"2023-04-11T02:48:47.488865Z","shell.execute_reply.started":"2023-04-11T02:47:59.577593Z","shell.execute_reply":"2023-04-11T02:48:47.487590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SubvolumeDataset(data.Dataset):\n    def __init__(self, image_stack, label, pixels):\n        self.image_stack = image_stack\n        self.label = label\n        self.pixels = pixels\n    def __len__(self):\n        return len(self.pixels)\n    def __getitem__(self, index):\n        y, x = self.pixels[index]\n        subvolume = self.image_stack[y-BUFFER:y+BUFFER+1, x-BUFFER:x+BUFFER+1, :].view(1, BUFFER*2+1, BUFFER*2+1, Z_DIM)\n        inklabel = self.label[y, x].view(1)\n        return subvolume, inklabel","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:48:47.490276Z","iopub.execute_input":"2023-04-11T02:48:47.491342Z","iopub.status.idle":"2023-04-11T02:48:49.104644Z","shell.execute_reply.started":"2023-04-11T02:48:47.491300Z","shell.execute_reply":"2023-04-11T02:48:49.103472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = nn.Sequential(\n    nn.Conv3d(in_channels=1, out_channels=16, kernel_size=3, stride=1, padding=1),\n    nn.MaxPool3d(kernel_size=2, stride=2),\n    nn.Conv3d(in_channels=16, out_channels=32, kernel_size=3, stride=1, padding=1),\n    nn.MaxPool3d(kernel_size=2, stride=2),\n    nn.Conv3d(in_channels=32, out_channels=64, kernel_size=3, stride=1, padding=1),\n    nn.MaxPool3d(kernel_size=2, stride=2),\n    #Flattens dimensions 1+ into a single dimension, so it is a dim=2 matrix afterwards\n    nn.Flatten(start_dim=1),\n    nn.LazyLinear(out_features=128),\n    nn.ReLU(),\n    nn.LazyLinear(out_features=1),\n    # Convert from (-inf, inf) to (0, 1) probability with a sigmoid\n    nn.Sigmoid()\n).to(DEVICE)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:48:49.108285Z","iopub.execute_input":"2023-04-11T02:48:49.109229Z","iopub.status.idle":"2023-04-11T02:48:49.145728Z","shell.execute_reply.started":"2023-04-11T02:48:49.109186Z","shell.execute_reply":"2023-04-11T02:48:49.144591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Training on fragment 1...\")\nTRAINING_STEPS = 30000\nLEARNING_RATE = 0.03\nBATCH_SIZE = 32\ntrain_dataset = SubvolumeDataset(image_stack, label, pixels_list)\ntrain_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)\ncriterion = nn.BCELoss()\noptimizer = optim.SGD(model.parameters(), lr=LEARNING_RATE)\nscheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=LEARNING_RATE, total_steps=TRAINING_STEPS)\nmodel.train()\nrunning_loss = 0.0\nfor i, (subvolumes, inklabels) in tqdm(enumerate(train_loader), total=TRAINING_STEPS):\n    if i >= TRAINING_STEPS:\n        break\n    optimizer.zero_grad()\n    outputs = model(subvolumes.to(DEVICE))\n    loss = criterion(outputs, inklabels.to(DEVICE))\n    loss.backward()\n    optimizer.step()\n    scheduler.step()\n    running_loss += loss.item()\n    if i % 3000 == 3000-1:\n        print(\"Loss:\", running_loss / 3000)\n        running_loss = 0.0","metadata":{"execution":{"iopub.status.busy":"2023-04-11T02:48:49.147119Z","iopub.execute_input":"2023-04-11T02:48:49.148049Z","iopub.status.idle":"2023-04-11T02:59:45.171766Z","shell.execute_reply.started":"2023-04-11T02:48:49.148008Z","shell.execute_reply":"2023-04-11T02:59:45.170566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is obviously non-kosher, but we're just going to evaluate the trained model on the training set itself to get some ML output data to play with.","metadata":{}},{"cell_type":"code","source":"eval_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=False)\noutput = torch.zeros_like(label).float()\nmodel.eval()\nwith torch.no_grad():\n    for i, (subvolumes, _) in enumerate(tqdm(eval_loader)):\n        for j, value in enumerate(model(subvolumes.to(DEVICE))):\n            output[pixels_list[i*BATCH_SIZE+j]] = value\n\nfig, (ax1, ax2) = plt.subplots(1, 2)\nax1.imshow(output.cpu(), cmap='gray')\nax2.imshow(label.cpu(), cmap='gray')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-11T03:01:58.063819Z","iopub.execute_input":"2023-04-11T03:01:58.064198Z","iopub.status.idle":"2023-04-11T04:48:31.101907Z","shell.execute_reply.started":"2023-04-11T03:01:58.064163Z","shell.execute_reply":"2023-04-11T04:48:31.100844Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And since it took a couple hours to do the inference here, I'm going to save the output array so I can use it again later without spending all this time.","metadata":{}},{"cell_type":"code","source":"np.save(\"/kaggle/working/frag1_results.array\", output.cpu().numpy())","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:50:05.666790Z","iopub.execute_input":"2023-04-11T04:50:05.667497Z","iopub.status.idle":"2023-04-11T04:50:06.004911Z","shell.execute_reply.started":"2023-04-11T04:50:05.667462Z","shell.execute_reply":"2023-04-11T04:50:06.003881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"OK.  Now let's sweep some threshold values and compare the F0.5 score at different thresholds.","metadata":{}},{"cell_type":"code","source":"def dice_coef(preds, targets, beta=0.5, smooth=1e-5):\n    preds = preds.view(-1).float()\n    targets = targets.view(-1).float()\n    \n    y_true_count = targets.sum()\n    ctp = preds[targets == 1].sum()\n    cfp = preds[targets == 0].sum()\n    beta_squared = beta * beta\n    \n    c_precision = ctp / (ctp + cfp + smooth)\n    c_recall = ctp / (y_true_count + smooth)\n    dice = (1 + beta_squared) * (c_precision * c_recall) / (beta_squared * c_precision + c_recall + smooth)\n    return dice","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:51:32.229286Z","iopub.execute_input":"2023-04-11T04:51:32.230359Z","iopub.status.idle":"2023-04-11T04:51:32.239102Z","shell.execute_reply.started":"2023-04-11T04:51:32.230303Z","shell.execute_reply":"2023-04-11T04:51:32.237537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"thresholds = np.linspace(0.2, 0.9, 25)\nscores = np.zeros_like(thresholds)\nfor i, t in enumerate(tqdm(thresholds)):\n    scores[i] = dice_coef(output > t, label)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:52:04.154141Z","iopub.execute_input":"2023-04-11T04:52:04.154505Z","iopub.status.idle":"2023-04-11T04:52:04.357851Z","shell.execute_reply.started":"2023-04-11T04:52:04.154471Z","shell.execute_reply":"2023-04-11T04:52:04.356760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nplt.grid()\nplt.plot(thresholds, scores)\nidx = np.argmax(scores)\nprint(thresholds[idx], scores[idx])","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:53:26.483721Z","iopub.execute_input":"2023-04-11T04:53:26.484161Z","iopub.status.idle":"2023-04-11T04:53:26.715701Z","shell.execute_reply.started":"2023-04-11T04:53:26.484120Z","shell.execute_reply":"2023-04-11T04:53:26.714591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Cool.  So the maximum score is about 0.713 at a threshold of 0.55.\nLet's take a quick look at what the image looks like.","metadata":{}},{"cell_type":"code","source":"plt.imshow(output.cpu() > 0.55, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:54:11.769879Z","iopub.execute_input":"2023-04-11T04:54:11.770471Z","iopub.status.idle":"2023-04-11T04:54:14.156808Z","shell.execute_reply.started":"2023-04-11T04:54:11.770422Z","shell.execute_reply":"2023-04-11T04:54:14.155667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, what I want to do is see how *large* the computed inked areas are for each threshold.\n`sklearn` has some tools that will make this fairly straightforward.","metadata":{}},{"cell_type":"code","source":"from skimage import measure, data\nimport pandas as pd","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:00:19.619526Z","iopub.execute_input":"2023-04-11T05:00:19.619978Z","iopub.status.idle":"2023-04-11T05:00:19.625228Z","shell.execute_reply.started":"2023-04-11T05:00:19.619934Z","shell.execute_reply":"2023-04-11T05:00:19.624256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_labels = measure.label(output.cpu() > 0.55, connectivity=2)\nprint(output_labels.max())","metadata":{"execution":{"iopub.status.busy":"2023-04-11T04:57:31.934320Z","iopub.execute_input":"2023-04-11T04:57:31.935119Z","iopub.status.idle":"2023-04-11T04:57:33.632987Z","shell.execute_reply.started":"2023-04-11T04:57:31.935075Z","shell.execute_reply":"2023-04-11T04:57:33.631716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"OK.  So we found 26,645 different contiguous regions of ink at this threshold.\nLet's see if we can break down some of the properties of these regions.","metadata":{}},{"cell_type":"code","source":"properties = [\n    \"label\",\n    \"area\",\n    \"area_bbox\",\n    \"area_convex\",\n    \"axis_major_length\",\n]\noutput_table = pd.DataFrame(\n    measure.regionprops_table(\n        output_labels, properties=properties\n    )\n)\noutput_table.head()","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:07:07.361363Z","iopub.execute_input":"2023-04-11T05:07:07.361739Z","iopub.status.idle":"2023-04-11T05:07:26.568316Z","shell.execute_reply.started":"2023-04-11T05:07:07.361703Z","shell.execute_reply":"2023-04-11T05:07:26.567118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I've pulled out some properties that sklearn can compute that I thought might be relevant:  `label` is just a range of integer labels for each distinct region, `area` is the number of pixels it covers, `area_bbox` is the area of the bounding box covering that region, `area_convex` is the area of the convex hull covering that region, and `axis_major_length` is the longer length of an ellipse fitted to that region.\n\nLet's look at the distribution of these properties across all my regions.","metadata":{}},{"cell_type":"code","source":"for prop in properties:\n    if prop == \"label\":\n        continue\n    plt.figure(figsize=(8, 4))\n    plt.grid()\n    plt.title(prop)\n    plt.hist(np.log1p(output_table[prop].values) / np.log(10), bins=100);\n    plt.yscale('log')","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:11:49.622532Z","iopub.execute_input":"2023-04-11T05:11:49.623456Z","iopub.status.idle":"2023-04-11T05:11:52.177534Z","shell.execute_reply.started":"2023-04-11T05:11:49.623401Z","shell.execute_reply":"2023-04-11T05:11:52.176371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks like we're more or less dominated by the small stuff here.\nThere's more single-pixel regions than any others!\nI suspect some post-processing will improve things a fair bit.\nLet's sweep over thresholds and cross compare the F0.5 score versus the sizes of the regions.","metadata":{}},{"cell_type":"code","source":"thresholds = np.linspace(0.2, 0.8, 20)\nscores = np.zeros_like(thresholds)\n\noutput_tables = {}\nfor i, t in enumerate(tqdm(thresholds)):\n    scores[i] = dice_coef(output > t, label)\n    output_labels = measure.label(output.cpu() > t, connectivity=2)\n    output_tables[t] = pd.DataFrame(\n        measure.regionprops_table(\n            output_labels, properties=properties\n        )\n    )","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:38:45.615685Z","iopub.execute_input":"2023-04-11T05:38:45.616732Z","iopub.status.idle":"2023-04-11T05:44:53.832284Z","shell.execute_reply.started":"2023-04-11T05:38:45.616680Z","shell.execute_reply":"2023-04-11T05:44:53.831141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_tables[thresholds[-1]]","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:44:53.835714Z","iopub.execute_input":"2023-04-11T05:44:53.836747Z","iopub.status.idle":"2023-04-11T05:44:53.849742Z","shell.execute_reply.started":"2023-04-11T05:44:53.836705Z","shell.execute_reply":"2023-04-11T05:44:53.848713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's look at the number of pixels in regions of size x.\nWe've got a number of very small regions that don't actually account for an enormous number of pixels.\nAdjusting by region size should scale up the larger regions and let us look at differences more profitably.","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:22:30.821924Z","iopub.execute_input":"2023-04-11T05:22:30.822307Z","iopub.status.idle":"2023-04-11T05:22:30.830151Z","shell.execute_reply.started":"2023-04-11T05:22:30.822271Z","shell.execute_reply":"2023-04-11T05:22:30.828992Z"}}},{"cell_type":"code","source":"for t in thresholds:\n    plt.figure(figsize=(10, 1.5))\n    plt.grid()\n    plt.title(f\"threshold {t}\")\n    plt.hist(output_tables[t][\"axis_major_length\"].values, weights=output_tables[t][\"area\"].values, bins=50)\n    plt.xlim((0, 3500))","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:45:26.523124Z","iopub.execute_input":"2023-04-11T05:45:26.523509Z","iopub.status.idle":"2023-04-11T05:45:33.769767Z","shell.execute_reply.started":"2023-04-11T05:45:26.523471Z","shell.execute_reply":"2023-04-11T05:45:33.768623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"OK.  So let's set pixel thresholds for the length of the regions, say 400-1000 pixels.  Then we'll plot what fraction of pixels are in regions of that length for each threshold.","metadata":{}},{"cell_type":"code","source":"MIN_LENGTH = 400\nMAX_LENGTH = 1000\n\nfrac_in_range = []\nfor t in thresholds:\n    results, _ = np.histogram(\n        output_tables[t][\"axis_major_length\"].values,\n        weights=output_tables[t][\"area\"].values,\n        bins=[0, MIN_LENGTH, MAX_LENGTH, 2 * max(output.shape)],\n    )\n    frac = results[1] / results.sum()\n    frac_in_range.append(frac)\n    \nplt.figure()\nplt.grid()\nplt.plot(thresholds, frac_in_range, label=\"Fraction of pixels in target\")\nplt.plot(thresholds, scores, label=\"F0.5 score\")\nplt.legend(loc=0)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:46:16.729983Z","iopub.execute_input":"2023-04-11T05:46:16.730927Z","iopub.status.idle":"2023-04-11T05:46:17.000643Z","shell.execute_reply.started":"2023-04-11T05:46:16.730871Z","shell.execute_reply":"2023-04-11T05:46:16.999631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Neat!  That works pretty well.  Let's tweak the thresholds a bit and see if we can do better.","metadata":{}},{"cell_type":"code","source":"MIN_LENGTH = 300\nMAX_LENGTH = 900\n\nfrac_in_range = []\nfor t in thresholds:\n    results, _ = np.histogram(\n        output_tables[t][\"axis_major_length\"].values,\n        weights=output_tables[t][\"area\"].values,\n        bins=[0, MIN_LENGTH, MAX_LENGTH, 2 * max(output.shape)],\n    )\n    frac = results[1] / results.sum()\n    frac_in_range.append(frac)\n    \nplt.figure()\nplt.grid()\nplt.plot(thresholds, frac_in_range, label=\"Fraction of pixels in target range\")\nplt.plot(thresholds, scores, label=\"F0.5 score\")\nplt.legend(loc=0)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T05:48:42.017223Z","iopub.execute_input":"2023-04-11T05:48:42.017588Z","iopub.status.idle":"2023-04-11T05:48:42.307611Z","shell.execute_reply.started":"2023-04-11T05:48:42.017553Z","shell.execute_reply":"2023-04-11T05:48:42.306623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Not bad!  We're not quite at the maximum score, but we're pretty damn close now.\nThe other thing I wanted to try here was for a given threshold, if we filter out regions that are very small or very large, we can improve our scores a bit.","metadata":{}},{"cell_type":"code","source":"local_output = output.cpu()\n\ndef max_intensity(region_mask):\n    return local_output[region_mask].max()\n\ndef correct_pixels(region_mask):\n    return inklabels[region_mask].sum()\n\noutput_labels = measure.label(local_output > 0.55, connectivity=2)\noutput_tables = pd.DataFrame(\n    measure.regionprops_table(\n        output_labels, properties=properties,\n    )\n)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T06:08:41.511394Z","iopub.execute_input":"2023-04-11T06:08:41.511789Z","iopub.status.idle":"2023-04-11T06:08:56.266980Z","shell.execute_reply.started":"2023-04-11T06:08:41.511753Z","shell.execute_reply":"2023-04-11T06:08:56.265960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"output_tables","metadata":{"execution":{"iopub.status.busy":"2023-04-11T06:12:00.543563Z","iopub.execute_input":"2023-04-11T06:12:00.544585Z","iopub.status.idle":"2023-04-11T06:12:00.559592Z","shell.execute_reply.started":"2023-04-11T06:12:00.544537Z","shell.execute_reply":"2023-04-11T06:12:00.558403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"inklabels = torch.from_numpy(\n    labelzarr[:,:,0].astype(np.float32)\n).to(DEVICE)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MIN_SIZE = 60\nMAX_SIZE = 1500\nlengths = output_tables[\"axis_major_length\"].values\nvalid_labels = output_tables[\"label\"][(MIN_SIZE <= lengths) & (lengths <= MAX_SIZE)].values\nfiltered_result = np.isin(output_labels, valid_labels)\nprint(output_labels[filtered_result].all())\nprint(dice_coef(output > 0.55, inklabels))\nprint(dice_coef(torch.from_numpy(filtered_result.astype(np.float32)).to(DEVICE), inklabels))","metadata":{"execution":{"iopub.status.busy":"2023-04-11T06:28:11.382577Z","iopub.execute_input":"2023-04-11T06:28:11.383304Z","iopub.status.idle":"2023-04-11T06:28:14.519598Z","shell.execute_reply.started":"2023-04-11T06:28:11.383267Z","shell.execute_reply":"2023-04-11T06:28:14.518348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like there's no benefit to filtering out big stuff, just to filtering out small stuff.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}