One of my first posts to this blog presented a Python implementation of the log-polar transform. Back then I was just beginning to learn NumPy, and didn't know many of its more useful tricks; in particular, I had no domain of fancy indexing, and did not know that it could be used to avoid the need for time-consuming Python loops. Besides that, I later realized that the algorithm I chose to implement, taken from an article I had read at the time, was unnecessarily complicated.

A couple of months ago I was considering the use of the log-polar transform for another project, and thought I should go back and rewrite my early implementation. I wanted to make use of the NumPy expertise I've earned since then, but also clean up my early blunder and present a simpler algorithm, that could be easily ported to other languages. The results are available below, and also for download in a nice package, with test code and thoroughly commented sources.

The first listing presents a "naive" implementation, which relies on Python loops:

def logpolar_naive(image, i_0, j_0, p_n=None, t_n=None): (i_n, j_n) = image.shape[:2] i_c = max(i_0, i_n - i_0) j_c = max(j_0, j_n - j_0) d_c = (i_c ** 2 + j_c ** 2) ** 0.5 if p_n == None: p_n = int(ceil(d_c)) if t_n == None: t_n = j_n p_s = log(d_c) / p_n t_s = 2.0 * pi / t_n transformed = zeros((p_n, t_n) + image.shape[2:], dtype=image.dtype) for p in range(0, p_n): p_exp = exp(p * p_s) for t in range(0, t_n): t_rad = t * t_s i = int(i_0 + p_exp * sin(t_rad)) j = int(j_0 + p_exp * cos(t_rad)) if 0 <= i < i_n and 0 <= j < j_n: transformed[p, t] = image[i, j] return transformedThe second version constructs a pair of NumPy fancy indices first – cached for later reference – and then uses them to perform the transform:

_transforms = {} def _get_transform(i_0, j_0, i_n, j_n, p_n, t_n, p_s, t_s): transform = _transforms.get((i_0, j_0, i_n, j_n, p_n, t_n)) if transform == None: i_k = [] j_k = [] p_k = [] t_k = [] for p in range(0, p_n): p_exp = exp(p * p_s) for t in range(0, t_n): t_rad = t * t_s i = int(i_0 + p_exp * sin(t_rad)) j = int(j_0 + p_exp * cos(t_rad)) if 0 <= i < i_n and 0 <= j < j_n: i_k.append(i) j_k.append(j) p_k.append(p) t_k.append(t) transform = ((array(p_k), array(t_k)), (array(i_k), array(j_k))) _transforms[i_0, j_0, i_n, j_n, p_n, t_n] = transform return transform def logpolar_fancy(image, i_0, j_0, p_n=None, t_n=None): (i_n, j_n) = image.shape[:2] i_c = max(i_0, i_n - i_0) j_c = max(j_0, j_n - j_0) d_c = (i_c ** 2 + j_c ** 2) ** 0.5 if p_n == None: p_n = int(ceil(d_c)) if t_n == None: t_n = j_n p_s = log(d_c) / p_n t_s = 2.0 * pi / t_n (pt, ij) = _get_transform(i_0, j_0, i_n, j_n, p_n, t_n, p_s, t_s) transformed = zeros((p_n, t_n) + image.shape[2:], dtype=image.dtype) transformed[pt] = image[ij] return transformedBoth implementations are based on a "reverse transform" algorithm: rather than iterate over the Cartesian image calculating where each pixel fits in the transform (which resulted in all sorts of complications), iterate over the

*transform's*cells and calculate which, if any, pixel from the Cartesian image should be put there.

According to my tests, the "fancy" version using a cached pair of indices runs about 10 times faster than the "naive" version – but even when it has to construct the indices before using them, it's still slightly faster.