{"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":"## Improving performance with L1/Hessian denoising\n**Brett Olsen, March 2023**\n\n### Introduction\n\nI'm going to demonstrate a potentially useful denoising tool for improving performance of ink detection results using any other approach that produces a continuous probability output.\nThe general approach is to exploit known properties of the ink distribution, in particular that:\n1. The ink is sparse; most regions of the papyrus will not contain ink.\n2. The ink is continuous; a single pixel is much more likely to contain ink if it is next to another ink-containing pixel.\n\nSuppose we have a noisy array $\\boldsymbol{A}$ that contains our initial estimates of where ink is present.\nAssuming it is a reasonable estimate, we can approximate it as a convolution of the truth set $\\boldsymbol{T}$ with an error distribution $\\epsilon$.\nWe would like to make a new estimate $\\boldsymbol{B}$ that minimizes this error term.\nWith no additional information, we would seek to find\n$$\n\\underset{\\boldsymbol{B}}{\\operatorname{argmin}} \\frac{1}{2}\\sum_{i=0}^N(\\boldsymbol{A}[i] - \\boldsymbol{B}[i])^2.\n$$\nThat is, find the minimum least-square fit of $\\boldsymbol{B}$ to $\\boldsymbol{A}$.\nThis is, of course, just $\\boldsymbol{A}$.\nWe can add additional terms to the minimization to account for our additional priors about the ink distribution.\nThe sparsity condition can be accounted for with a scaled L1-norm regularization term:  we simultaneously wish to minimize the least squared distance between $\\boldsymbol{A}$ and $\\boldsymbol{B}$ as well as the L1-norm of the denoised signal $\\boldsymbol{B}$:\n$$\n\\underset{\\boldsymbol{B}}{\\operatorname{argmin}} \\frac{1}{2}\\sum_{i=0}^N(\\boldsymbol{A}[i] - \\boldsymbol{B}[i])^2 + \\lambda \\sum_{i=0}^N | \\boldsymbol{B}[i] |.\n$$","metadata":{}},{"cell_type":"markdown","source":"With the single regularization term, this is analytically solvable (proof is left as an exercise for the reader) as the so-called soft-thresholding function:\n$$\n\\boldsymbol{B}[i] = \n\\begin{cases}\n    |\\boldsymbol{A}[i]| \\leq \\lambda & 0 \\\\\n    \\boldsymbol{A}[i] > \\lambda & \\boldsymbol{A}[i] - \\lambda \\\\\n    \\boldsymbol{A}[i] < -\\lambda & \\boldsymbol{A}[i] + \\lambda \\\\\n\\end{cases}\n$$","metadata":{}},{"cell_type":"markdown","source":"The continuity condition can be accounted for by adding scaled terms for all the entries in the $2 \\times 2$ Hessian matrix of second-order partial derivatives of $\\boldsymbol{B}$:\n\n$$\n\\underset{\\boldsymbol{B}}{\\operatorname{argmin}} \\frac{1}{2}\\sum_{i=0}^N(\\boldsymbol{A}[i] - \\boldsymbol{B}[i])^2 + \\lambda \\sum_{i=0}^N | \\boldsymbol{B}[i] | + \\zeta \\left [ \\frac{\\partial^2}{\\partial x^2} + \\frac{\\partial^2}{\\partial y^2} + 2 \\frac{\\partial^2}{\\partial x \\partial y} \\right ] \\boldsymbol{B}[i].\n$$\n\nThat is, we are penalizing $\\boldsymbol{B}$ for having strongly varying values.\nUnfortunately, this is no longer analytically solvable, so we will use iterative methods to approach the solution.\nSee references:\n- https://eeweb.engineering.nyu.edu/iselesni/lecture_notes/SALSA/SALSA.pdf\n- https://github.com/WeisongZhao/Sparse-SIM","metadata":{}},{"cell_type":"markdown","source":"### Data generation and scoring\n\nLet us begin by generating some synthetic data.\nWe will take the truth set for one of our fragments and convert it to a continuous probability array by adding scaled Gaussian noise.\nThe level of noise will be adjusted so as to achieve similar performance on the [F0.5 scale](https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection/overview/evaluation) used for scoring as the `ink-id` results, of 0.48.","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\n\nfrom PIL import Image\nimport matplotlib.pyplot as plt\nimport torch\nfrom tqdm import tqdm\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(DEVICE)\n\nINPUT_FOLDER = \"/kaggle/input/vesuvius-challenge-ink-detection/train/1/\"","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:45:38.692852Z","iopub.execute_input":"2023-04-11T00:45:38.693332Z","iopub.status.idle":"2023-04-11T00:45:38.704469Z","shell.execute_reply.started":"2023-04-11T00:45:38.693288Z","shell.execute_reply":"2023-04-11T00:45:38.703213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Attempt to import GPU-accelerated numpy routines\ntry:\n    import cupy as cp\n    xp = cp\nexcept ImportError:\n    xp = np","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:41:05.978321Z","iopub.execute_input":"2023-04-11T00:41:05.981038Z","iopub.status.idle":"2023-04-11T00:41:07.347603Z","shell.execute_reply.started":"2023-04-11T00:41:05.980987Z","shell.execute_reply":"2023-04-11T00:41:07.346364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mask = xp.array(Image.open(os.path.join(INPUT_FOLDER, \"mask.png\")), dtype=bool)\ntruth = xp.array(Image.open(os.path.join(INPUT_FOLDER, \"inklabels.png\")), dtype=bool)\n\nxslice = slice(3000, 4500, None)\nyslice = slice(1000, 2500, None)\n\nmask = mask[xslice, yslice]\ntruth = truth[xslice, yslice]\n\ndef generate_noisy_data(truth, mask, scale=0.25):\n    results = xp.zeros_like(truth, dtype=float)\n    results[truth] = 1\n    results += xp.random.normal(loc=0, scale=scale, size=results.shape)\n    results -= results[mask].min()\n    results /= results[mask].max()\n    results[mask == False] = 0\n    return results","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:17.944714Z","iopub.execute_input":"2023-04-11T00:43:17.945627Z","iopub.status.idle":"2023-04-11T00:43:18.352995Z","shell.execute_reply.started":"2023-04-11T00:43:17.945576Z","shell.execute_reply":"2023-04-11T00:43:18.351981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"noisy = generate_noisy_data(truth, mask, 1.0)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:22.919634Z","iopub.execute_input":"2023-04-11T00:43:22.919997Z","iopub.status.idle":"2023-04-11T00:43:23.579583Z","shell.execute_reply.started":"2023-04-11T00:43:22.919966Z","shell.execute_reply":"2023-04-11T00:43:23.578608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def to_numpy(array):\n    if xp == cp:\n        return array.get()\n    return array","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:28.682382Z","iopub.execute_input":"2023-04-11T00:43:28.683080Z","iopub.status.idle":"2023-04-11T00:43:28.688022Z","shell.execute_reply.started":"2023-04-11T00:43:28.683045Z","shell.execute_reply":"2023-04-11T00:43:28.687012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.subplots(1, 3, sharey=True)\nplt.subplot(1, 3, 1)\nplt.imshow(to_numpy(mask), cmap=\"gray\")\nplt.subplot(1, 3, 2)\nplt.imshow(to_numpy(truth), cmap=\"gray\")\nplt.subplot(1, 3, 3)\nplt.imshow(to_numpy(noisy), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:29.835859Z","iopub.execute_input":"2023-04-11T00:43:29.836942Z","iopub.status.idle":"2023-04-11T00:43:30.688874Z","shell.execute_reply.started":"2023-04-11T00:43:29.836896Z","shell.execute_reply":"2023-04-11T00:43:30.687744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can do a bit of clever array work to rapidly score our guesses.","metadata":{}},{"cell_type":"code","source":"def score_guess(guess, truth, beta=0.5, debug=False):\n    results = np.zeros_like(guess, dtype=int)\n    results[truth] += 1\n    results[guess] += 2\n    fn = xp.count_nonzero(results == 1)\n    fp = xp.count_nonzero(results == 2)\n    tp = xp.count_nonzero(results == 3)\n    precision = tp / (tp + fp)\n    recall = tp / (tp + fn)\n    if debug:\n        print(f\"{tp:,} Correct Pixels\")\n        print(f\"{fp:,} False Positives\")\n        print(f\"{fn:,} False Negatives\")\n    return (1 + beta ** 2) * precision * recall / (precision * beta ** 2 + recall)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:35.697567Z","iopub.execute_input":"2023-04-11T00:43:35.698195Z","iopub.status.idle":"2023-04-11T00:43:35.705035Z","shell.execute_reply.started":"2023-04-11T00:43:35.698155Z","shell.execute_reply":"2023-04-11T00:43:35.703927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look at the scoring of our noisy data as we sweep the calling threshold.","metadata":{}},{"cell_type":"code","source":"thresholds = xp.linspace(0.4, 0.7, 10)\nnoisy_masked = noisy[mask]\ntruth_masked = truth[mask]\nscores = xp.array([score_guess(noisy_masked > t, truth_masked) for t in thresholds])","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:37.578962Z","iopub.execute_input":"2023-04-11T00:43:37.579669Z","iopub.status.idle":"2023-04-11T00:43:37.861887Z","shell.execute_reply.started":"2023-04-11T00:43:37.579631Z","shell.execute_reply":"2023-04-11T00:43:37.860873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nplt.grid()\nplt.plot(to_numpy(thresholds), to_numpy(scores))\nplt.xlabel(\"Signal Threshold\")\nplt.ylabel(\"F0.5 Score\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:39.298337Z","iopub.execute_input":"2023-04-11T00:43:39.299404Z","iopub.status.idle":"2023-04-11T00:43:39.537941Z","shell.execute_reply.started":"2023-04-11T00:43:39.299346Z","shell.execute_reply":"2023-04-11T00:43:39.537064Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that we're peaking out here around a score of 0.45, pretty close to the `ink-id` score on the leaderboard.\nLet's look at the calls we make with the signal threshold of 0.6 here.\n\n","metadata":{}},{"cell_type":"code","source":"best_guess = noisy > 0.55\nbest_guess_masked = best_guess[mask]\nscore_guess(best_guess_masked, truth_masked, debug=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:49:47.781848Z","iopub.execute_input":"2023-04-11T00:49:47.782249Z","iopub.status.idle":"2023-04-11T00:49:47.800426Z","shell.execute_reply.started":"2023-04-11T00:49:47.782214Z","shell.execute_reply":"2023-04-11T00:49:47.799052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So most of our calls are right, but we're missing almost three times as many as we're calling.","metadata":{}},{"cell_type":"code","source":"plt.subplots(1, 2, sharey=True)\nplt.subplot(1, 2, 1)\nplt.imshow(to_numpy(truth), cmap=\"gray\")\nplt.subplot(1, 2, 2)\nplt.imshow(to_numpy(noisy) > 0.6, cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:43:48.361376Z","iopub.execute_input":"2023-04-11T00:43:48.362493Z","iopub.status.idle":"2023-04-11T00:43:49.016586Z","shell.execute_reply.started":"2023-04-11T00:43:48.362445Z","shell.execute_reply":"2023-04-11T00:43:49.015405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can see that what is happening is we're calling many, but not all, pixels inside the inked area, while fewer but not none outside, due to the random noise.\n\nLet's see if we can improve performance on this noisy data.\n\n### Denoising\n\nOK, let's look at the actual denoising code.","metadata":{}},{"cell_type":"code","source":"delta_lookup = {\n    \"xx\": xp.array([[1, -2, 1]], dtype=float),\n    \"yy\": xp.array([[1], [-2], [1]], dtype=float),\n    \"xy\": xp.array([[1, -1], [-1, 1]], dtype=float),\n}\n\ndef operate_derivative(img_shape, pair):\n    assert len(img_shape) == 2\n    delta = delta_lookup[pair]\n    fft = xp.fft.fftn(delta, img_shape)\n    return fft * xp.conj(fft)\n\ndef soft_threshold(vector, threshold):\n    return xp.sign(vector) * xp.maximum(xp.abs(vector) - threshold, 0)\n\ndef back_diff(input_image, dim):\n    assert dim in (0, 1)\n    r, n = xp.shape(input_image)\n    size = xp.array((r, n))\n    position = xp.zeros(2, dtype=int)\n    temp1 = xp.zeros((r+1, n+1), dtype=float)\n    temp2 = xp.zeros((r+1, n+1), dtype=float)\n    \n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] += 1\n    position[dim] += 1\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] -= 1\n    return temp1[0:size[0], 0:size[1]]\n\ndef forward_diff(input_image, dim):\n    assert dim in (0, 1)\n    r, n = xp.shape(input_image)\n    size = xp.array((r, n))\n    position = xp.zeros(2, dtype=int)\n    temp1 = xp.zeros((r+1, n+1), dtype=float)\n    temp2 = xp.zeros((r+1, n+1), dtype=float)\n        \n    size[dim] += 1\n    position[dim] += 1\n\n    temp1[position[0]:size[0], position[1]:size[1]] = input_image\n    temp2[position[0]:size[0], position[1]:size[1]] = input_image\n    \n    size[dim] -= 1\n    temp2[0:size[0], 0:size[1]] = input_image\n    temp1 -= temp2\n    size[dim] += 1\n    return -temp1[position[0]:size[0], position[1]:size[1]]\n\ndef iter_deriv(input_image, b, scale, mu, dim1, dim2):\n    g = back_diff(forward_diff(input_image, dim1), dim2)\n    d = soft_threshold(g + b, 1 / mu)\n    b = b + (g - d)\n    L = scale * back_diff(forward_diff(d - b, dim2), dim1)\n    return L, b\n\ndef iter_xx(*args):\n    return iter_deriv(*args, dim1=1, dim2=1)\n\ndef iter_yy(*args):\n    return iter_deriv(*args, dim1=0, dim2=0)\n\ndef iter_xy(*args):\n    return iter_deriv(*args, dim1=0, dim2=1)\n\ndef iter_sparse(input_image, bsparse, scale, mu):\n    d = soft_threshold(input_image + bsparse, 1 / mu)\n    bsparse = bsparse + (input_image - d)\n    Lsparse = scale * (d - bsparse)\n    return Lsparse, bsparse\n\ndef denoise_image(input_image, iter_num=100, fidelity=150, sparsity_scale=10, continuity_scale=0.5, mu=1):\n    image_size = xp.shape(input_image)\n    #print(\"Initialize denoising\")\n    norm_array = (\n        operate_derivative(image_size, \"xx\") + \n        operate_derivative(image_size, \"yy\") + \n        2 * operate_derivative(image_size, \"xy\")\n    )\n    norm_array += (fidelity / mu) + sparsity_scale ** 2\n    b_arrays = {\n        \"xx\": xp.zeros(image_size, dtype=float),\n        \"yy\": xp.zeros(image_size, dtype=float),\n        \"xy\": xp.zeros(image_size, dtype=float),\n        \"L1\": xp.zeros(image_size, dtype=float),\n    }\n    g_update = xp.multiply(fidelity / mu, input_image)\n    for i in tqdm(range(iter_num), total=iter_num):\n        #print(f\"Starting iteration {i+1}\")\n        g_update = xp.fft.fftn(g_update)\n        if i == 0:\n            g = xp.fft.ifftn(g_update / (fidelity / mu)).real\n        else:\n            g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real\n        g_update = xp.multiply((fidelity / mu), input_image)\n        \n        #print(\"XX update\")\n        L, b_arrays[\"xx\"] = iter_xx(g, b_arrays[\"xx\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"YY update\")\n        L, b_arrays[\"yy\"] = iter_yy(g, b_arrays[\"yy\"], continuity_scale, mu)\n        g_update += L\n        \n        #print(\"XY update\")\n        L, b_arrays[\"xy\"] = iter_xy(g, b_arrays[\"xy\"], 2 * continuity_scale, mu)\n        g_update += L\n        \n        #print(\"L1 update\")\n        L, b_arrays[\"L1\"] = iter_sparse(g, b_arrays[\"L1\"], sparsity_scale, mu)\n        g_update += L\n        \n    g_update = xp.fft.fftn(g_update)\n    g = xp.fft.ifftn(xp.divide(g_update, norm_array)).real\n    \n    g[g < 0] = 0\n    g -= g.min()\n    g /= g.max()\n    return g","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-11T00:46:01.872723Z","iopub.execute_input":"2023-04-11T00:46:01.873233Z","iopub.status.idle":"2023-04-11T00:46:01.903288Z","shell.execute_reply.started":"2023-04-11T00:46:01.873187Z","shell.execute_reply":"2023-04-11T00:46:01.902265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"denoised = denoise_image(noisy, iter_num=250)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:46:08.201903Z","iopub.execute_input":"2023-04-11T00:46:08.202287Z","iopub.status.idle":"2023-04-11T00:46:12.020177Z","shell.execute_reply.started":"2023-04-11T00:46:08.202252Z","shell.execute_reply":"2023-04-11T00:46:12.019120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 10))\nplt.subplots(1, 2, sharey=True)\nplt.subplot(1, 2, 1)\nplt.imshow(to_numpy(noisy), cmap=\"gray\")\nplt.subplot(1, 2, 2)\nplt.imshow(to_numpy(denoised), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:46:28.167713Z","iopub.execute_input":"2023-04-11T00:46:28.168348Z","iopub.status.idle":"2023-04-11T00:46:29.381217Z","shell.execute_reply.started":"2023-04-11T00:46:28.168309Z","shell.execute_reply":"2023-04-11T00:46:29.380298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_guess = noisy > 0.6\nbest_denoised = denoised > 0.6\n\nplt.subplots(1, 3, sharey=True, figsize=(10, 10))\nplt.subplot(1, 3, 1)\nplt.imshow(to_numpy(truth), cmap=\"gray\")\nplt.subplot(1, 3, 2)\nplt.imshow(to_numpy(best_guess), cmap=\"gray\")\nplt.subplot(1, 3, 3)\nplt.imshow(to_numpy(best_denoised), cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:47:14.540212Z","iopub.execute_input":"2023-04-11T00:47:14.540585Z","iopub.status.idle":"2023-04-11T00:47:15.740550Z","shell.execute_reply.started":"2023-04-11T00:47:14.540552Z","shell.execute_reply":"2023-04-11T00:47:15.739567Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"denoised_masked = denoised[mask]\ndenoised_scores = xp.array([score_guess(denoised_masked > t, truth_masked) for t in thresholds])","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:47:31.365910Z","iopub.execute_input":"2023-04-11T00:47:31.366288Z","iopub.status.idle":"2023-04-11T00:47:31.432974Z","shell.execute_reply.started":"2023-04-11T00:47:31.366255Z","shell.execute_reply":"2023-04-11T00:47:31.431988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nplt.grid()\nplt.plot(to_numpy(thresholds), to_numpy(scores), label=\"Original\")\nplt.plot(to_numpy(thresholds), to_numpy(denoised_scores), label=\"Denoised\")\nplt.xlabel(\"Signal Threshold\")\nplt.ylabel(\"F0.5 Score\")\nplt.legend(loc=0)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:48:07.493472Z","iopub.execute_input":"2023-04-11T00:48:07.494381Z","iopub.status.idle":"2023-04-11T00:48:07.735999Z","shell.execute_reply.started":"2023-04-11T00:48:07.494332Z","shell.execute_reply":"2023-04-11T00:48:07.735066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_guess = noisy > 0.53\nbest_guess_masked = best_guess[mask]\nscore_guess(best_guess_masked, truth_masked, debug=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:48:51.189384Z","iopub.execute_input":"2023-04-11T00:48:51.190695Z","iopub.status.idle":"2023-04-11T00:48:51.208452Z","shell.execute_reply.started":"2023-04-11T00:48:51.190656Z","shell.execute_reply":"2023-04-11T00:48:51.207309Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_guess = denoised > 0.53\nbest_guess_masked = best_guess[mask]\nscore_guess(best_guess_masked, truth_masked, debug=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T00:48:55.476570Z","iopub.execute_input":"2023-04-11T00:48:55.477162Z","iopub.status.idle":"2023-04-11T00:48:55.495721Z","shell.execute_reply.started":"2023-04-11T00:48:55.477117Z","shell.execute_reply":"2023-04-11T00:48:55.494766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Results","metadata":{}},{"cell_type":"markdown","source":"Results look promising, at least on the kind of completely random noise I'm modeling here.\nWe're seeing a significant improvement in score just from this denoising application.\nMoreover, we're seeing the improvements across a wide range of thresholds.\nWe can see that the improvements are across the board:  an increased number of true positives, fewer false positives, and fewer false negatives.\n\nCompared to the previous version, we've significantly improved performance by putting most of the work onto the GPU.\nHowever, the limited memory available on the GPU meant we needed to denoise only a smaller section of the image.\nTo denoise the entire region, we'd need to chunk and denoise in pieces, but this is perfectly feasible with the ~100x performance improvement from running on the GPU.\n\nFor future work:\n1. Tuning the parameters for better performance would be useful.  In particular, I think we're probably doing too many iterations and could cut back a bit.  The scale factors for the three components have not been tuned at all, nor has the size of the identity regularization (`fidelity`).\n2. The noise model needs to be tested.  Ideally, I'd like to get ahold of actual outputs from various people's machine learning approaches and see if a denoising pass would improve performance.  I'll ask around on the discord.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}