diff options
Diffstat (limited to 'pkg/images/immatch/src/xregister/rgxfft.x')
| -rw-r--r-- | pkg/images/immatch/src/xregister/rgxfft.x | 179 |
1 files changed, 179 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/xregister/rgxfft.x b/pkg/images/immatch/src/xregister/rgxfft.x new file mode 100644 index 00000000..8847cf56 --- /dev/null +++ b/pkg/images/immatch/src/xregister/rgxfft.x @@ -0,0 +1,179 @@ +# RG_FFTCOR -- Compute the FFT of the reference and image data, take their +# product, and compute the inverse transform to get the cross-correlation +# function. The reference and input image are loaded into alternate memory +# locations. + +procedure rg_fftcor (fft, nxfft nyfft) + +real fft[ARB] #I/O array to be fft'd +int nxfft, nyfft #I dimensions of the fft + +pointer sp, dim + +begin + call smark (sp) + call salloc (dim, 2, TY_INT) + + # Fourier transform the two arrays. + Memi[dim] = nxfft + Memi[dim+1] = nyfft + if (Memi[dim+1] == 1) + call rg_fourn (fft, Memi[dim], 1, 1) + else + call rg_fourn (fft, Memi[dim], 2, 1) + + # Compute the product of the two transforms. + call rg_mulfft (fft, fft, 2 * nxfft, nyfft) + + # Shift the array to center the transform. + call rg_fshift (fft, fft, 2 * nxfft, nyfft) + + # Normalize the transform. + call adivkr (fft, real (nxfft * nyfft), fft, 2 * nxfft * nyfft) + + # Compute the inverse transform. + if (Memi[dim+1] == 1) + call rg_fourn (fft, Memi[dim], 1, -1) + else + call rg_fourn (fft, Memi[dim], 2, -1) + + call sfree (sp) +end + + +# RG_MULFFT -- Unpack the two individual ffts and compute their product. + +procedure rg_mulfft (fft1, fft2, nxfft, nyfft) + +real fft1[nxfft,nyfft] #I array containing 2 ffts of 2 real functions +real fft2[nxfft,nyfft] #O fft of correlation function +int nxfft, nyfft #I dimensions of fft + +int i,j, nxd2p2, nxp2, nxp3, nyd2p1, nyp2 +real c1, c2, h1r, h1i, h2r, h2i + +begin + c1 = 0.5 + c2 = -0.5 + + nxd2p2 = nxfft / 2 + 2 + nxp2 = nxfft + 2 + nxp3 = nxfft + 3 + nyd2p1 = nyfft / 2 + 1 + nyp2 = nyfft + 2 + + # Compute the 0 frequency point. + h1r = fft1[1,1] + h1i = 0.0 + h2r = fft1[2,1] + h2i = 0.0 + fft2[1,1] = h1r * h2r + fft2[2,1] = 0.0 + + # Compute the x axis points. + do i = 3, nxd2p2, 2 { + h2r = c1 * (fft1[i,1] + fft1[nxp2-i,1]) + h2i = c1 * (fft1[i+1,1] - fft1[nxp3-i,1]) + h1r = -c2 * (fft1[i+1,1] + fft1[nxp3-i,1]) + h1i = c2 * (fft1[i,1] - fft1[nxp2-i,1]) + fft2[i,1] = (h1r * h2r + h1i * h2i) + fft2[i+1,1] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,1] = fft2[i,1] + fft2[nxp3-i,1] = - fft2[i+1,1] + } + + # Quit if the transform is 1D. + if (nyfft < 2) + return + + # Compute the y axis points. + do i = 2, nyd2p1 { + h2r = c1 * (fft1[1,i] + fft1[1, nyp2-i]) + h2i = c1 * (fft1[2,i] - fft1[2,nyp2-i]) + h1r = -c2 * (fft1[2,i] + fft1[2,nyp2-i]) + h1i = c2 * (fft1[1,i] - fft1[1,nyp2-i]) + fft2[1,i] = (h1r * h2r + h1i * h2i) + fft2[2,i] = (h1i * h2r - h2i * h1r) + fft2[1,nyp2-i] = fft2[1,i] + fft2[2,nyp2-i] = - fft2[2,i] + } + + # Compute along the axis of symmetry. + do i = 3, nxd2p2, 2 { + h2r = c1 * (fft1[i,nyd2p1] + fft1[nxp2-i, nyd2p1]) + h2i = c1 * (fft1[i+1,nyd2p1] - fft1[nxp3-i,nyd2p1]) + h1r = -c2 * (fft1[i+1,nyd2p1] + fft1[nxp3-i,nyd2p1]) + h1i = c2 * (fft1[i,nyd2p1] - fft1[nxp2-i,nyd2p1]) + fft2[i,nyd2p1] = (h1r * h2r + h1i * h2i) + fft2[i+1,nyd2p1] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,nyd2p1] = fft2[i,nyd2p1] + fft2[nxp3-i,nyd2p1] = - fft2[i+1,nyd2p1] + } + + # Compute the remainder of the transform. + do j = 2, nyd2p1 - 1 { + do i = 3, nxfft, 2 { + h2r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j]) + h2i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j]) + h1r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j]) + h1i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j]) + fft2[i,j] = (h1r * h2r + h1i * h2i) + fft2[i+1,j] = (h1i * h2r - h2i * h1r) + fft2[nxp2-i,nyp2-j] = fft2[i,j] + fft2[nxp3-i,nyp2-j] = - fft2[i+1,j] + } + } +end + + +# RG_FNORM -- Normalize the reference and image data before computing +# the fft's. + +procedure rg_fnorm (array, ncols, nlines, nxfft, nyfft) + +real array[ARB] #I/O the input/output data array +int ncols, nlines #I dimensions of the input data array +int nxfft, nyfft #I dimensions of the fft + +int i, j, index +real sumr, sumi, meanr, meani + +begin + # Compute the mean. + sumr = 0.0 + sumi = 0.0 + index = 0 + do j = 1, nlines { + do i = 1, ncols { + sumr = sumr + array[index+2*i-1] + sumi = sumi + array[index+2*i] + } + index = index + 2 * nxfft + } + meanr = sumr / (ncols * nlines) + meani = sumi / (ncols * nlines) + + # Compute the sigma. + sumr = 0.0 + sumi = 0.0 + index = 0 + do j = 1, nlines { + do i = 1, ncols { + sumr = sumr + (array[index+2*i-1] - meanr) ** 2 + sumi = sumi + (array[index+2*i] - meani) ** 2 + } + index = index + 2 * nxfft + } + sumr = sqrt (sumr) + sumi = sqrt (sumi) + + # Normalize the data. + index = 0 + do j = 1, nlines { + do i = 1, ncols { + array[index+2*i-1] = (array[index+2*i-1] - meanr) / sumr + array[index+2*i] = (array[index+2*i] - meani) / sumi + } + index = index + 2 * nxfft + } +end |