{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Using Polar and Log-Polar Transformations for Registration\n\nPhase correlation (``registration.phase_cross_correlation``) is an efficient\nmethod for determining translation offset between pairs of similar images.\nHowever this approach relies on a near absence of rotation/scaling differences\nbetween the images, which are typical in real-world examples.\n\nTo recover rotation and scaling differences between two images, we can take\nadvantage of two geometric properties of the log-polar transform and the\ntranslation invariance of the frequency domain. First, rotation in Cartesian\nspace becomes translation along the angular coordinate ($\\theta$) axis\nof log-polar space. Second, scaling in Cartesian space becomes translation\nalong the radial coordinate ($\\rho = \\ln\\sqrt{x^2 + y^2}$) of log-polar\nspace. Finally, differences in translation in the spatial domain do not impact\nmagnitude spectrum in the frequency domain.\n\nIn this series of examples, we build on these concepts to show how the\nlog-polar transform ``transform.warp_polar`` can be used in conjunction with\nphase correlation to recover rotation and scaling differences between two\nimages that also have a translation offset.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recover rotation difference with a polar transform\n\nIn this first example, we consider the simple case of two images that only\ndiffer with respect to rotation around a common center point. By remapping\nthese images into polar space, the rotation difference becomes a simple\ntranslation difference that can be recovered by phase correlation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom skimage import data\nfrom skimage.registration import phase_cross_correlation\nfrom skimage.transform import warp_polar, rotate, rescale\nfrom skimage.util import img_as_float\n\nradius = 705\nangle = 35\nimage = data.retina()\nimage = img_as_float(image)\nrotated = rotate(image, angle)\nimage_polar = warp_polar(image, radius=radius, channel_axis=-1)\nrotated_polar = warp_polar(rotated, radius=radius, channel_axis=-1)\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 8))\nax = axes.ravel()\nax[0].set_title(\"Original\")\nax[0].imshow(image)\nax[1].set_title(\"Rotated\")\nax[1].imshow(rotated)\nax[2].set_title(\"Polar-Transformed Original\")\nax[2].imshow(image_polar)\nax[3].set_title(\"Polar-Transformed Rotated\")\nax[3].imshow(rotated_polar)\nplt.show()\n\nshifts, error, phasediff = phase_cross_correlation(image_polar, rotated_polar)\nprint(f'Expected value for counterclockwise rotation in degrees: '\n f'{angle}')\nprint(f'Recovered value for counterclockwise rotation: '\n f'{shifts[0]}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recover rotation and scaling differences with log-polar transform\n\nIn this second example, the images differ by both rotation and scaling (note\nthe axis tick values). By remapping these images into log-polar space,\nwe can recover rotation as before, and now also scaling, by phase\ncorrelation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# radius must be large enough to capture useful info in larger image\nradius = 1500\nangle = 53.7\nscale = 2.2\nimage = data.retina()\nimage = img_as_float(image)\nrotated = rotate(image, angle)\nrescaled = rescale(rotated, scale, channel_axis=-1)\nimage_polar = warp_polar(image, radius=radius,\n scaling='log', channel_axis=-1)\nrescaled_polar = warp_polar(rescaled, radius=radius,\n scaling='log', channel_axis=-1)\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 8))\nax = axes.ravel()\nax[0].set_title(\"Original\")\nax[0].imshow(image)\nax[1].set_title(\"Rotated and Rescaled\")\nax[1].imshow(rescaled)\nax[2].set_title(\"Log-Polar-Transformed Original\")\nax[2].imshow(image_polar)\nax[3].set_title(\"Log-Polar-Transformed Rotated and Rescaled\")\nax[3].imshow(rescaled_polar)\nplt.show()\n\n# setting `upsample_factor` can increase precision\nshifts, error, phasediff = phase_cross_correlation(image_polar, rescaled_polar,\n upsample_factor=20)\nshiftr, shiftc = shifts[:2]\n\n# Calculate scale factor from translation\nklog = radius / np.log(radius)\nshift_scale = 1 / (np.exp(shiftc / klog))\n\nprint(f'Expected value for cc rotation in degrees: {angle}')\nprint(f'Recovered value for cc rotation: {shiftr}')\nprint()\nprint(f'Expected value for scaling difference: {scale}')\nprint(f'Recovered value for scaling difference: {shift_scale}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Register rotation and scaling on a translated image - Part 1\n\nThe above examples only work when the images to be registered share a\ncenter. However, it is more often the case that there is also a translation\ncomponent to the difference between two images to be registered. One\napproach to register rotation, scaling and translation is to first correct\nfor rotation and scaling, then solve for translation. It is possible to\nresolve rotation and scaling differences for translated images by working on\nthe magnitude spectra of the Fourier-transformed images.\n\nIn this next example, we first show how the above approaches fail when two\nimages differ by rotation, scaling, and translation.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from skimage.color import rgb2gray\nfrom skimage.filters import window, difference_of_gaussians\nfrom scipy.fft import fft2, fftshift\n\nangle = 24\nscale = 1.4\nshiftr = 30\nshiftc = 15\n\nimage = rgb2gray(data.retina())\ntranslated = image[shiftr:, shiftc:]\nrotated = rotate(translated, angle)\nrescaled = rescale(rotated, scale)\nsizer, sizec = image.shape\nrts_image = rescaled[:sizer, :sizec]\n\n# When center is not shared, log-polar transform is not helpful!\nradius = 705\nwarped_image = warp_polar(image, radius=radius, scaling=\"log\")\nwarped_rts = warp_polar(rts_image, radius=radius, scaling=\"log\")\nshifts, error, phasediff = phase_cross_correlation(warped_image, warped_rts,\n upsample_factor=20)\nshiftr, shiftc = shifts[:2]\nklog = radius / np.log(radius)\nshift_scale = 1 / (np.exp(shiftc / klog))\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 8))\nax = axes.ravel()\nax[0].set_title(\"Original Image\")\nax[0].imshow(image, cmap='gray')\nax[1].set_title(\"Modified Image\")\nax[1].imshow(rts_image, cmap='gray')\nax[2].set_title(\"Log-Polar-Transformed Original\")\nax[2].imshow(warped_image)\nax[3].set_title(\"Log-Polar-Transformed Modified\")\nax[3].imshow(warped_rts)\nfig.suptitle('log-polar-based registration fails when no shared center')\nplt.show()\n\nprint(f'Expected value for cc rotation in degrees: {angle}')\nprint(f'Recovered value for cc rotation: {shiftr}')\nprint()\nprint(f'Expected value for scaling difference: {scale}')\nprint(f'Recovered value for scaling difference: {shift_scale}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Register rotation and scaling on a translated image - Part 2\n\nWe next show how rotation and scaling differences, but not translation\ndifferences, are apparent in the frequency magnitude spectra of the images.\nThese differences can be recovered by treating the magnitude spectra as\nimages themselves, and applying the same log-polar + phase correlation\napproach taken above.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# First, band-pass filter both images\nimage = difference_of_gaussians(image, 5, 20)\nrts_image = difference_of_gaussians(rts_image, 5, 20)\n\n# window images\nwimage = image * window('hann', image.shape)\nrts_wimage = rts_image * window('hann', image.shape)\n\n# work with shifted FFT magnitudes\nimage_fs = np.abs(fftshift(fft2(wimage)))\nrts_fs = np.abs(fftshift(fft2(rts_wimage)))\n\n# Create log-polar transformed FFT mag images and register\nshape = image_fs.shape\nradius = shape[0] // 8 # only take lower frequencies\nwarped_image_fs = warp_polar(image_fs, radius=radius, output_shape=shape,\n scaling='log', order=0)\nwarped_rts_fs = warp_polar(rts_fs, radius=radius, output_shape=shape,\n scaling='log', order=0)\n\nwarped_image_fs = warped_image_fs[:shape[0] // 2, :] # only use half of FFT\nwarped_rts_fs = warped_rts_fs[:shape[0] // 2, :]\nshifts, error, phasediff = phase_cross_correlation(warped_image_fs,\n warped_rts_fs,\n upsample_factor=10)\n\n# Use translation parameters to calculate rotation and scaling parameters\nshiftr, shiftc = shifts[:2]\nrecovered_angle = (360 / shape[0]) * shiftr\nklog = shape[1] / np.log(radius)\nshift_scale = np.exp(shiftc / klog)\n\nfig, axes = plt.subplots(2, 2, figsize=(8, 8))\nax = axes.ravel()\nax[0].set_title(\"Original Image FFT\\n(magnitude; zoomed)\")\ncenter = np.array(shape) // 2\nax[0].imshow(image_fs[center[0] - radius:center[0] + radius,\n center[1] - radius:center[1] + radius],\n cmap='magma')\nax[1].set_title(\"Modified Image FFT\\n(magnitude; zoomed)\")\nax[1].imshow(rts_fs[center[0] - radius:center[0] + radius,\n center[1] - radius:center[1] + radius],\n cmap='magma')\nax[2].set_title(\"Log-Polar-Transformed\\nOriginal FFT\")\nax[2].imshow(warped_image_fs, cmap='magma')\nax[3].set_title(\"Log-Polar-Transformed\\nModified FFT\")\nax[3].imshow(warped_rts_fs, cmap='magma')\nfig.suptitle('Working in frequency domain can recover rotation and scaling')\nplt.show()\n\nprint(f'Expected value for cc rotation in degrees: {angle}')\nprint(f'Recovered value for cc rotation: {recovered_angle}')\nprint()\nprint(f'Expected value for scaling difference: {scale}')\nprint(f'Recovered value for scaling difference: {shift_scale}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Some notes on this approach\n\nIt should be noted that this approach relies on a couple of parameters\nthat have to be chosen ahead of time, and for which there are no clearly\noptimal choices:\n\n1. The images should have some degree of bandpass filtering\napplied, particularly to remove high frequencies, and different choices here\nmay impact outcome. The bandpass filter also complicates matters because,\nsince the images to be registered will differ in scale and these scale\ndifferences are unknown, any bandpass filter will necessarily attenuate\ndifferent features between the images. For example, the log-polar transformed\nmagnitude spectra don't really look \"alike\" in the last example here. Yet if\nyou look closely, there are some common patterns in those spectra, and they\ndo end up aligning well by phase correlation as demonstrated.\n\n2. Images must be windowed using windows with circular symmetry, to remove\nthe spectral leakage coming from image borders. There is no clearly optimal\nchoice of window.\n\nFinally, we note that large changes in scale will dramatically alter the\nmagnitude spectra, especially since a big change in scale will usually be\naccompanied by some cropping and loss of information content. The literature\nadvises staying within 1.8-2x scale change [1]_ [2]_. This is fine for most\nbiological imaging applications.\n\n### References\n\n.. [1] B.S. Reddy and B.N. Chatterji. An FFT-based technique for translation,\n rotation and scale- invariant image registration. IEEE Trans. Image\n Processing, 5(8):1266\u20131271, 1996. :DOI:`10.1109/83.506761`\n\n.. [2] Tzimiropoulos, Georgios, and Tania Stathaki. \"Robust FFT-based\n scale-invariant image registration.\" In 4th SEAS DTC Technical\n Conference. 2009. :DOI:`10.1109/TPAMI.2010.107`\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 0
}