aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/lib/rgbckgrd.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/lib/rgbckgrd.x')
-rw-r--r--pkg/images/lib/rgbckgrd.x661
1 files changed, 661 insertions, 0 deletions
diff --git a/pkg/images/lib/rgbckgrd.x b/pkg/images/lib/rgbckgrd.x
new file mode 100644
index 00000000..feef4fd4
--- /dev/null
+++ b/pkg/images/lib/rgbckgrd.x
@@ -0,0 +1,661 @@
+include <mach.h>
+include <math.h>
+include <math/gsurfit.h>
+
+
+# RG_BORDER -- Fetch the border pixels from a 2D subraster.
+
+int procedure rg_border (buf, nx, ny, pnx, pny, ptr)
+
+real buf[nx,ARB] #I the input data subraster
+int nx, ny #I the dimensions of the input subraster
+int pnx, pny #I the size of the data region
+pointer ptr #I the pointer to the output buffer
+
+int j, nborder, wxborder, wyborder, index
+
+begin
+ # Compute the size of the array
+ nborder = nx * ny - pnx * pny
+ if (nborder <= 0) {
+ ptr = NULL
+ return (0)
+ } else if (nborder >= nx * ny) {
+ call malloc (ptr, nx * ny, TY_REAL)
+ call amovr (buf, Memr[ptr], nx * ny)
+ return (nx * ny)
+ } else
+ call malloc (ptr, nborder, TY_REAL)
+
+ # Fill the array.
+ wxborder = (nx - pnx) / 2
+ wyborder = (ny - pny) / 2
+ index = ptr
+ do j = 1, wyborder {
+ call amovr (buf[1,j], Memr[index], nx)
+ index = index + nx
+ }
+ do j = wyborder + 1, ny - wyborder {
+ call amovr (buf[1,j], Memr[index], wxborder)
+ index = index + wxborder
+ call amovr (buf[nx-wxborder+1,j], Memr[index], wxborder)
+ index = index + wxborder
+ }
+ do j = ny - wyborder + 1, ny {
+ call amovr (buf[1,j], Memr[index], nx)
+ index = index + nx
+ }
+
+ return (nborder)
+end
+
+
+# RG_SUBTRACT -- Subtract a plane from the data.
+
+procedure rg_subtract (data, nx, ny, zero, xslope, yslope)
+
+real data[nx,ARB] #I/O the input/output data array
+int nx, ny #I the dimensions of the input data array
+real zero #I the input zero point
+real xslope #I the input x slope
+real yslope #I the input y slope
+
+int i, j
+real ydelta
+
+begin
+ do j = 1, ny {
+ ydelta = yslope * j
+ do i = 1, nx
+ data[i,j] = data[i,j] - zero - xslope * i - ydelta
+ }
+end
+
+
+# RG_APODIZE -- Apply a cosine bell to the data. The operation can be
+# performed in place
+
+procedure rg_apodize (data, nx, ny, apodize, forward)
+
+real data[nx,ARB] #I the input data array
+int nx, ny #I the size of the input array
+real apodize #I the percentage of the end to apodize
+int forward #I YES for forward, NO for reverse
+
+int i, j, nxpercent, nypercent, iindex, jindex
+real f
+
+begin
+ nxpercent = apodize * nx
+ nypercent = apodize * ny
+
+ if (forward == YES) {
+ do j = 1, ny {
+ do i = 1, nxpercent {
+ iindex = nx - i + 1
+ f = (1.0 - cos (PI * real (i-1) / real(nxpercent))) / 2.0
+ data[i,j] = f * data[i,j]
+ data[iindex,j] = f * data[iindex,j]
+ }
+ }
+ do i = 1, nx {
+ do j = 1, nypercent {
+ jindex = ny - j + 1
+ f = (1.0 - cos (PI * real (j-1) / real(nypercent))) / 2.0
+ data[i,j] = f * data[i,j]
+ data[i,jindex] = f * data[i,jindex]
+ }
+ }
+ } else {
+ do j = 1, ny {
+ do i = 1, nxpercent {
+ iindex = nx - i + 1
+ f = (1.0 - cos (PI * real (i-1) / real(nxpercent))) / 2.0
+ if (f < 1.0e-3)
+ f = 1.0e-3
+ data[i,j] = data[i,j] / f
+ data[iindex,j] = data[iindex,j] / f
+ }
+ }
+ do i = 1, nx {
+ do j = 1, nypercent {
+ jindex = ny - j + 1
+ f = (1.0 - cos (PI * real (j-1) / real(nypercent))) / 2.0
+ if (f < 1.0e-3)
+ f = 1.0e-3
+ data[i,j] = data[i,j] / f
+ data[i,jindex] = data[i,jindex] / f
+ }
+ }
+ }
+end
+
+
+# RG_ZNSUM -- Compute the mean and number of good points in the array with
+# one optional level of rejections.
+
+int procedure rg_znsum (data, npts, mean, lcut, hcut)
+
+real data[ARB] #I the input data array
+int npts #I the number of data points
+real mean #O the mean of the data
+real lcut, hcut #I the good data limits
+
+int i, ngpts
+real dif, sigma, sum, sumsq, lo, hi
+real asumr(), assqr()
+
+begin
+ # Get the mean.
+ if (npts == 0) {
+ mean = INDEFR
+ return (0)
+ } else if (npts == 1) {
+ mean = data[1]
+ return (1)
+ } else {
+ sum = asumr (data, npts)
+ mean = sum / npts
+ }
+
+ # Quit if the rejection flags are not set.
+ if (IS_INDEFR(lcut) && IS_INDEFR(hcut))
+ return (npts)
+
+ # Compute sigma
+ sumsq = assqr (data, npts)
+ sigma = sumsq / (npts - 1) - mean * sum / (npts - 1)
+ if (sigma <= 0.0)
+ sigma = 0.0
+ else
+ sigma = sqrt (sigma)
+ if (sigma <= 0.0)
+ return (npts)
+
+ # Do the k-sigma rejection.
+ if (IS_INDEF(lcut))
+ lo = -MAX_REAL
+ else
+ lo = -lcut * sigma
+ if (IS_INDEFR(hcut))
+ hi = MAX_REAL
+ else
+ hi = hcut * sigma
+
+ # Reject points.
+ ngpts = npts
+ do i = 1, npts {
+ dif = (data[i] - mean)
+ if (dif >= lo && dif <= hi)
+ next
+ ngpts = ngpts - 1
+ sum = sum - data[i]
+ sumsq = sumsq - data[i] ** 2
+ }
+
+ # Get the final mean.
+ if (ngpts == 0) {
+ mean = INDEFR
+ return (0)
+ } else if (ngpts == 1) {
+ mean = sum
+ return (1)
+ } else
+ mean = sum / ngpts
+
+ return (ngpts)
+end
+
+
+# RG_ZNMEDIAN -- Compute the median and number of good points in the array
+# with one level of rejection.
+
+int procedure rg_znmedian (data, npts, median, lcut, hcut)
+
+real data[ARB] #I the input data array
+int npts #I the number of data points
+real median #O the median of the data
+real lcut, hcut #I the good data limits
+
+int i, ngpts, lindex, hindex
+pointer sp, sdata
+real mean, sigma, dif, lo, hi
+real amedr()
+
+begin
+ if (IS_INDEFR (lcut) && IS_INDEFR(hcut)) {
+ median = amedr (data, npts)
+ return (npts)
+ }
+
+ # Allocate working space.
+ call smark (sp)
+ call salloc (sdata, npts, TY_REAL)
+ call asrtr (data, Memr[sdata], npts)
+ if (mod (npts, 2) == 0)
+ median = (Memr[sdata+(1+npts)/2-1] + Memr[sdata+(1+npts)/2]) / 2.0
+ else
+ median = Memr[sdata+(1+npts)/2-1]
+
+ # Compute the sigma.
+ call aavgr (Memr[sdata], npts, mean, sigma)
+ if (sigma <= 0.0) {
+ call sfree (sp)
+ return (npts)
+ }
+
+ # Do rejection.
+ ngpts = npts
+ if (IS_INDEFR(lo))
+ lo = -MAX_REAL
+ else
+ lo = -lcut * sigma
+ if (IS_INDEFR(hi))
+ hi = MAX_REAL
+ else
+ hi = hcut * sigma
+
+ do i = 1, npts {
+ lindex = i
+ dif = Memr[sdata+i-1] - median
+ if (dif >= lo)
+ break
+ }
+ do i = npts, 1, -1 {
+ hindex = i
+ dif = Memr[sdata+i-1] - median
+ if (dif <= hi)
+ break
+ }
+
+ ngpts = hindex - lindex + 1
+ if (ngpts <= 0)
+ median = INDEFR
+ else if (mod (ngpts, 2) == 0)
+ median = (Memr[sdata+lindex-1+(ngpts+1)/2-1] + Memr[sdata+lindex-1+
+ (ngpts+1)/2]) / 2.0
+ else
+ median = Memr[sdata+lindex-1+(ngpts+1)/2-1]
+
+ call sfree (sp)
+
+ return (ngpts)
+end
+
+
+# RG_SLOPE -- Subtract a slope from the data to be psf matched.
+
+int procedure rg_slope (gs, data, npts, nx, ny, wxborder, wyborder, loreject,
+ hireject)
+
+pointer gs #I the pointer to surfit structure
+real data[ARB] #I/O the input/output data
+int npts #I the number of points
+int nx, ny #I dimensions of the original data
+int wxborder, wyborder #I the x and y width of the border
+real loreject, hireject #I the rejection criteria
+
+int i, stat, ier
+pointer sp, x, y, w, zfit
+real lcut, hcut, sigma
+int rg_reject(), rg_breject()
+real rg_sigma(), rg_bsigma()
+
+begin
+ # Initialize.
+ call smark (sp)
+ call salloc (x, nx, TY_REAL)
+ call salloc (y, nx, TY_REAL)
+ call salloc (w, nx, TY_REAL)
+ call salloc (zfit, nx, TY_REAL)
+ do i = 1, nx
+ Memr[x+i-1] = i
+ call amovkr (1.0, Memr[w], nx)
+
+ # Accumulate the fit.
+ call gszero (gs)
+ if (npts >= nx * ny)
+ call rg_gsaccum (gs, Memr[x], Memr[y], Memr[w], data, nx, ny)
+ else
+ call rg_gsborder (gs, Memr[x], Memr[y], Memr[w], data, nx, ny,
+ wxborder, wyborder)
+
+ # Solve the surface.
+ call gssolve (gs, ier)
+ if (ier == NO_DEG_FREEDOM) {
+ call sfree (sp)
+ return (ERR)
+ }
+
+ # Perform the rejection.
+ if (! IS_INDEFR(loreject) || ! IS_INDEFR(hireject)) {
+ if (npts >= nx * ny)
+ sigma = rg_sigma (gs, Memr[x], Memr[y], Memr[w], Memr[zfit],
+ data, nx, ny)
+ else
+ sigma = rg_bsigma (gs, Memr[x], Memr[y], Memr[w], Memr[zfit],
+ data, nx, ny, wxborder, wyborder)
+ if (sigma <= 0.0) {
+ call sfree (sp)
+ return (OK)
+ }
+ if (! IS_INDEFR(loreject))
+ lcut = -loreject * sigma
+ else
+ lcut = -MAX_REAL
+ if (! IS_INDEFR(hireject))
+ hcut = hireject * sigma
+ else
+ hcut = MAX_REAL
+ if (npts >= nx * ny)
+ stat = rg_reject (gs, Memr[x], Memr[y], Memr[w], Memr[zfit],
+ data, nx, ny, lcut, hcut)
+ else
+ stat = rg_breject (gs, Memr[x], Memr[y], Memr[w], Memr[zfit],
+ data, nx, ny, wxborder, wyborder, lcut, hcut)
+ }
+
+ call sfree (sp)
+ return (stat)
+end
+
+
+# RG_GSACCUM -- Accumulate the points into the fits assuming the data is in the
+# form of a two-dimensional subraster.
+
+procedure rg_gsaccum (gs, x, y, w, data, nx, ny)
+
+pointer gs #I pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the input y array
+real w[ARB] #I the input weight array
+real data[ARB] #I the input data array
+int nx, ny #I the size of the input data array
+
+int i, index
+
+begin
+ index = 1
+ do i = 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsacpts (gs, x, y, data[index], w, nx, WTS_USER)
+ index = index + nx
+ }
+end
+
+
+# RG_GSBORDER -- Procedure to accumulate the points into the fit assuming
+# that a border has been extracted
+
+procedure rg_gsborder (gs, x, y, w, data, nx, ny, wxborder, wyborder)
+
+pointer gs #I pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the input y array
+real w[ARB] #I the input weight array
+real data[ARB] #I the input data array
+int nx, ny #I the dimensions of the input data
+int wxborder, wyborder #I the width of the border
+
+int i, index, nborder
+
+begin
+ nborder = nx * ny - (nx - wxborder) * (ny - wyborder)
+
+ index = 1
+ do i = 1, wyborder {
+ call amovkr (real (i), y, nx)
+ call gsacpts (gs, x, y, data[index], w, nx, WTS_USER)
+ index = index + nx
+ }
+
+ index = nx * wyborder + 1
+ do i = wyborder + 1, ny - wyborder {
+ call amovkr (real (i), y, nx)
+ call gsacpts (gs, x, y, data[index], w, wxborder, WTS_USER)
+ index = index + wxborder
+ call gsacpts (gs, x[1+nx-wxborder], y[1+nx-wxborder],
+ data[index], w[1+nx-wxborder], wxborder, WTS_USER)
+ index = index + wxborder
+ }
+
+ index = 1 + nborder - nx * wyborder
+ do i = ny - wyborder + 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsacpts (gs, x, y, data[index], w, nx, WTS_USER)
+ index = index + nx
+ }
+
+end
+
+
+# RG_SIGMA -- Compute sigma assuming the data is in the form of a
+# two-dimensional subraster.
+
+real procedure rg_sigma (gs, x, y, w, zfit, data, nx, ny)
+
+pointer gs #I the pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the input y array
+real w[ARB] #I the input w array
+real zfit[ARB] #O the output fitted data
+real data[ARB] #I/O the input/output data array
+int nx, ny #I the dimensions of the output data
+
+int i, j, index, npts
+real sum
+
+begin
+ npts = 0
+ index = 1
+ sum = 0.0
+
+ do i = 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (w[j] > 0.0) {
+ sum = sum + zfit[j] ** 2
+ npts = npts + 1
+ }
+ }
+ index = index + nx
+ }
+
+ return (sqrt (sum / npts))
+end
+
+
+# RG_BSIGMA -- Procedure to compute sigma assuming a border has been
+# extracted from a subraster.
+
+real procedure rg_bsigma (gs, x, y, w, zfit, data, nx, ny, wxborder, wyborder)
+
+pointer gs #I the pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the output y array
+real w[ARB] #I the output weight array
+real zfit[ARB] #O the fitted z array
+real data[ARB] #I/O the input/output data array
+int nx, ny #I the dimensions of original subraster
+int wxborder, wyborder #I the width of the border
+
+int i, j, npts, nborder, index
+real sum
+
+begin
+ nborder = nx * ny - (nx - wxborder) * (ny - wyborder)
+ npts = 0
+ index = 1
+ sum = 0.0
+
+ do i = 1, wyborder {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (w[j] > 0.0) {
+ npts = npts + 1
+ sum = sum + zfit[j] ** 2
+ }
+ }
+ index = index + nx
+ }
+
+ index = nx * wyborder + 1
+ do i = wyborder + 1, ny - wyborder {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, wxborder)
+ call asubr (data[index], zfit, zfit, wxborder)
+ do j = 1, wxborder {
+ if (w[j] > 0.0) {
+ npts = npts + 1
+ sum = sum + zfit[j] ** 2
+ }
+ }
+ index = index + wxborder
+ call gsvector (gs, x[1+nx-wxborder], y[1+nx-wxborder], zfit,
+ wxborder)
+ call asubr (data[index], zfit, zfit, wxborder)
+ do j = 1, wxborder {
+ if (w[j] > 0.0) {
+ npts = npts + 1
+ sum = sum + zfit[j] ** 2
+ }
+ }
+ index = index + wxborder
+ }
+
+ index = 1 + nborder - nx * wyborder
+ do i = ny - wyborder + 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (w[j] > 0.0) {
+ npts = npts + 1
+ sum = sum + zfit[j] ** 2
+ }
+ }
+ index = index + nx
+ }
+
+ return (sqrt (sum / npts))
+end
+
+
+# RG_REJECT -- Reject points from the fit assuming the data is in the form of a
+# two-dimensional subraster.
+
+int procedure rg_reject (gs, x, y, w, zfit, data, nx, ny, lcut, hcut)
+
+pointer gs #I the pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the input y array
+real w[ARB] #I the input w array
+real zfit[ARB] #O the fitted data
+real data[ARB] #I/O the input/output data array
+int nx, ny #I the dimensions of the data
+real lcut, hcut #I the lo and high side rejection criteria
+
+int i, j, index, ier
+
+begin
+ index = 1
+
+ do i = 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (zfit[j] < lcut || zfit[j] > hcut)
+ call gsrej (gs, x[j], y[j], data[index+j-1], w[j], WTS_USER)
+ }
+ index = index + nx
+ }
+
+ call gssolve (gs, ier)
+ if (ier == NO_DEG_FREEDOM)
+ return (ERR)
+ else
+ return (OK)
+end
+
+
+# RG_BREJECT -- Reject deviant points from the fits assuming a border has
+# been extracted from the subraster.
+
+int procedure rg_breject (gs, x, y, w, zfit, data, nx, ny, wxborder,
+ wyborder, lcut, hcut)
+
+pointer gs #I the pointer to the surface fitting structure
+real x[ARB] #I the input x array
+real y[ARB] #I the input y array
+real w[ARB] #I the input weight array
+real zfit[ARB] #O the fitted z array
+real data[ARB] #I/O the input/output data array
+int nx, ny #I the dimensions of the original subraster
+int wxborder, wyborder #I the width of the border
+real lcut, hcut #I the low and high rejection criteria
+
+int i, j, nborder, index, ier
+
+begin
+ nborder = nx * ny - (nx - wxborder) * (ny - wyborder)
+ index = 1
+
+ do i = 1, wyborder {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (zfit[j] < lcut || zfit[j] > hcut)
+ call gsrej (gs, x[j], y[j], data[index+j-1], w[j],
+ WTS_USER)
+ }
+ index = index + nx
+ }
+
+ index = nx * wyborder + 1
+ do i = wyborder + 1, ny - wyborder {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, wxborder)
+ call asubr (data[index], zfit, zfit, wxborder)
+ do j = 1, wxborder {
+ if (zfit[j] < lcut || zfit[j] > hcut)
+ call gsrej (gs, x[j], y[j], data[index+j-1], w[j],
+ WTS_USER)
+ }
+ index = index + wxborder
+ call gsvector (gs, x[1+nx-wxborder], y[1+nx-wxborder], zfit,
+ wxborder)
+ call asubr (data[index], zfit, zfit, wxborder)
+ do j = 1, wxborder {
+ if (zfit[j] < lcut || zfit[j] > hcut)
+ call gsrej (gs, x[j], y[j], data[index+j-1], w[j],
+ WTS_USER)
+ }
+ index = index + wxborder
+ }
+
+ index = 1 + nborder - nx * wyborder
+ do i = ny - wyborder + 1, ny {
+ call amovkr (real (i), y, nx)
+ call gsvector (gs, x, y, zfit, nx)
+ call asubr (data[index], zfit, zfit, nx)
+ do j = 1, nx {
+ if (zfit[j] < lcut || zfit[j] > hcut)
+ call gsrej (gs, x[j], y[j], data[index+j-1], w[j],
+ WTS_USER)
+ }
+ index = index + nx
+ }
+
+ call gssolve (gs, ier)
+ if (ier == NO_DEG_FREEDOM)
+ return (ERR)
+ else
+ return (OK)
+end
+
generated by cgit (git 2.34.1) at 2025-09-20 05:24:37 -0400