adds project_polar_to_cartesian function, with tests

this function takes an image that was created with polar coordinates
and remaps it to cartesian coordinates. it's currently fairly limited,
only working on square images, but that's what we need for now

Author: billbrod

TESTS/unitTests.py

@@ -1160,6 +1160,19 @@ def test2(self):
         res = pt.skew(disc, mn, v)
         self.assertTrue(np.absolute(res - mres) <= np.power(10.0,-11))
 
+
+class ProjectPolarTests(unittest.TestCase):
+    def test0(self):
+        # check that it runs, even though here the output will be nonsensical
+        img = pt.synthetic_images.disk(256)
+        proj_img = pt.project_polar_to_cartesian(img)
+    def test1(self):
+        # currently only works for square images
+        img = pt.synthetic_images.disk((256, 512))
+        with self.assertRaises(Exception):
+            pt.project_polar_to_cartesian(img)
+
+
 # class cconv2Tests(unittest.TestCase):
 #     def test0(self):
 #         matPyr = scipy.io.loadmat(op.join(matfiles_path, 'cconv2_0.mat'))

pyrtools/__init__.py

@@ -7,7 +7,7 @@
 from .tools.convolutions import blurDn, blur, upBlur, image_gradient, rconv2
 from .tools.display import imshow, animshow, pyrshow
 from .tools.image_stats import image_compare, image_stats, range, skew, var, entropy
-from .tools.utils import rcosFn, matlab_histo, matlab_round
+from .tools.utils import rcosFn, matlab_histo, matlab_round, project_polar_to_cartesian
 from .tools.compare_matpyrtools import comparePyr, compareRecon
 
 from .version import version as __version__

pyrtools/tools/utils.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy import ndimage
 import warnings
 
 
@@ -119,6 +120,65 @@ def rcosFn(width=1, position=0, values=(0, 1)):
     return (X, Y)
 
 
+def project_polar_to_cartesian(data):
+    """Take a function defined in polar coordinates and project it into Cartesian coordinates
+
+    Inspired by https://pyabel.readthedocs.io/en/latest/_modules/abel/tools/polar.html, which went
+    the other way. Note that we currently don't implement the Cartesian to polar projection, but
+    could do so based on this code fairly simply if it's necessary.
+
+    Currently, this only works for square images and we require that the original image and the
+    reprojected image are the same size. There should be a way to avoid both of these issues, but I
+    can't think of a way to do that right now.
+
+    Parameters
+    ----------
+    data : array_like
+        The 2d array to convert from polar to Cartesian coordinates. We assume the first dimension
+        is the polar radius and the second is the polar angle.
+
+    Returns
+    -------
+    output : np.array
+        The 2d array in Cartesian coordinates.
+
+    """
+    if np.isnan(data).any():
+        data[np.isnan(data)] = 0
+        warnings.warn("project_polar_to_cartesian won't work if there are any NaNs in the array, "
+                      "so we've replaced all NaNs with 0s")
+    nx = data.shape[1]
+    ny = data.shape[0]
+    if nx != ny:
+        raise Exception("There's an occasional bug where we don't wrap the angle correctly if nx "
+                        "and ny aren't equal, so we don't support this for now!")
+
+    max_radius = data.shape[0]
+    x_i = np.linspace(-max_radius, max_radius, nx, endpoint=False)
+    y_i = np.linspace(-max_radius, max_radius, ny, endpoint=False)
+    x_grid, y_grid = np.meshgrid(x_i, y_i)
+    # need to flip the y indices so that negative is at the bottom (to correspond with how we have
+    # the polar angle -- 0 on the right)
+    y_grid = np.flipud(y_grid)
+
+    r = np.sqrt(x_grid**2 + y_grid**2)
+
+    theta = np.arctan2(y_grid, x_grid)
+    # having the angle run from 0 to 2 pi seems to avoid most of the discontinuities
+    theta = np.mod(theta, 2*np.pi)
+    # need to convert from 2pi to pixel values
+    theta *= nx/(2*np.pi)
+
+    r_i, theta_i = r.flatten(), theta.flatten()
+    # map_coordinates requires a 2xn array
+    coords = np.vstack((r_i, theta_i))
+    # we use mode="nearest" to deal with weird discontinuities that may pop up near the theta=0
+    # line
+    zi = ndimage.map_coordinates(data, coords, mode='nearest')
+    output = zi.reshape((ny, nx))
+    return output
+
+
 if __name__ == "__main__":
     X, Y = rcosFn(width=1, position=0, values=(0, 1))