From 4351d31cd0a94bdb64c68d34645a42a5f6102e6d Mon Sep 17 00:00:00 2001 From: Ziv Yaniv Date: Mon, 8 Apr 2024 11:36:45 -0400 Subject: [PATCH] Adding Laplacian pyramid based blending code. Adding code and example of Laplacian pyramid based image blending. The classic Burt and Adelson approach. --- Python/05_Results_Visualization.ipynb | 145 ++++++++++++++++++++++++-- tests/additional_dictionary.txt | 4 + 2 files changed, 143 insertions(+), 6 deletions(-) diff --git a/Python/05_Results_Visualization.ipynb b/Python/05_Results_Visualization.ipynb index 61785957..bdb74d2a 100644 --- a/Python/05_Results_Visualization.ipynb +++ b/Python/05_Results_Visualization.ipynb @@ -6,12 +6,12 @@ "source": [ "# Visualization of Segmentation and Registration Results \n", "\n", - "In this notebook we illustrate various ways one can display the results of segmentation and registration algorithms so that they can be easily incorporated into a manuscript. For interactive data exploration we recommend using dedicated programs (e.g. 3D slicer). \n", + "In this notebook we illustrate various ways one can display the results of segmentation and registration algorithms so that they can be easily incorporated into a manuscript or presentation. For interactive data exploration we recommend using dedicated programs (e.g. 3D slicer). \n", "\n", "Two key points to remember when working with bio-medical images:\n", "\n", - "1. Most often images have a high dynamic range. Thus, to write them to file in a format appropriate for use in a manuscript we will need to map the intensities to a low dynamic range (e.g. [0,255]). In SimpleITK this is readily done with the [IntensityWindowingImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1IntensityWindowingImageFilter.html).\n", - "2. Images may have non-isotropic spacing between pixels. The file formats appropriate for use in a manuscript (e.g. png, jpg) assume isotropic pixel spacing. This requires that we resample the image before writing to disk. The function `make_isotropic` in the code cell bellow resolves this issue. \n", + "1. Most often images have a high dynamic range. Thus, to write them to file in a format appropriate for use in a manuscript or presentation we will need to map the intensities to a low dynamic range (e.g. [0,255]). In SimpleITK this is readily done with the [IntensityWindowingImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1IntensityWindowingImageFilter.html).\n", + "2. Images may have non-isotropic spacing between pixels. The file formats appropriate for use in manuscripts and presentations (e.g. png, jpg) assume isotropic pixel spacing. This requires that we resample the image before writing to disk. The function `make_isotropic` in the code cell bellow resolves this issue. \n", "\n", "The following filters and their procedural counterparts are useful for various image creation tasks, as illustrated in this notebook:\n", " * [CastImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1CastImageFilter.html)\n", @@ -177,7 +177,7 @@ "## Combining two images\n", "\n", "There are a variety of ways we can overlay two (partially) overlapping images onto each other. The common approaches include:\n", - "1. Use of alpha blending.\n", + "1. Use of blending, alpha or Laplacian pyramid based.\n", "2. Use of a checkerboard pattern with the pixel values in adjacent squares/boxes taken from each of the images.\n", "3. When the pixel values are scalars (gray scale images), combine the two images in different channels, resulting in a color image.\n", "\n", @@ -234,9 +234,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Alpha blending\n", + "### Blending - alpha blending and Laplacian pyramids\n", "\n", - "Alpha blending combines the pixels from the two images as follows:\n", + "**Alpha blending** combines the pixels from the two images as follows:\n", "$$\n", "I_{output} = \\alpha I_1 + (1-\\alpha)I_2,\\;\\;\\; \\alpha \\in[0.0,1.0]\n", "$$\n", @@ -395,6 +395,139 @@ " )" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Laplacian pyramid** based blending results in smooth transitions between the constituent images. This classic approach was described in P. J. Burt, E. H. Adelson, \"A Multiresolution Spline With Application to Image Mosaics\", ACM Transactions on Graphics, 2(4):217-236, 1983. \n", + "\n", + "Using this form of blending will most often yield more visually pleasing results than alpha blending. It should be noted though that in the context of medical image analysis, this may not always be what you want. If there is an alignment error we possibly want it to be clearly visible and not smoothed out by the visualization approach." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def multiresolution_blend(\n", + " image1, image2, mask2=None, smoothing_variances=1, level_num=4\n", + "):\n", + " \"\"\"\n", + " Classic multi-resolution Laplacian pyramid based blending due to Burt and Adelson:\n", + " P. J. Burt, E. H. Adelson, \"A Multiresolution Spline With Application to Image\n", + " Mosaics\", ACM Transactions on Graphics, 2(4):217-236, 1983.\n", + " \"\"\"\n", + " components_per_pixel = image1.GetNumberOfComponentsPerPixel()\n", + " if components_per_pixel > 1:\n", + " image1 = sitk.Cast(image1, sitk.sitkVectorFloat32)\n", + " image2 = sitk.Cast(image2, sitk.sitkVectorFloat32)\n", + " else:\n", + " image1 = sitk.Cast(image1, sitk.sitkFloat32)\n", + " image2 = sitk.Cast(image2, sitk.sitkFloat32)\n", + "\n", + " if not mask2:\n", + " mask2 = sitk.Image(image2.GetSize(), sitk.sitkFloat32) + 1\n", + " mask2.CopyInformation(image2)\n", + " else:\n", + " mask2 = sitk.Cast(mask2, sitk.sitkFloat32)\n", + "\n", + " # Create Laplacian pyramids for images and Gaussian pyramid for mask\n", + " laplacian_pyramid_image1 = create_laplacian_pyramid(\n", + " image1, smoothing_variances, level_num\n", + " )\n", + " laplacian_pyramid_image2 = create_laplacian_pyramid(\n", + " image2, smoothing_variances, level_num\n", + " )\n", + " gaussian_pyramid_mask2 = [mask2]\n", + " for l in range(level_num):\n", + " gaussian_pyramid_mask2.append(\n", + " reduce(gaussian_pyramid_mask2[-1], smoothing_variances)\n", + " )\n", + " # Combine the Laplacian pyramid data using the mask and collapse the combined Laplacian pyramid\n", + " res = (1 - gaussian_pyramid_mask2[-1]) * laplacian_pyramid_image1[\n", + " -1\n", + " ] + gaussian_pyramid_mask2[-1] * laplacian_pyramid_image2[-1]\n", + " for l1, l2, g2 in zip(\n", + " laplacian_pyramid_image1[-2::-1],\n", + " laplacian_pyramid_image2[-2::-1],\n", + " gaussian_pyramid_mask2[-2::-1],\n", + " ):\n", + " res = sitk.Resample(res, l1) + (1 - g2) * l1 + g2 * l2\n", + " return res\n", + "\n", + "\n", + "def create_laplacian_pyramid(image, smoothing_variances, level_num=4):\n", + " gaussian_pyramid = [image]\n", + " for l in range(level_num):\n", + " gaussian_pyramid.append(reduce(gaussian_pyramid[-1], smoothing_variances))\n", + "\n", + " laplacian_pyramid = [\n", + " gaussian_pyramid[l]\n", + " - sitk.Resample(gaussian_pyramid[l + 1], gaussian_pyramid[l])\n", + " for l in range(level_num)\n", + " ]\n", + " laplacian_pyramid.append(gaussian_pyramid[-1])\n", + " return laplacian_pyramid\n", + "\n", + "\n", + "def reduce(image, smoothing_variances):\n", + " if np.isscalar(smoothing_variances):\n", + " smoothing_variances = [smoothing_variances] * image.GetDimension()\n", + " smoothed_image = sitk.DiscreteGaussian(\n", + " image,\n", + " variance=smoothing_variances,\n", + " maximumKernelWidth=32,\n", + " maximumError=0.01,\n", + " useImageSpacing=True,\n", + " )\n", + "\n", + " original_spacing = image.GetSpacing()\n", + " original_size = image.GetSize()\n", + " new_size = [int(sz / 2 + 0.5) for sz in original_size]\n", + " new_spacing = [\n", + " ((original_sz - 1) * original_spc) / (new_sz - 1)\n", + " for original_sz, original_spc, new_sz in zip(\n", + " original_size, original_spacing, new_size\n", + " )\n", + " ]\n", + " return sitk.Resample(\n", + " smoothed_image,\n", + " new_size,\n", + " sitk.Transform(),\n", + " sitk.sitkLinear,\n", + " image.GetOrigin(),\n", + " new_spacing,\n", + " image.GetDirection(),\n", + " 0.0,\n", + " image.GetPixelID(),\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now create 3D images using both the alpha-blending and Laplacian pyramids approaches and visually compare the results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gui.MultiImageDisplay(\n", + " image_list=[\n", + " alpha_blend(img1_255, img2_255, mask2=msk2),\n", + " multiresolution_blend(img1_255, img2_255, msk2),\n", + " ],\n", + " title_list=[\"alpha_blend_mask2\", \"multiresolution_blend_mask2\"],\n", + " figure_size=(9, 3),\n", + " shared_slider=True,\n", + ");" + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/tests/additional_dictionary.txt b/tests/additional_dictionary.txt index 1603e4f8..438cc5be 100644 --- a/tests/additional_dictionary.txt +++ b/tests/additional_dictionary.txt @@ -1,7 +1,9 @@ +ACM ANTSNeighborhoodCorrelation API Acknowledgements AddTransform +Adelson Affine AffineTransform Annulus @@ -164,6 +166,7 @@ LabelOverlayImageFilter LabelShapeStatisticsImageFilter LabelToRGBImageFilter LandmarkBasedTransformInitializer +Laplacian LaplacianSegmentation Lasser Léon @@ -192,6 +195,7 @@ MeanSquares MetaDataDictionaryArrayUpdateOn MetaImageIO MetricEvaluate +Multiresolution NamedTemporaryFile Nanometers Narayanan