diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py index 4cf1ac42..6eaf7dee 100644 --- a/scikits/image/transform/__init__.py +++ b/scikits/image/transform/__init__.py @@ -1 +1,2 @@ from hough_transform import * +from finite_radon_transform import * diff --git a/scikits/image/transform/finite_radon_transform.py b/scikits/image/transform/finite_radon_transform.py new file mode 100755 index 00000000..d2a69df8 --- /dev/null +++ b/scikits/image/transform/finite_radon_transform.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +:author: Gary Ruben, 2009 +:license: modified BSD +""" + +__all__ = ["frt"] +__docformat__ = "restructuredtext en" + +import numpy as np +from numpy import roll, newaxis + +def frt2(a): + """Compute the 2-dimensional finite radon transform (FRT) for an n x n + integer array. + + Parameters + ---------- + a : array_like + A 2-D square n x n integer array. + + Returns + ------- + FRT : 2-D ndarray + Finite Radon Transform array of (n+1) x n integer coefficients. + + See Also + -------- + ifrt2 : The two-dimensional inverse FRT. + + Notes + ----- + The FRT has a unique inverse iff n is prime. [FRT] + The idea for this algorithm is due to Vlad Negnevitski. + + Examples + -------- + + Generate a test image: + Use a prime number for the array dimensions + + >>> SIZE = 59 + >>> img = np.tri(SIZE, dtype=np.int32) + + Apply the Finite Radon Transform: + + >>> f = frt2(img) + + Plot the results: + + >>> import matplotlib.pyplot as plt + >>> plt.imshow(f, interpolation='nearest', cmap=plt.cm.gray) + >>> plt.xlabel('Angle') + >>> plt.ylabel('Translation') + >>> plt.show() + + References + ---------- + .. [FRT] A. Kingston and I. Svalbe, "Projective transforms on periodic + discrete image arrays," in P. Hawkes (Ed), Advances in Imaging + and Electron Physics, 139 (2006) + + """ + if a.ndim != 2 or a.shape[0] != a.shape[1]: + raise ValueError("Input must be a square, 2-D array") + + ai = a.copy() + n = ai.shape[0] + f = np.empty((n+1, n), np.uint32) + f[0] = ai.sum(axis=0) + for m in range(1, n): + # Roll the pth row of ai left by p places + for row in range(1, n): + ai[row] = roll(ai[row], -row) + f[m] = ai.sum(axis=0) + f[n] = ai.sum(axis=1) + return f + + +def ifrt2(a): + """Compute the 2-dimensional inverse finite radon transform (iFRT) for + an (n+1) x n integer array. + + Parameters + ---------- + a : array_like + A 2-D (n+1) row x n column integer array. + + Returns + ------- + iFRT : 2-D n x n ndarray + Inverse Finite Radon Transform array of n x n integer coefficients. + + See Also + -------- + frt2 : The two-dimensional FRT + + Notes + ----- + The FRT has a unique inverse iff n is prime. [FRT] + The idea for this algorithm is due to Vlad Negnevitski. + + Examples + -------- + + >>> SIZE = 59 + >>> img = np.tri(SIZE, dtype=np.int32) + + Apply the Finite Radon Transform: + + >>> f = frt2(img) + + Apply the Inverse Finite Radon Transform to recover the input + >>> fi = ifrt2(f) + + Check that it's identical to the original + + >>> assert len(np.nonzero(img-fi)[0]) == 0 + + References + ---------- + .. [FRT] A. Kingston and I. Svalbe, "Projective transforms on periodic + discrete image arrays," in P. Hawkes (Ed), Advances in Imaging + and Electron Physics, 139 (2006) + + """ + if a.ndim != 2 or a.shape[0] != a.shape[1]+1: + raise ValueError("Input must be an (n+1) row x n column, 2-D array") + + ai = a.copy()[:-1] + n = ai.shape[1] + f = np.empty((n, n), np.uint32) + f[0] = ai.sum(axis=0) + for m in range(1, n): + # Rolls the pth row of ai right by p places. + for row in range(1, ai.shape[0]): + ai[row] = roll(ai[row], row) + f[m] = ai.sum(axis=0) + f += a[-1][newaxis].T + f = (f - ai[0].sum())/n + return f diff --git a/scikits/image/transform/tests/test_finite_radon_transform.py b/scikits/image/transform/tests/test_finite_radon_transform.py new file mode 100755 index 00000000..822e2e26 --- /dev/null +++ b/scikits/image/transform/tests/test_finite_radon_transform.py @@ -0,0 +1,18 @@ +import numpy as np +from numpy.testing import * + +from scikits.image.transform import * + +def test_frt(): + SIZE = 59 + try: + import sympy.ntheory as sn + assert sn.isprime(SIZE) == True + except ImportError: + sn = False + + # Generate a test image + L = np.tri(SIZE, dtype=np.int32) + np.tri(SIZE, dtype=np.int32)[::-1] + f = frt2(L) + fi = ifrt2(f) + assert len(np.nonzero(L-fi)[0]) == 0