Range Flag for MatrixOP::Correlate

It is not clear what your end goal is. Do you want a single, small average translation vector between the pair of images?

If so, how about creating a difference image (with suitable boundary adjustments), then take the MatrixOP FFT of the difference image. Suitably scaling the dominant low-frequency FFT component should provide the translation vector.

The real goal is SPM drift analysis or lateral displacement (or actually it's derivative) vs time. Ususally this is done by cross correlation.

I'll give your idea a shot (tomorrow morning local time), it might work. I'm not 100% sure how this procedure will handle drift, i.e., distortions, within the images. But I'll see and let you know. 

Hello HJ,

I confess that I do not understand your code.  First, you should never include DoUpdate in code that you are trying to optimize.  Even if it is not costing you much execution time, it introduces other variability (possibly updating other windows) that does not help you focus on the important code.

The second point has to do with your multithread calls.  The second call:

multithread result[idxP+range][idxQ+range]=sum(temp) 

evaluates a single wave point on the left hand side and sets it to the sum of the wave temp.  Since the right hand side is constant with respect to the left hand side, there is no point in multithread here.

MatrixOP supported (since 2014) the function subRange(w,rs,re,cs,ce) which allows you to extract (re-rs)x(ce-cs) pixels from anywhere in the image which you could then pass to correlate (in a single command line).

It is also worth mentioning that depending on your application, you may be able to use ImageRegistration to directly obtain the relative offset between two images.

AG

Hi AG, 

thank you for your input. The DoUpdate is really not harmful and shows that 'something' is going on (provided the result image is displayed in advance). It's one update per 2 secs or so (2000² px², range 50 px) and the final code would certainly not include such things. (The function would also return a wave reference to a free result wave and not use a wave in root, etc).

Optimization on for-loops over p/q indices is kind of futile anyway... 

If I understand your remark towards multithread correctly, Igor is parallelizing the left hand side of the equation (runing p1, p2, ... q1, q2,... parallel) and not running over 'segments' on the right hand side (like a sum of quartes (on e.g. 4 kernels) of subranges of a wave). Then MT is indeed useless here. I didn't time it, but it felt a bit faster (could be a personal bias thing).

Subrange is unfortunately not doing what I need. I need to use all pixels in the original images yet with limited shift to the original. Maybe the ppt-squares-image illustrates it a little bit. Referece in blue, shifted dataset in green (for inital and final offsets); left full correlation, right correlation around the center with a certain range (light blue). If range amounts to half the image size (in pnts/px; square) it should turn into a regular correlation. 

My code is padding the waves with zeros on all sides to avoid out of index errors and then calculates the central area of a correlation manually (at least that's the idea; the results were consistent in a few quick tests I did, see rainbow images; z values also match).

ImageRegistration might be worth a revisit (esp. sub-pixel resolution). I'll need to make sure that I convert my data to single presicion preserving all relevant information. Or is ImageRegistration also accepting DP waves, despite the manual states "All the input waves are expected to be single precision float"?

Cheers
HJ

Hello HJ,

ImageRegistration works with SP only.

If your calculation involves a kernel that is more than 13x13 then the FFT approach is computationally more efficient.

The images you attached suggest a high degree of symmetry which makes me wonder if you can't reduce this to a few 1D calculations.

AG

Hi all

SRC, the idea seems to work after offset subtraction (but that step is fine). Zapping [0,0] does not help, since its bleeding still dominates the spectrum.

AG, ImageRegistration works fast on the 'kind' test data I have here at the moment. I need to see how it behaves on a larger data set. Thank you for pushing me into the Image-funktions. One question: the flag /CONV=0 results in a "One or more parameters are inappropriate" error message. Does it need another flag in combination? (The flag /CON mentionend in first paragraph in details seems to be unknown). High symmetry is not guaranteed and can be an issue, actually (e.g., shifts of more than a unit cell). 

JJW, thank you for your idea. One image is around 10 MB, range will probably be useful around 50 (to cover the large inital displacements in the image sereies as well). Thus, the image stack will amount to ~25 GB. I see issues coming there. However, you just inspired me for a similar approach. Padding inside the loop and MatrixOP for the math might work as well.

Multithread for the sum is indeed not helpful: there is a jitter of ~1sec on the runtime (20 sec / 60sec / 250sec total as function of range). Sometimes the code with MT is faster, sometimes without (-> it's something else). 

Thank you for the discussion, I'll implement the FFT and ImageReg, see how they perfrom on a full data set, and come back with results after implementation.

Cheers, 
HJ

HJ, I'm glad that FFT was useful. I had neglected to mention DC offset adjustment. It does remove 0-frequency content, but some remaining low-frequency content may obscure signal discrimination. I suspect your signal "bleeding" from various irregularities will be a general problem, no matter what computational approach is used. You might try various types of FFT window filters, but the MatrixOP help file gives no indication that the "standard" FFT window-functions are supported.