speed of iradon transform

[Matthew Newville]

Hi,

Some time ago, I raised Issue #929 on github about the interpolation step in skimage.transforms.radon_transform.iradon being too slow.

After some on-and-off investigation, I found two opportunities to speed this up.    First, numpy's linear interpolation routine interp() was very slow for repeated interpolations of well-ordered input -- this is now fixed in the github master repo for numy, and gives roughly a 3x to 4x speed-up. 

Second, the trigonometric calculations can be cached when doing repeated calls to iradon() for the same geometry (image size, number of angles).  This can give another factor of 1.5 to 2.0x speed-up.    This is PR #1474.

Some preliminary timing results:

Python2.7.8 (linux), with an image shape of (500, 500), and sinogram shape of (500, 361). I used iradon() options of  dict(filter='shepp-logan', interpolation='linear') and tried both circle=True and False. 

The "workspace?" column here indicates whether the iradon_workspace() function introduced in PR1474 was used.

with circle=True:

   numpy      skimage    workspace?   Best, Worst of 5 (s)
 |---------+-----------+------------+---------------------|
   1.9.2      master      N/A            6.085  6.132
   master     master      N/A            1.784, 1.813
   master     PR1474      No             1.788, 1.809
   master     PR1474      Yes            1.103, 1.160
 |---------+-----------+------------+---------------------|

with circle=False:

   numpy      skimage    workspace?   Best, Worst of 5 (s)
 |---------+-----------+------------+---------------------|
   1.9.2      master      N/A            2.736, 2.767
   master     master      N/A            0.820, 0.831
   master     PR1474      No             0.802, 0.814
   master     PR1474      Yes            0.535, 0.540
 |---------+-----------+------------+---------------------|  

That is, the cumulative improvement is about 5x.    Also note that not using the workspace reverts to the older behavior, and that using he workspace for 1 run of iradon() has basically no improvement over not using the workspace.

Oddly, iradon() is now faster than radon() on this machine/image size (radon() takes about 3.5 sec).   I don't understand why that would be.  Anybody understand why radon() is so slow?

Cheers,

--Matt Newville

[Stefan van der Walt]

Hi Matt

On 2015-04-20 04:59:06, Matthew Newville <matt.n...@gmail.com>
wrote:
> Some preliminary timing results:

Thank you very much for your detailed investigation!

> numpy skimage workspace? Best, Worst of 5 (s)
> |---------+-----------+------------+---------------------|
> 1.9.2 master N/A 6.085 6.132 master
> master N/A 1.784, 1.813 master PR1474
> No 1.788, 1.809 master PR1474 Yes
> 1.103, 1.160
> |---------+-----------+------------+---------------------|

It looks like one gets about 1.5x for using the workspace. I am
always careful about added code complexity for a relatively small
gain, but in this case your refactoring *improves* legibility of
the code--so +1 from me.

> Oddly, iradon() is now faster than radon() on this machine/image
> size (radon() takes about 3.5 sec). I don't understand why
> that would be. Anybody understand why radon() is so slow?

I'm afraid we don't implement an optimized version of the forward
radon transform, such as

Brady, "A Fast Discrete Approximation Algorithm for the Radon
Transform", SIAM J. Comput., 27(1), 2006

A fast digital Radon transform--An efficient means for evaluating
the Hough transform WA Götz, HJ Druckmüller - Pattern Recognition,
1996

This overview from William Press of Numerical Recipes fame is
helpful:

http://www.pnas.org/content/103/51/19249.full

Regards
Stéfan