Merge pull request #442 from zivy/laplacianPyramidBlending

Adding laplacian pyramid based blending code.

Python/05_Results_Visualization.ipynb

@@ -6,12 +6,12 @@
    "source": [
     "# Visualization of Segmentation and Registration Results <a href=\"https://mybinder.org/v2/gh/InsightSoftwareConsortium/SimpleITK-Notebooks/master?filepath=Python%2F05_Results_Visualization.ipynb\"><img style=\"float: right;\" src=\"https://mybinder.org/badge_logo.svg\"></a>\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": {},

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