aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/immatch/src/xregister/rgxfit.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/immatch/src/xregister/rgxfit.x')
-rw-r--r--pkg/images/immatch/src/xregister/rgxfit.x814
1 files changed, 814 insertions, 0 deletions
diff --git a/pkg/images/immatch/src/xregister/rgxfit.x b/pkg/images/immatch/src/xregister/rgxfit.x
new file mode 100644
index 00000000..34e6398c
--- /dev/null
+++ b/pkg/images/immatch/src/xregister/rgxfit.x
@@ -0,0 +1,814 @@
+include <mach.h>
+include <math/iminterp.h>
+include <math/nlfit.h>
+include "xregister.h"
+
+define NL_MAXITER 10
+define NL_TOL 0.001
+
+# RG_FIT -- Fit the peak of the cross-correlation function using one of the
+# fitting functions.
+
+procedure rg_fit (xc, nreg, gd, xshift, yshift)
+
+pointer xc #I the pointer to the cross-corrrelation structure
+int nreg #I the current region
+pointer gd #I the pointer to the graphics stream
+real xshift, yshift #O the computed shifts
+
+int nrlines, xwindow, ywindow, xcbox, ycbox, xlag, ylag
+real xin, yin, xout, yout
+int rg_xstati()
+pointer rg_xstatp()
+
+begin
+ # Check the window and centering box sizes.
+ nrlines = Memi[rg_xstatp(xc,RL2)+nreg-1] -
+ Memi[rg_xstatp(xc,RL1)+nreg-1] + 1
+ xwindow = rg_xstati (xc, XWINDOW)
+ if (nrlines == 1)
+ ywindow = 1
+ else
+ ywindow = rg_xstati (xc, YWINDOW)
+ xcbox = rg_xstati (xc, XCBOX)
+ if (nrlines == 1)
+ ycbox = 1
+ else
+ ycbox = rg_xstati (xc, YCBOX)
+
+ # Do the centering.
+ switch (rg_xstati (xc, PFUNC)) {
+ case XC_PNONE:
+ call rg_maxmin (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow,
+ xshift, yshift)
+ case XC_CENTROID:
+ call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow,
+ ywindow, xcbox, ycbox, xshift, yshift)
+ case XC_SAWTOOTH:
+ call rg_sawtooth (Memr[rg_xstatp(xc,XCOR)], xwindow,
+ ywindow, xcbox, ycbox, xshift, yshift)
+ case XC_PARABOLA:
+ call rg_iparabolic (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow,
+ xcbox, ycbox, xshift, yshift)
+ case XC_MARK:
+ if (gd == NULL)
+ call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow,
+ ywindow, xcbox, ycbox, xshift, yshift)
+ else
+ call rg_xmkpeak (gd, xwindow, ywindow, xshift, yshift)
+ default:
+ call rg_imean (Memr[rg_xstatp(xc,XCOR)], xwindow, ywindow,
+ xcbox, ycbox, xshift, yshift)
+ }
+
+ # Store the shifts.
+ if (rg_xstati (xc, NREFPTS) > 0) {
+ xin = (Memi[rg_xstatp(xc,RC1)+nreg-1] +
+ Memi[rg_xstatp(xc,RC2)+nreg-1]) / 2.0
+ yin = (Memi[rg_xstatp(xc,RL1)+nreg-1] +
+ Memi[rg_xstatp(xc,RL2)+nreg-1]) / 2.0
+ call rg_etransform (xc, xin, yin, xout, yout)
+ xlag = xout - xin
+ ylag = yout - yin
+ } else {
+ xlag = rg_xstati (xc, XLAG)
+ ylag = rg_xstati (xc, YLAG)
+ }
+ xshift = - (xshift + xlag)
+ yshift = - (yshift + ylag)
+ Memr[rg_xstatp(xc,XSHIFTS)+nreg-1] = xshift
+ Memr[rg_xstatp(xc,YSHIFTS)+nreg-1] = yshift
+end
+
+
+# RG_MAXMIN -- Procedure to compute the peak of the cross-correlation function
+# by determining the maximum point.
+
+procedure rg_maxmin (xcor, xwindow, ywindow, xshift, yshift)
+
+real xcor[xwindow,ywindow] #I the cross-correlation function
+int xwindow, ywindow #I dimensions of cross-correlation function
+real xshift, yshift #O x and shift of the peak
+
+int xindex, yindex
+
+begin
+ # Locate the maximum point.
+ call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+ xshift = xindex - (1.0 + xwindow) / 2.0
+ yshift = yindex - (1.0 + ywindow) / 2.0
+end
+
+
+# RG_IMEAN -- Compute the peak of the cross-correlation function using the
+# intensity weighted mean of the marginal distributions in x and y.
+
+procedure rg_imean (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real xcor[xwindow,ARB] #I the cross-correlation function
+int xwindow, ywindow #I dimensions of the cross-correlation function
+int xcbox, ycbox #I dimensions of the centering box
+real xshift, yshift #O x and y shift of cross-correlation function
+
+int xindex, yindex, xlo, xhi, ylo, yhi, nx, ny
+pointer sp, xmarg, ymarg
+
+begin
+ call smark (sp)
+ call salloc (xmarg, xcbox, TY_REAL)
+ call salloc (ymarg, ycbox, TY_REAL)
+
+ # Locate the maximum point and normalize.
+ call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+ # Compute the limits of the centering box.
+ xlo = max (1, xindex - xcbox / 2)
+ xhi = min (xwindow, xindex + xcbox / 2)
+ nx = xhi - xlo + 1
+ ylo = max (1, yindex - ycbox / 2)
+ yhi = min (ywindow, yindex + ycbox / 2)
+ ny = yhi - ylo + 1
+
+ # Accumulate the marginals.
+ call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi,
+ Memr[xmarg], Memr[ymarg])
+
+ # Compute the shifts.
+ call rg_centroid (Memr[xmarg], nx, xshift)
+ xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+ call rg_centroid (Memr[ymarg], ny, yshift)
+ yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+ call sfree (sp)
+end
+
+
+# RG_IPARABOLIC -- Computer the peak of the cross-correlation function by
+# doing parabolic interpolation around the peak.
+
+procedure rg_iparabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real xcor[xwindow,ARB] #I the cross-correlation function
+int xwindow, ywindow #I dimensions of the cross-correlation fucntion
+int xcbox, ycbox #I the dimensions of the centering box
+real xshift, yshift #O the x and y shift of the peak
+
+int i, j, xindex, yindex, xlo, xhi, nx, ylo, yhi, ny
+pointer sp, x, y, c, xfit, yfit
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (x, 3, TY_REAL)
+ call salloc (y, 3, TY_REAL)
+ call salloc (c, 3, TY_REAL)
+ call salloc (xfit, 3, TY_REAL)
+ call salloc (yfit, 3, TY_REAL)
+
+ # Locate the maximum point.
+ call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+ xlo = max (1, xindex - 1)
+ xhi = min (xwindow, xindex + 1)
+ nx = xhi - xlo + 1
+ ylo = max (1, yindex - 1)
+ yhi = min (ywindow, yindex + 1)
+ ny = yhi - ylo + 1
+
+ # Initialize.
+ do i = 1, 3
+ Memr[x+i-1] = i
+
+ # Fit the x shift.
+ if (nx >= 3) {
+ do j = ylo, yhi {
+ do i = xlo, xhi
+ Memr[y+i-xlo] = xcor[i,j]
+ call rg_iparab (Memr[x], Memr[y], Memr[c])
+ Memr[xfit+j-ylo] = - Memr[c+1] / (2.0 * Memr[c+2])
+ Memr[yfit+j-ylo] = Memr[c] + Memr[c+1] * Memr[xfit+j-ylo] +
+ Memr[c+2] * Memr[xfit+j-ylo] ** 2
+ }
+ if (ny >= 3)
+ call rg_iparab (Memr[xfit], Memr[yfit], Memr[c])
+ xshift = - Memr[c+1] / (2.0 * Memr[c+2])
+ } else
+ xshift = xindex - xlo + 1
+
+ # Fit the y shift.
+ if (ny >= 3) {
+ do i = xlo, xhi {
+ do j = ylo, yhi
+ Memr[y+j-ylo] = xcor[i,j]
+ call rg_iparab (Memr[x], Memr[y], Memr[c])
+ Memr[xfit+i-xlo] = - Memr[c+1] / (2.0 * Memr[c+2])
+ Memr[yfit+i-xlo] = Memr[c] + Memr[c+1] * Memr[xfit+i-xlo] +
+ Memr[c+2] * Memr[xfit+i-xlo] ** 2
+ }
+ call rg_iparab (Memr[xfit], Memr[yfit], Memr[c])
+ yshift = - Memr[c+1] / (2.0 * Memr[c+2])
+ } else
+ yshift = yindex - ylo + 1
+
+ xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+ yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+ call sfree (sp)
+end
+
+
+define NPARS_PARABOLA 3
+
+# RG_PARABOLIC -- Compute the peak of the cross-correlation function by fitting
+# a parabola to the peak.
+
+procedure rg_parabolic (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real xcor[xwindow,ARB] #I the cross-correlation function
+int xwindow, ywindow #I dimensions of the cross-correlation fucntion
+int xcbox, ycbox #I the dimensions of the centering box
+real xshift, yshift #O the x and y shift of the peak
+
+extern rg_polyfit, rg_dpolyfit
+int i, xindex, yindex, xlo, xhi, ylo, yhi, nx, ny, npar, ier
+pointer sp, x, w, xmarg, ymarg, params, eparams, list, nl
+int locpr()
+
+begin
+ call smark (sp)
+ call salloc (x, max (xwindow, ywindow), TY_REAL)
+ call salloc (w, max (xwindow, ywindow), TY_REAL)
+ call salloc (xmarg, max (xwindow, ywindow), TY_REAL)
+ call salloc (ymarg, max (xwindow, ywindow), TY_REAL)
+ call salloc (params, NPARS_PARABOLA, TY_REAL)
+ call salloc (eparams, NPARS_PARABOLA, TY_REAL)
+ call salloc (list, NPARS_PARABOLA, TY_INT)
+
+ # Locate the maximum point.
+ call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+ xlo = max (1, xindex - xcbox / 2)
+ xhi = min (xwindow, xindex + xcbox / 2)
+ nx = xhi - xlo + 1
+ ylo = max (1, yindex - ycbox / 2)
+ yhi = min (ywindow, yindex + ycbox / 2)
+ ny = yhi - ylo + 1
+
+ # Accumulate the marginals.
+ call rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi,
+ Memr[xmarg], Memr[ymarg])
+
+ # Compute the x shift.
+ if (nx >= 3) {
+ do i = 1, nx
+ Memr[x+i-1] = i
+ do i = 1, nx
+ Memr[w+i-1] = Memr[xmarg+i-1]
+ call rg_iparab (Memr[x+xindex-xlo-1], Memr[xmarg+xindex-xlo-1],
+ Memr[params])
+ xshift = - Memr[params+1] / (2.0 * Memr[params+2])
+ call eprintf ("\txshift=%g\n")
+ call pargr (xshift)
+ call aclrr (Memr[eparams], NPARS_PARABOLA)
+ do i = 1, NPARS_PARABOLA
+ Memi[list+i-1] = i
+ call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit),
+ Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list],
+ NPARS_PARABOLA, .0001, NL_MAXITER)
+ call nlfitr (nl, Memr[x], Memr[xmarg], Memr[w], nx, 1, WTS_USER,
+ ier)
+ call nlvectorr (nl, Memr[x], Memr[w], nx, 1)
+ do i = 1, nx {
+ call eprintf ("x=%g y=%g yfit=%g\n")
+ call pargr (Memr[x+i-1])
+ call pargr (Memr[xmarg+i-1])
+ call pargr (Memr[w+i-1])
+ }
+ if (ier != NO_DEG_FREEDOM) {
+ call nlpgetr (nl, Memr[params], npar)
+ if (Memr[params+2] != 0)
+ xshift = - Memr[params+1] / (2.0 * Memr[params+2])
+ else
+ xshift = xindex - xlo + 1
+ } else
+ xshift = xindex - xlo + 1
+ call nlfreer (nl)
+ } else
+ xshift = xindex - xlo + 1
+
+ # Compute the y shift.
+ if (ny >= 3) {
+ do i = 1, ny
+ Memr[x+i-1] = i
+ do i = 1, ny
+ Memr[w+i-1] = Memr[ymarg+i-1]
+ call rg_iparab (Memr[x+yindex-ylo-1], Memr[ymarg+yindex-ylo-1],
+ Memr[params])
+ yshift = - Memr[params+1] / (2.0 * Memr[params+2])
+ call eprintf ("\tyshift=%g\n")
+ call pargr (yshift)
+ call aclrr (Memr[eparams], NPARS_PARABOLA)
+ do i = 1, NPARS_PARABOLA
+ Memi[list+i-1] = i
+ call nlinitr (nl, locpr (rg_polyfit), locpr (rg_dpolyfit),
+ Memr[params], Memr[eparams], NPARS_PARABOLA, Memi[list],
+ NPARS_PARABOLA, 0.0001, NL_MAXITER)
+ call nlfitr (nl, Memr[x], Memr[ymarg], Memr[w], ny, 1, WTS_USER,
+ ier)
+ call nlvectorr (nl, Memr[x], Memr[w], ny, 1)
+ do i = 1, ny {
+ call eprintf ("x=%g y=%g yfit=%g\n")
+ call pargr (Memr[x+i-1])
+ call pargr (Memr[ymarg+i-1])
+ call pargr (Memr[w+i-1])
+ }
+ if (ier != NO_DEG_FREEDOM) {
+ call nlpgetr (nl, Memr[params], npar)
+ if (Memr[params+2] != 0)
+ yshift = -Memr[params+1] / (2.0 * Memr[params+2])
+ else
+ yshift = yindex - ylo + 1
+ } else
+ yshift = yindex - ylo + 1
+ call nlfreer (nl)
+ } else
+ yshift = yindex - ylo + 1
+
+ xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+ yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+ call sfree (sp)
+end
+
+define EMISSION 1 # emission features
+define ABSORPTION 2 # emission features
+
+# RG_SAWTOOTH -- Compute the the x and y centers using a sawtooth
+# convolution function.
+
+procedure rg_sawtooth (xcor, xwindow, ywindow, xcbox, ycbox, xshift, yshift)
+
+real xcor[xwindow,ARB] #I the cross-correlation function
+int xwindow, ywindow #I the dimensions of the cross-correlation
+int xcbox, ycbox #I the dimensions of the centering box
+real xshift, yshift #O the x and y shifts
+
+int i, j, xindex, yindex, xlo, xhi, ylo, yhi, nx, ny
+pointer sp, data, xfit, yfit, yclean
+real ic
+
+begin
+ call smark (sp)
+ call salloc (data, max (xwindow, ywindow), TY_REAL)
+ call salloc (xfit, max (xwindow, ywindow), TY_REAL)
+ call salloc (yfit, max (xwindow, ywindow), TY_REAL)
+ call salloc (yclean, max (xwindow, ywindow), TY_REAL)
+
+ # Locate the maximum point and normalize.
+ call rg_alim2r (xcor, xwindow, ywindow, xindex, yindex)
+
+ xlo = max (1, xindex - xcbox)
+ xhi = min (xwindow, xindex + xcbox)
+ nx = xhi - xlo + 1
+ ylo = max (1, yindex - ycbox)
+ yhi = min (ywindow, yindex + ycbox)
+ ny = yhi - ylo + 1
+
+ # Compute the y shift.
+ if (ny >= 3) {
+ do j = ylo, yhi {
+ do i = xlo, xhi
+ Memr[data+i-xlo] = xcor[i,j]
+ call rg_x1dcenter (real (xindex - xlo + 1), Memr[data], nx,
+ Memr[xfit+j-ylo], Memr[yfit+j-ylo], real (nx / 2.0),
+ EMISSION, real (nx / 2.0), 0.0)
+ }
+ call arbpix (Memr[yfit], Memr[yclean], ny, II_SPLINE3,
+ II_BOUNDARYEXT)
+ call rg_x1dcenter (real (yindex - ylo + 1), Memr[yclean], ny,
+ yshift, ic, real (ny / 2.0), EMISSION, real (ny / 2.0), 0.0)
+ if (IS_INDEFR(yshift))
+ yshift = yindex - ylo + 1
+ } else
+ yshift = yindex - ylo + 1
+ yshift = yshift + ylo - 1 - (1.0 + ywindow) / 2.0
+
+ # Compute the x shift.
+ if (nx >= 3) {
+ if (ny >= 3) {
+ do i = xlo, xhi {
+ do j = ylo, yhi
+ Memr[data+j-ylo] = xcor[i,j]
+ call rg_x1dcenter (real (yindex - ylo + 1), Memr[data], ny,
+ Memr[xfit+i-xlo], Memr[yfit+i-xlo], real (ny / 2.0),
+ EMISSION, real (ny / 2.0), 0.0)
+ }
+ call arbpix (Memr[yfit], Memr[yclean], nx, II_SPLINE3,
+ II_BOUNDARYEXT)
+ call rg_x1dcenter (real (xindex - xlo + 1), Memr[yclean], nx,
+ xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0)
+ } else {
+ call rg_x1dcenter (real (xindex - xlo + 1), xcor[xlo,1], nx,
+ xshift, ic, real (nx / 2.0), EMISSION, real (nx / 2.0), 0.0)
+ }
+ if (IS_INDEFR(xshift))
+ xshift = xindex - xlo + 1
+ } else
+ xshift = xindex - xlo + 1
+ xshift = xshift + xlo - 1 - (1.0 + xwindow) / 2.0
+
+ call sfree (sp)
+end
+
+
+# RG_ALIM2R -- Determine the pixel position of the data maximum.
+
+procedure rg_alim2r (data, nx, ny, i, j)
+
+real data[nx,ARB] #I the input data
+int nx, ny #I the dimensions of the input array
+int i, j #O the indices of the maximum pixel
+
+int ii, jj
+real datamax
+
+begin
+ datamax = -MAX_REAL
+ do jj = 1, ny {
+ do ii = 1, nx {
+ if (data[ii,jj] > datamax) {
+ datamax = data[ii,jj]
+ i = ii
+ j = jj
+ }
+ }
+ }
+end
+
+
+# RG_XMKMARG -- Acumulate the marginal arrays in x and y.
+
+procedure rg_xmkmarg (xcor, xwindow, ywindow, xlo, xhi, ylo, yhi, xmarg,
+ ymarg)
+
+real xcor[xwindow,ARB] #I the cross-correlation function
+int xwindow, ywindow #I dimensions of cross-correlation function
+int xlo, xhi #I the x limits for centering
+int ylo, yhi #I the y limits for centering
+real xmarg[ARB] #O the output x marginal array
+real ymarg[ARB] #O the output y marginal array
+
+int i, j, index, nx, ny
+
+begin
+ nx = xhi - xlo + 1
+ ny = yhi - ylo + 1
+
+ # Compute the x marginal.
+ index = 1 - xlo
+ do i = xlo, xhi {
+ xmarg[index+i] = 0.0
+ do j = ylo, yhi
+ xmarg[index+i] = xmarg[index+i] + xcor[i,j]
+ }
+
+ # Normalize the x marginal.
+ call adivkr (xmarg, real (ny), xmarg, nx)
+
+ # Compute the y marginal.
+ index = 1 - ylo
+ do j = ylo, yhi {
+ ymarg[index+j] = 0.0
+ do i = xlo, xhi
+ ymarg[index+j] = ymarg[index+j] + xcor[i,j]
+ }
+
+ # Normalize the ymarginal.
+ call adivkr (ymarg, real (nx), ymarg, ny)
+end
+
+
+# RG_CENTROID -- Compute the intensity weighted maximum of an array.
+
+procedure rg_centroid (a, npts, shift)
+
+real a[ARB] #I the input array
+int npts #I the number of points
+real shift #O the position of the maximum
+
+int i
+real mean, dif, sumi, sumix
+bool fp_equalr()
+real asumr()
+
+begin
+ sumi = 0.0
+ sumix = 0.0
+ mean = asumr (a, npts) / npts
+
+ do i = 1, npts {
+ dif = a[i]
+ dif = a[i] - mean
+ if (dif < 0.0)
+ next
+ sumi = sumi + dif
+ sumix = sumix + i * dif
+ }
+
+ if (fp_equalr (sumi, 0.0))
+ shift = (1.0 + npts) / 2.0
+ else
+ shift = sumix / sumi
+end
+
+
+define MIN_WIDTH 3. # minimum centering width
+define EPSILON 0.001 # accuracy of centering
+define EPSILON1 0.005 # tolerance for convergence check
+define ITERATIONS 100 # maximum number of iterations
+define MAX_DXCHECK 3 # look back for failed convergence
+define INTERPTYPE II_SPLINE3 # image interpolation type
+
+
+# RG_X1DCENTER -- Locate the center of a one dimensional feature.
+# A value of INDEF is returned in the centering fails for any reason.
+# This procedure just sets up the data and adjusts for emission or
+# absorption features. The actual centering is done by C1D_CENTER.
+
+procedure rg_x1dcenter (x, data, npts, xc, ic, width, type, radius, threshold)
+
+real x #I initial guess
+real data[npts] #I data points
+int npts #I number of data points
+real xc #O computed center
+real ic #O intensity at computed center
+real width #I feature width
+int type #I feature type
+real radius #I centering radius
+real threshold #I minimum range in feature
+
+int x1, x2, nx
+real a, b, rad, wid
+pointer sp, data1
+
+begin
+ # Check starting value.
+ if (IS_INDEF(x) || (x < 1) || (x > npts)) {
+ xc = INDEF
+ ic = INDEF
+ return
+ }
+
+ # Set minimum width and error radius. The minimum in the error radius
+ # is for defining the data window. The user error radius is used to
+ # check for an error in the derived center at the end of the centering.
+
+ wid = max (width, MIN_WIDTH)
+ rad = max (2., radius)
+
+ # Determine the pixel value range around the initial center, including
+ # the width and error radius buffer. Check for a minimum range.
+
+ x1 = max (1., x - wid / 2 - rad - wid)
+ x2 = min (real (npts), x + wid / 2 + rad + wid + 1)
+ nx = x2 - x1 + 1
+ call alimr (data[x1], nx, a, b)
+ if (b - a < threshold) {
+ xc = INDEF
+ ic = INDEF
+ return
+ }
+
+ # Allocate memory for the continuum subtracted data vector. The X
+ # range is just large enough to include the error radius and the
+ # half width.
+
+ x1 = max (1., x - wid / 2 - rad)
+ x2 = min (real (npts), x + wid / 2 + rad + 1)
+ nx = x2 - x1 + 1
+
+ call smark (sp)
+ call salloc (data1, nx, TY_REAL)
+ call amovr (data[x1], Memr[data1], nx)
+
+ # Make the centering data positive, subtract the continuum, and
+ # apply a threshold to eliminate noise spikes.
+
+ switch (type) {
+ case EMISSION:
+ a = 0.
+ call asubkr (data[x1], a + threshold, Memr[data1], nx)
+ call amaxkr (Memr[data1], 0., Memr[data1], nx)
+ case ABSORPTION:
+ call anegr (data[x1], Memr[data1], nx)
+ call asubkr (Memr[data1], threshold - b, Memr[data1], nx)
+ call amaxkr (Memr[data1], 0., Memr[data1], nx)
+ default:
+ call error (0, "Unknown feature type")
+ }
+
+ # Determine the center.
+ call rg_xcenter (x - x1 + 1, Memr[data1], nx, xc, ic, wid)
+
+ # Check user centering error radius.
+ if (!IS_INDEF(xc)) {
+ xc = xc + x1 - 1
+ if (abs (x - xc) > radius) {
+ xc = INDEF
+ ic = INDEF
+ }
+ }
+
+ # Free memory and return the center position.
+ call sfree (sp)
+end
+
+
+# RG_XCENTER -- One dimensional centering algorithm.
+
+procedure rg_xcenter (x, data, npts, xc, ic, width)
+
+real x #I starting guess
+int npts #I number of points in data vector
+real data[npts] #I data vector
+real xc #O computed xc
+real ic #O computed intensity at xc
+real width #I centering width
+
+int i, j, iteration, dxcheck
+real hwidth, dx, dxabs, dxlast
+real a, b, sum1, sum2, intgrl1, intgrl2
+pointer asi1, asi2, sp, data1
+
+real asigrl(), asieval()
+
+define done_ 99
+
+begin
+ # Find the nearest local maxima as the starting point.
+ # This is required because the threshold limit may have set
+ # large regions of the data to zero and without a gradient
+ # the centering will fail.
+
+ i = x
+ for (i=x+.5; (i<npts) && (data[i]<=data[i+1]); i=i+1)
+ ;
+ for (j=x+.5; (j>1) && (data[j]<=data[j-1]); j=j-1)
+ ;
+
+ if (i-x < x-j)
+ xc = i
+ else
+ xc = j
+
+ # Check data range.
+ hwidth = width / 2
+ if ((xc - hwidth < 1) || (xc + hwidth > npts)) {
+ xc = INDEF
+ ic = INDEF
+ return
+ }
+
+ # Set interpolation functions.
+ call asiinit (asi1, INTERPTYPE)
+ call asiinit (asi2, INTERPTYPE)
+ call asifit (asi1, data, npts)
+
+ # Allocate, compute, and interpolate the x*y values.
+ call smark (sp)
+ call salloc (data1, npts, TY_REAL)
+ do i = 1, npts
+ Memr[data1+i-1] = data[i] * i
+ call asifit (asi2, Memr[data1], npts)
+ call sfree (sp)
+
+ # Iterate to find center. This loop exits when 1) the maximum
+ # number of iterations is reached, 2) the delta is less than
+ # the required accuracy (criterion for finding a center), 3)
+ # there is a problem in the computation, 4) successive steps
+ # continue to exceed the minimum delta.
+
+ dxlast = 1.
+ do iteration = 1, ITERATIONS {
+
+ # Triangle centering function.
+ a = xc - hwidth
+ b = xc - hwidth / 2
+ intgrl1 = asigrl (asi1, a, b)
+ intgrl2 = asigrl (asi2, a, b)
+ sum1 = (xc - hwidth) * intgrl1 - intgrl2
+ sum2 = -intgrl1
+ a = b
+ b = xc + hwidth / 2
+ intgrl1 = asigrl (asi1, a, b)
+ intgrl2 = asigrl (asi2, a, b)
+ sum1 = sum1 - xc * intgrl1 + intgrl2
+ sum2 = sum2 + intgrl1
+ a = b
+ b = xc + hwidth
+ intgrl1 = asigrl (asi1, a, b)
+ intgrl2 = asigrl (asi2, a, b)
+ sum1 = sum1 + (xc + hwidth) * intgrl1 - intgrl2
+ sum2 = sum2 - intgrl1
+
+ # Return no center if sum2 is zero.
+ if (sum2 == 0.)
+ break
+
+ # Limit dx change in one iteration to 1 pixel.
+ dx = max (-1., min (1., sum1 / abs (sum2)))
+ dxabs = abs (dx)
+ xc = xc + dx
+ ic = asieval (asi1, xc)
+
+ # Check data range. Return no center if at edge of data.
+ if ((xc - hwidth < 1) || (xc + hwidth > npts))
+ break
+
+ # Convergence tests.
+ if (dxabs < EPSILON)
+ goto done_
+ if (dxabs > dxlast + EPSILON1) {
+ dxcheck = dxcheck + 1
+ if (dxcheck > MAX_DXCHECK)
+ break
+ } else {
+ dxcheck = 0
+ dxlast = dxabs
+ }
+ }
+
+ # If we get here then no center was found.
+ xc = INDEF
+ ic = INDEF
+
+done_ call asifree (asi1)
+ call asifree (asi2)
+end
+
+
+# RG_IPARAB -- Compute the coefficients of the parabola through three
+# evenly spaced points.
+
+procedure rg_iparab (x, y, c)
+
+real x[NPARS_PARABOLA] #I input x values
+real y[NPARS_PARABOLA] #I input y values
+real c[NPARS_PARABOLA] #O computed coefficients
+
+begin
+ c[3] = (y[1]-y[2]) * (x[2]-x[3]) / (x[1]-x[2]) - (y[2]-y[3])
+ c[3] = c[3] / ((x[1]**2-x[2]**2) * (x[2]-x[3]) / (x[1]-x[2]) -
+ (x[2]**2-x[3]**2))
+
+ c[2] = (y[1] - y[2]) - c[3] * (x[1]**2 - x[2]**2)
+ c[2] = c[2] / (x[1] - x[2])
+
+ c[1] = y[1] - c[2] * x[1] - c[3] * x[1]**2
+end
+
+
+# RG_POLYFIT -- Evaluate an nth order polynomial.
+
+procedure rg_polyfit (x, nvars, p, np, z)
+
+real x #I position coordinate
+int nvars #I number of variables
+real p[ARB] #I coefficients of polynomial
+int np #I number of parameters
+real z #O function return
+
+int i
+real r
+
+begin
+ r = 0.0
+ do i = 2, np
+ r = r + x**(i-1) * p[i]
+ z = p[1] + r
+end
+
+
+# RG_DPOLYFIT -- Evaluate an nth order polynomial and its derivatives.
+
+procedure rg_dpolyfit (x, nvars, p, dp, np, z, der)
+
+real x #I position coordinate
+int nvars #I number of variables
+real p[ARB] #I coefficients of polynomial
+real dp[ARB] #I parameter derivative increments
+int np #I number of parameters
+real z #O function value
+real der[ARB] #O derivatives
+
+int i
+
+begin
+ der[1] = 1.0
+ z = 0.0
+ do i = 2, np {
+ der[i] = x ** (i-1)
+ z = z + x**(i-1) * p[i]
+ }
+ z = p[1] + z
+end
generated by cgit (git 2.34.1) at 2025-09-23 08:21:41 -0400