Computational cost of Voronoi image interpolation
Thomas.Dane
Thu, 01/05/2012 - 06:35 am
1. regularly-spaced 2D matrix of pixels to XYZ waves
2. apply equations to get the true values or scaling for the X&Y waves - it is at this point that I generate irregularly spaced 2D data.
3. concatenated the XYZ waves to a triplet wave
4. use Voronoi image interpolation to generate the regularly spaced output image
The problem I have is that my data consist of 2048x2048 pixels. In writing my code for this operation I have been working with a segment of data to test which is approx. 250 x 300 pixels. The Interpolation step here takes too long based on how many images I will have to process and this is only with ~80k pixels rather than 4M pixels for a full image!
The Interpolate I have used is as follows:
ImageInterpolate/Dest=$matrixName/S={(xmin),(dx),(xmax),(ymin),(dy),(ymax)} Voronoi, tripletwhere the min/max values are obtained from ImageStats on the triplet and dx=(xmax-xmin)/(cols-1). After getting this code to work I am now scared to try it on a full data image! Is there any way to speed this up or an alternative Interpolation method? The Kriging method is the only other where the input is a triplet XYZ but this actually explicitly states that it is computationally expensive.
Additionally in the wavemetrics macros for doing this I noticed the PFTL=(pftl) flag where pftl = 1e-5, what does this mean? It is one of a few secret flags I have found that aren't listed in the manual!
Many thanks,
Tom
P.S. In case it is important I am running Igor on an iMac with Intel core duo 3.06 GHz & 4Gb of RAM
Voronoi interpolation performance depends on two factors. First is the number of XYZ triplets (N) and second is the density of the output interpolation. Since the triangulation step is O(N^2) you'd want to do what you can to reduce this number (see below). The final density of the interpolation in IP6 is still single threaded but you could compromise (depending on your data) by making it less dense and follow up with simple bi-linear interpolation.
I would start looking for details at your step (2). If the perturbation of the X, Y locations is not too large it would make sense to split your 2048x2048 into smaller rectangular patches (say 128x128), run the calculation (2) on them and follow with the interpolation. At the end you can combine the output patches using ImageTransform or similar.
I can't think of any general methods for interpolating such data. Kriging involves solving a matrix of NxN which would be (2048x2048)^2 and so probably beyond the bounds of 32 bit LAPACK. Loess is another option but it usually does not offer better computational efficiency than Voronoi.
If you want us to look further into this feel free to send us an experiment containing the data (and in particular the computations involved in your step 2) to support@wavemetrics.com so we can be in a better position to advise.
A.G.
WaveMetrics, Inc.
January 5, 2012 at 10:47 am - Permalink
The flag is designed to help with some pathologies and should not be relevant to this problem. The full documentation (which will appear with the next release) is:
/PFTL=tol use this flag with Voronoi interpolation when you want to control the rectangular neighborhood about each datum position XY where the algorithm returns the original Z-value. The tolerance value 0<=tol <1 is tested separately in X and Y (i.e., it is not a geometric distance) when both X and Y are normalized to the range [0,1]. By default tol=1e-15.
A.G.
WaveMetrics, Inc.
January 5, 2012 at 11:16 am - Permalink
As a separate line of thought, I spoke to a someone who also does this kind of task but using python. There is a series of functions within the matlibplot library for numerical operations. See the "matplotlib.mlab.griddata" function about half way down this page http://matplotlib.sourceforge.net/api/mlab_api.html . The person I spoke to uses this function to interpolate after preparing the regular grid in advance. The two functions that cab be used for the interpolation are called "matplotlib.delaunay" and " mpl_tookits.natgrid". I bring this up because this is how they do it for multiple images and I know that the processing does not take a very long time. I wonder if anyone reading this who understands python would be able to have a look at the code (see the link) and see if it would be possible to implement those algorithms in Igor.
Tom
January 6, 2012 at 04:56 am - Permalink
I am not familiar with the details of this library but what you are describing sounds equivalent to using /STW and /PTW flags with ImageInterpolate Voronoi.
A.G.
WaveMetrics, Inc.
January 6, 2012 at 11:33 am - Permalink
I have since gained a tip from someone else who does this processing - their solution was essentially to work backwards. In this way I take the limits of the data to generate the output grid which has the axes scaled with the units I want my data to have (i.e. x axis = Qxy; y axis = Qz in reciprocal angstroms). I then perform a pixel by pixel operation on each pixel in the output grid to find the corresponding pixel of the input image. The main body of the functional code is shown below:
Variable Qz, Qxy, alpha_f, tth_f, x_pix, y_pix, x_pnt, y_pnt
For ( i=0 ; i<DimSize(TransformOutput,0) ; i+=1 )
For ( ii=0 ; ii<DimSize(TransformOutput,1) ; ii+=1 )
// Q-values taken from output grid
Qz = DimOffset(TransformOutput,1) + ii*DimDelta(TransformOutput,1)
Qxy = DimOffset(TransformOutput,0) + i*DimDelta(TransformOutput,0)
// Q-vaules calculated as scattering angles
alpha_f = asin((Qz/waveVector)-(sin(alpha_i)*pi/180))*180/pi
tth_f = acos(((cos(alpha_i*pi/180)^2+cos(alpha_f*pi/180)^2-(Qxy/waveVector)^2))/(2*cos(alpha_i*pi/180)*cos(alpha_f*pi/180)))*180/pi
alpha_f = alpha_f + alpha_i*cos(tth_f*pi/180)
// Corresponding pixel location of scattering angles
x_pix = sign(Qxy)*round((tan(tth_f*pi/180))*LSD)
y_pix = round((tan(alpha_f*pi/180))*(sqrt(x_pix^2 + LSD^2)))
// If Q-values are not recorded by detector i.e. as a consequence of planar detector intersecting
// Ewald sphere then NaN is returned - Do not look for NaN pixel in original image:
if ( numtype(x_pix) != 2 && numtype(y_pix) != 2 )
x_pnt = (x_pix - DimOffset(InputMatrix, 0))/DimDelta(InputMatrix,0)
y_pnt = (y_pix - DimOffset(InputMatrix, 1))/DimDelta(InputMatrix,1)
TransformOutput[i][ii] = InputMatrix[x_pnt][y_pnt]
endif
EndFor
See also the attached images where graph0 is uncorrected and grap0_corr is the corrected image. At the moment this takes approximately 15-20s per image. I will eventually have to perform this on many, many images and so would ideally like to perform this a little more rapidly. Any suggestions for speeding up the processing time?
Cheers,
Tom
February 27, 2012 at 02:50 am - Permalink
Here are some hints (in arbitrary order):
1) Compute values only once if possible. For example,
sin(alpha_i)*pi/180is calculated repeatedly, but as far I can tell it always evaluates to the same value. Just dovariable mySinAlpha = sin(alpha_i)*pi/180and use that in the code.2) Consider using multithreading since the iterations of your loops appear to be parallel.
However, in your case I think you'll get a very significant speedup by doing...
3) Try to reduce explicit looping as much as possible, striving to replace them with calls to operations that will handle multiple datapoints 'at once'. MatrixOP/ImageTransform/WaveTransform are your friends. These built-in operations will execute much faster.
For example, you could precompute Qz and Qxy as 2048x2048 matrices instead of calculating updated values every iteration. Then, once you have that you can calculate all alpha_f in a single step using
MatrixOP /FREE M_AlphaF = asin(M_Qz / waveVector) - some_constant_factorAnd try to do that for the other calculations well.
The downside to this approach is that you'll consume more memory. However, a 2048x2048 matrix of doubles works out to about 33 MB, so you can have a couple of those around without problems.
Note that I didn't consider the algorithm itself in detail, I assumed that the code you provided is correct and efficient.
February 27, 2012 at 04:26 am - Permalink
What you describe is an approach that I use in various places in IGOR. It is fairly standard in Image Processing, e.g., in image rotation algorithms. If you are going to use this approach you might as well take advantage of a few tricks.
1. Eliminate all calls to DimOffset() from your loops. Execute them once outside the loops and store them in variables. Do the same to DimSize() and DimDelta(), even the calls in the for() statement.
2. Eliminate all pi/180. Define a constant that has this value once outside the loops.
3. Qxy should be factored out of the inner loop -- it does not depend on ii.
4. sin(alpha_i)*pi/180 should be evaluated outside the loops as are cos(alpha_i*pi/180)^2 and tan(alpha_f*pi/180), LSD^2, etc.
AG
February 27, 2012 at 10:45 am - Permalink
February 27, 2012 at 01:13 pm - Permalink
Re: 741
I did have a quick go a this approach but performing calculations on the full 2048x2048 data files also seemed to take quite a while. I will have a more thorough play around with this though.
I have another question about this: As the last step involves finding the pixel in the original image by rounding the calculated pixel value, I am worried that I may lose some data points from the original image. Is this likely? It would not be too much of a problem if only a few were not sampled as the the resolution of the image is far greater than the resolution of diffraction rings/spots I have observed.
Correct!
Thanks again.
March 1, 2012 at 01:27 am - Permalink
This is handled differently depending on the particular operation in which it is used. For example, in image rotation I use bilinear interpolation to determine the source pixel. In ImageInterpolate Resample the user can specify the sampling method via the /FUNC flag. If you try some examples you may be able to determine which method provides satisfactory results in your application. Clearly, if you have enough time on your hands, the spline method is better but then again, if your computation time is not an issue than the voronoi approach is superior to everything else.
A.G.
WaveMetrics, Inc.
March 1, 2012 at 10:12 am - Permalink