diff options
Diffstat (limited to 'pkg/images/immatch/src/psfmatch/rgpfft.x')
| -rw-r--r-- | pkg/images/immatch/src/psfmatch/rgpfft.x | 443 |
1 files changed, 443 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/psfmatch/rgpfft.x b/pkg/images/immatch/src/psfmatch/rgpfft.x new file mode 100644 index 00000000..b5f36375 --- /dev/null +++ b/pkg/images/immatch/src/psfmatch/rgpfft.x @@ -0,0 +1,443 @@ + +# RG_PG10F -- Fetch the 0 component of the fft. + +real procedure rg_pg10f (fft, nxfft, nyfft) + +real fft[nxfft,nyfft] #I array containing 2 real ffts +int nxfft #I x dimension of complex array +int nyfft #I y dimension of complex array + +int xcen, ycen + +begin + xcen = nxfft / 2 + 1 + ycen = nyfft / 2 + 1 + + return (fft[xcen,ycen]) +end + + +# RG_PG1NORM -- Estimate the normalization factor by computing the amplitude +# of the best fitting Gaussian. This routine may eventually be replaced by +# on which does a complete Gaussian fit. The Gaussian is assumed to be +# of the form g = a * exp (b * r * r). The input array is a 2D real array +# storing 1 fft of dimension nxfft by nyfft in complex order with the +# zero frequency in the center. + +real procedure rg_pg1norm (fft, nxfft, nyfft) + +real fft[nxfft,nyfft] #I array containing 2 real ffts +int nxfft #I x dimension of complex array +int nyfft #I y dimension of complex array + +int xcen, ycen +real ln1, ln2, cx, cy + +begin + xcen = nxfft / 2 + 1 + ycen = nyfft / 2 + 1 + + if (nxfft >= 8) { + ln1 = log (sqrt (fft[xcen-2,ycen] ** 2 + fft[xcen-1,ycen] ** 2)) + ln2 = log (sqrt (fft[xcen-4,ycen] ** 2 + fft[xcen-3,ycen] ** 2)) + cx = exp ((4.0 * ln1 - ln2) / 3.0) + } else + cx = 0.0 + + if (nyfft >= 4) { + ln1 = log (sqrt (fft[xcen,ycen-1] ** 2 + fft[xcen+1,ycen-1] ** 2)) + ln2 = log (sqrt (fft[xcen,ycen-2] ** 2 + fft[xcen+1,ycen-2] ** 2)) + cy = exp ((4.0 * ln1 - ln2) / 3.0) + } else + cy = 0.0 + + if (cx <= 0.0) + return (cy) + else if (cy <= 0.0) + return (cx) + else + return (0.5 * (cx + cy)) +end + + +# RG_PG20F -- Fetch the 0 component of the fft. + +real procedure rg_pg20f (fft, nxfft, nyfft) + +real fft[nxfft,nyfft] #I array containing 2 real ffts +int nxfft #I x dimension of complex array +int nyfft #I y dimension of complex array + +int xcen, ycen + +begin + xcen = nxfft / 2 + 1 + ycen = nyfft / 2 + 1 + + return (fft[xcen,ycen] / fft[xcen+1,ycen]) +end + + +# RG_PG2NORM -- Estimate the normalization factor by computing the amplitude +# of the best fitting Gaussian. This routine may eventually be replaced by +# on which does a complete Gaussian fit. The Gaussian is assumed to be +# of the form g = a * exp (b * r * r). The input array is a 2D real array +# storing 2 2D ffts of dimension nxfft by nyfft in complex order with the +# zero frequency in the center. + +real procedure rg_pg2norm (fft, nxfft, nyfft) + +real fft[nxfft,nyfft] #I array containing 2 real ffts +int nxfft #I x dimension of complex array +int nyfft #I y dimension of complex array + +int xcen, ycen +real fftr, ffti, ln1r, ln2r, ln1i, ln2i, cxr, cyr, cxi, cyi, ampr, ampi + +begin + + xcen = nxfft / 2 + 1 + ycen = nyfft / 2 + 1 + + # Compute the x amplitude for the first fft. + if (nxfft >= 8) { + + fftr = 0.5 * (fft[xcen+2,ycen] + fft[xcen-2,ycen]) + ffti = 0.5 * (fft[xcen+3,ycen] - fft[xcen-1,ycen]) + ln1r = log (sqrt (fftr ** 2 + ffti ** 2)) + fftr = 0.5 * (fft[xcen+4,ycen] + fft[xcen-4,ycen]) + ffti = 0.5 * (fft[xcen+5,ycen] - fft[xcen-3,ycen]) + ln2r = log (sqrt (fftr ** 2 + ffti ** 2)) + + fftr = 0.5 * (fft[xcen+3,ycen] + fft[xcen-1,ycen]) + ffti = -0.5 * (fft[xcen+2,ycen] - fft[xcen-2,ycen]) + ln1i = log (sqrt (fftr ** 2 + ffti ** 2)) + fftr = 0.5 * (fft[xcen+5,ycen] + fft[xcen-3,ycen]) + ffti = -0.5 * (fft[xcen+4,ycen] - fft[xcen-4,ycen]) + ln2i = log (sqrt (fftr ** 2 + ffti ** 2)) + + cxr = exp ((4.0 * ln1r - ln2r) / 3.0) + cxi = exp ((4.0 * ln1i - ln2i) / 3.0) + + } else { + + cxr = 0.0 + cxi = 0.0 + + } + + # Compute the y ratio. + if (nyfft >= 4) { + + fftr = 0.5 * (fft[xcen,ycen+1] + fft[xcen,ycen-1]) + ffti = 0.5 * (fft[xcen+1,ycen+1] - fft[xcen+1,ycen-1]) + ln1r = log (sqrt (fftr ** 2 + ffti ** 2)) + fftr = 0.5 * (fft[xcen,ycen+2] + fft[xcen,ycen-2]) + ffti = 0.5 * (fft[xcen+1,ycen+2] - fft[xcen+1,ycen-2]) + ln2r = log (sqrt (fftr ** 2 + ffti ** 2)) + + fftr = 0.5 * (fft[xcen+1,ycen+1] + fft[xcen+1,ycen-1]) + ffti = -0.5 * (fft[xcen,ycen+1] - fft[xcen,ycen-1]) + ln1i = log (sqrt (fftr ** 2 + ffti ** 2)) + fftr = 0.5 * (fft[xcen+1,ycen+2] + fft[xcen+1,ycen-2]) + ffti = -0.5 * (fft[xcen,ycen+2] - fft[xcen,ycen-2]) + ln2i = log (sqrt (fftr ** 2 + ffti ** 2)) + + cyr = exp ((4.0 * ln1r - ln2r) / 3.0) + cyi = exp ((4.0 * ln1i - ln2i) / 3.0) + + } else { + + cyr = 0.0 + cyi = 0.0 + + } + + if (cxr <= 0.0) + ampr = cyr + else if (cyr <= 0.0) + ampr = cxr + else + ampr = 0.5 * (cxr + cyr) + + if (cxi <= 0.0) + ampi = cyi + else if (cyi <= 0.0) + ampi = cxi + else + ampi = 0.5 * (cxi + cyi) + + if (ampi <= 0.0) + return (INDEFR) + else + return (ampr /ampi) +end + + +# RG_PDIVFFT -- Unpack the two fft's, save the first fft, and compute the +# quotient of the two ffts. + +procedure rg_pdivfft (fft1, fftnum, fftdenom, fft2, nxfft, nyfft) + +real fft1[nxfft,nyfft] # array containing 2 ffts of 2 real functions +real fftnum[nxfft,nyfft] # the numerator fft +real fftdenom[nxfft,nyfft] # the denominator fft +real fft2[nxfft,nyfft] # fft of psf matching function +int nxfft, nyfft # dimensions of fft + +int i, j, xcen, ycen, nxp2, nxp3, nyp2 +real c1, c2, h1r, h1i, h2r, h2i, denom + +begin + c1 = 0.5 + c2 = -0.5 + xcen = nxfft / 2 + 1 + ycen = nyfft / 2 + 1 + nxp2 = nxfft + 2 + nxp3 = nxfft + 3 + nyp2 = nyfft + 2 + + # Compute the 0 frequency point. + h1r = fft1[xcen,ycen] + h1i = 0.0 + h2r = fft1[xcen+1,ycen] + h2i = 0.0 + fftnum[xcen,ycen] = h1r + fftnum[xcen+1,ycen] = 0.0 + fftdenom[xcen,ycen] = h2r + fftdenom[xcen+1,ycen] = 0.0 + fft2[xcen,ycen] = h1r / h2r + fft2[xcen+1,ycen] = 0.0 + + #call eprintf ("fft11=%g fft21=%g\n") + #call pargr (fft1[1,1]) + #call pargr (fft1[2,1]) + + # Compute the first point. + h1r = c1 * (fft1[1,1] + fft1[1,1]) + h1i = 0.0 + h2r = -c2 * (fft1[2,1] + fft1[2,1]) + h2i = 0.0 + + fftnum[1,1] = h1r + fftnum[2,1] = h1i + fftdenom[1,1] = h2r + fftdenom[2,1] = h2i + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[1,1] = 1.0 + fft2[2,1] = 0.0 + } else { + fft2[1,1] = (h1r * h2r + h1i * h2i) / denom + fft2[2,1] = (h1i * h2r - h2i * h1r) / denom + } + + # Compute the x symmetry axis points. + do i = 3, xcen - 1, 2 { + + h1r = c1 * (fft1[i,ycen] + fft1[nxp2-i,ycen]) + h1i = c1 * (fft1[i+1,ycen] - fft1[nxp3-i,ycen]) + h2r = -c2 * (fft1[i+1,ycen] + fft1[nxp3-i,ycen]) + h2i = c2 * (fft1[i,ycen] - fft1[nxp2-i,ycen]) + + fftnum[i,ycen] = h1r + fftnum[i+1,ycen] = h1i + fftnum[nxp2-i,ycen] = h1r + fftnum[nxp3-i,ycen] = -h1i + + fftdenom[i,ycen] = h2r + fftdenom[i+1,ycen] = h2i + fftdenom[nxp2-i,ycen] = h2r + fftdenom[nxp3-i,ycen] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[i,ycen] = 1.0 + fft2[i+1,ycen] = 0.0 + } else { + fft2[i,ycen] = (h1r * h2r + h1i * h2i) / denom + fft2[i+1,ycen] = (h1i * h2r - h2i * h1r) / denom + } + fft2[nxp2-i,ycen] = fft2[i,ycen] + fft2[nxp3-i,ycen] = -fft2[i+1,ycen] + + } + + # Quit if the transform is 1D. + if (nyfft < 2) + return + + # Compute the x axis points. + do i = 3, xcen + 1, 2 { + + h1r = c1 * (fft1[i,1] + fft1[nxp2-i,1]) + h1i = c1 * (fft1[i+1,1] - fft1[nxp3-i,1]) + h2r = -c2 * (fft1[i+1,1] + fft1[nxp3-i,1]) + h2i = c2 * (fft1[i,1] - fft1[nxp2-i,1]) + + fftnum[i,1] = h1r + fftnum[i+1,1] = h1i + fftnum[nxp2-i,1] = h1r + fftnum[nxp3-i,1] = -h1i + + fftdenom[i,1] = h2r + fftdenom[i+1,1] = h2i + fftdenom[nxp2-i,1] = h2r + fftdenom[nxp3-i,1] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0) { + fft2[i,1] = 1.0 + fft2[i+1,1] = 0.0 + } else { + fft2[i,1] = (h1r * h2r + h1i * h2i) / denom + fft2[i+1,1] = (h1i * h2r - h2i * h1r) / denom + } + fft2[nxp2-i,1] = fft2[i,1] + fft2[nxp3-i,1] = -fft2[i+1,1] + } + + # Compute the y symmetry axis points. + do i = 2, ycen - 1 { + + h1r = c1 * (fft1[xcen,i] + fft1[xcen, nyp2-i]) + h1i = c1 * (fft1[xcen+1,i] - fft1[xcen+1,nyp2-i]) + h2r = -c2 * (fft1[xcen+1,i] + fft1[xcen+1,nyp2-i]) + h2i = c2 * (fft1[xcen,i] - fft1[xcen,nyp2-i]) + + fftnum[xcen,i] = h1r + fftnum[xcen+1,i] = h1i + fftnum[xcen,nyp2-i] = h1r + fftnum[xcen+1,nyp2-i] = -h1i + + fftdenom[xcen,i] = h2r + fftdenom[xcen+1,i] = h2i + fftdenom[xcen,nyp2-i] = h2r + fftdenom[xcen+1,nyp2-i] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[xcen,i] = 1.0 + fft2[xcen+1,i] = 0.0 + } else { + fft2[xcen,i] = (h1r * h2r + h1i * h2i) / denom + fft2[xcen+1,i] = (h1i * h2r - h2i * h1r) / denom + } + fft2[xcen,nyp2-i] = fft2[xcen,i] + fft2[xcen+1,nyp2-i] = -fft2[xcen+1,i] + + } + + # Compute the y axis points. + do i = 2, ycen { + + h1r = c1 * (fft1[1,i] + fft1[1,nyp2-i]) + h1i = c1 * (fft1[2,i] - fft1[2,nyp2-i]) + h2r = -c2 * (fft1[2,i] + fft1[2,nyp2-i]) + h2i = c2 * (fft1[1,i] - fft1[1,nyp2-i]) + + fftnum[1,i] = h1r + fftnum[2,i] = h1i + fftnum[1,nyp2-i] = h1r + fftnum[2,nyp2-i] = -h1i + + fftdenom[1,i] = h2r + fftdenom[2,i] = h2i + fftdenom[1,nyp2-i] = h2r + fftdenom[2,nyp2-i] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[1,i] = 1.0 + fft2[2,i] = 0.0 + } else { + fft2[1,i] = (h1r * h2r + h1i * h2i) / denom + fft2[2,i] = (h1i * h2r - h2i * h1r) / denom + } + fft2[1,nyp2-i] = fft2[1,i] + fft2[2,nyp2-i] = -fft2[2,i] + } + + # Compute the remainder of the transform. + do j = 2, ycen - 1 { + + do i = 3, xcen - 1, 2 { + + h1r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j]) + h1i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j]) + h2r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j]) + h2i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j]) + + fftnum[i,j] = h1r + fftnum[i+1,j] = h1i + fftnum[nxp2-i,nyp2-j] = h1r + fftnum[nxp3-i,nyp2-j] = -h1i + + fftdenom[i,j] = h2r + fftdenom[i+1,j] = h2i + fftdenom[nxp2-i,nyp2-j] = h2r + fftdenom[nxp3-i,nyp2-j] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[i,j] = 1.0 + fft2[i+1,j] = 0.0 + } else { + fft2[i,j] = (h1r * h2r + h1i * h2i) / denom + fft2[i+1,j] = (h1i * h2r - h2i * h1r) / denom + } + fft2[nxp2-i,nyp2-j] = fft2[i,j] + fft2[nxp3-i,nyp2-j] = - fft2[i+1,j] + } + + do i = xcen + 2, nxfft, 2 { + + h1r = c1 * (fft1[i,j] + fft1[nxp2-i, nyp2-j]) + h1i = c1 * (fft1[i+1,j] - fft1[nxp3-i,nyp2-j]) + h2r = -c2 * (fft1[i+1,j] + fft1[nxp3-i,nyp2-j]) + h2i = c2 * (fft1[i,j] - fft1[nxp2-i,nyp2-j]) + + fftnum[i,j] = h1r + fftnum[i+1,j] = h1i + fftnum[nxp2-i,nyp2-j] = h1r + fftnum[nxp3-i,nyp2-j] = -h1i + + fftdenom[i,j] = h2r + fftdenom[i+1,j] = h2i + fftdenom[nxp2-i,nyp2-j] = h2r + fftdenom[nxp3-i,nyp2-j] = -h2i + + denom = h2r * h2r + h2i * h2i + if (denom == 0.0) { + fft2[i,j] = 1.0 + fft2[i+1,j] = 0.0 + } else { + fft2[i,j] = (h1r * h2r + h1i * h2i) / denom + fft2[i+1,j] = (h1i * h2r - h2i * h1r) / denom + } + fft2[nxp2-i,nyp2-j] = fft2[i,j] + fft2[nxp3-i,nyp2-j] = - fft2[i+1,j] + + } + } +end + + +# RG_PNORM -- Insert the normalization value into the 0 frequency of the +# fft. The fft is a 2D fft stored in a real array in complex order. +# The fft is assumed to be centered. + +procedure rg_pnorm (fft, nxfft, nyfft, norm) + +real fft[ARB] #I the input fft +int nxfft #I the x dimension of fft (complex storage) +int nyfft #I the y dimension of the fft +real norm #I the flux ratio + +int index + +begin + index = nxfft + 1 + 2 * (nyfft / 2) * nxfft + fft[index] = norm + fft[index+1] = 0.0 +end |