diff options
Diffstat (limited to 'pkg/images/imfit/src')
| -rw-r--r-- | pkg/images/imfit/src/fit1d.x | 597 | ||||
| -rw-r--r-- | pkg/images/imfit/src/imsurfit.h | 40 | ||||
| -rw-r--r-- | pkg/images/imfit/src/imsurfit.x | 1172 | ||||
| -rw-r--r-- | pkg/images/imfit/src/mkpkg | 15 | ||||
| -rw-r--r-- | pkg/images/imfit/src/pixlist.h | 11 | ||||
| -rw-r--r-- | pkg/images/imfit/src/pixlist.x | 369 | ||||
| -rw-r--r-- | pkg/images/imfit/src/ranges.x | 524 | ||||
| -rw-r--r-- | pkg/images/imfit/src/t_imsurfit.x | 400 | ||||
| -rw-r--r-- | pkg/images/imfit/src/t_lineclean.x | 270 |
9 files changed, 3398 insertions, 0 deletions
diff --git a/pkg/images/imfit/src/fit1d.x b/pkg/images/imfit/src/fit1d.x new file mode 100644 index 00000000..84ffddf1 --- /dev/null +++ b/pkg/images/imfit/src/fit1d.x @@ -0,0 +1,597 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <pkg/gtools.h> +include <error.h> + +define MAXBUF (512*100) # Maximum number of pixels per block + + +# FIT1D -- Fit a function to image lines or columns and output an image +# consisting of the fit, the difference, or the ratio. The fitting parameters +# may be set interactively using the icfit package. + +procedure t_fit1d () + +int listin # Input image list +int listout # Output image list +int listbpm # Bad pixel mask list +bool interactive # Interactive? + +char sample[SZ_LINE] # Sample ranges +int naverage # Sample averaging size +char function[SZ_LINE] # Curve fitting function +int order # Order of curve fitting function +real low_reject, high_reject # Rejection thresholds +int niterate # Number of rejection iterations +real grow # Rejection growing radius + +int axis # Image axis to fit +int ntype # Type of output +char input[SZ_LINE] # Input image +char output[SZ_FNAME] # Output image +char bpm[SZ_FNAME] # Bad pixel mask +pointer in, out, bp # IMIO pointers +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer + +bool same, clgetb() +int imtopen(), imtgetim(), imtlen(), strdic(), gt_init() +int clgeti() +real clgetr() + +begin + # Get input and output lists and check that the number of images + # are the same. + + call clgstr ("input", input, SZ_LINE) + listin = imtopen (input) + call clgstr ("output", input, SZ_LINE) + listout = imtopen (input) + if (imtlen (listin) != imtlen (listout)) { + call imtclose (listin) + call imtclose (listout) + call error (0, "Input and output image lists do not match") + } + call clgstr ("bpm", input, SZ_LINE) + listbpm = imtopen (input) + if (imtlen (listbpm) > 1 && imtlen (listin) != imtlen (listbpm)) { + call imtclose (listin) + call imtclose (listout) + call imtclose (listbpm) + call error (0, "Input and mask lists do not match") + } + + # Get task parameters. + + axis = clgeti ("axis") + call clgstr ("type", input, SZ_LINE) + call clgstr ("sample", sample, SZ_LINE) + naverage = clgeti ("naverage") + call clgstr ("function", function, SZ_LINE) + order = clgeti ("order") + low_reject = clgetr ("low_reject") + high_reject = clgetr ("high_reject") + niterate = clgeti ("niterate") + grow = clgetr ("grow") + interactive = clgetb ("interactive") + + # Decode the output type and initialize the curve fitting package. + + ntype = strdic (input, input, SZ_LINE, "|fit|difference|ratio|") + if (ntype == 0) + call error (0, "Unknown output type") + + # Set the ICFIT pointer structure. + call ic_open (ic) + call ic_pstr (ic, "sample", sample) + call ic_puti (ic, "naverage", naverage) + call ic_pstr (ic, "function", function) + call ic_puti (ic, "order", order) + call ic_putr (ic, "low", low_reject) + call ic_putr (ic, "high", high_reject) + call ic_puti (ic, "niterate", niterate) + call ic_putr (ic, "grow", grow) + call ic_pstr (ic, "ylabel", "") + + gt = gt_init() + call gt_sets (gt, GTTYPE, "line") + + # Fit the lines in each input image. + + bpm[1] = EOS + while ((imtgetim (listin, input, SZ_LINE) != EOF) && + (imtgetim (listout, output, SZ_FNAME) != EOF)) { + if (imtgetim (listbpm, bpm, SZ_FNAME) == EOF) + ; + + iferr (call f1d_immap (input,output,bpm,ntype,in,out,bp,same)) { + call erract (EA_WARN) + next + } + call f1d_fit1d (in,out,bp,ic,gt,input,axis,ntype,interactive) + call imunmap (in) + if (!same) + call imunmap (out) + if (bp != NULL) + call yt_pmunmap (bp) + } + + call ic_closer (ic) + call gt_free (gt) + call imtclose (listin) + call imtclose (listout) +end + + +# F1D_FIT1D -- Given the image descriptor determine the fitting function +# for each line or column and create an output image. If the interactive flag +# is set then set the fitting parameters interactively. + +procedure f1d_fit1d (in, out, bp, ic, gt, title, axis, ntype, interactive) + +pointer in # IMIO pointer for input image +pointer out # IMIO pointer for output image +pointer bp # IMIO pointer for bad pixel mask +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer +char title[ARB] # Title +int axis # Image axis to fit +int ntype # Type of output +bool interactive # Interactive? + +char graphics[SZ_FNAME] # Graphics device +int i, nx, new +real div +pointer cv, gp, sp, x, wts, indata, outdata + +int f1d_getline(), f1d_getdata(), strlen() +real cveval() +pointer gopen() + +begin + # Error check. + + if (IM_NDIM (in) > 2) + call error (0, "Image dimensions > 2 are not implemented") + if (axis > IM_NDIM (in)) + call error (0, "Axis exceeds image dimension") + + # Allocate memory for curve fitting. + + nx = IM_LEN (in, axis) + call smark (sp) + call salloc (x, nx, TY_REAL) + + do i = 1, nx + Memr[x+i-1] = i + + call ic_putr (ic, "xmin", Memr[x]) + call ic_putr (ic, "xmax", Memr[x+nx-1]) + + # If the interactive flag is set then use icg_fit to set the + # fitting parameters. Get_fitline returns EOF when the user + # is done. The weights are reset since the user may delete + # points. + + if (interactive) { + call clgstr ("graphics", graphics, SZ_FNAME) + gp = gopen (graphics, NEW_FILE, STDGRAPH) + + i = strlen (title) + indata = NULL + while (f1d_getline (ic,gt,in,bp,axis,title,indata,wts) != EOF) { + title[i + 1] = EOS + call icg_fit (ic, gp, "cursor", gt, cv, Memr[x], Memr[indata], + Memr[wts], nx) + } + call mfree (indata, TY_REAL) + call mfree (wts, TY_REAL) + call gclose (gp) + } + + # Loop through the input image and create an output image. + + new = YES + + while (f1d_getdata (in,out,bp,axis,MAXBUF,indata,outdata,wts) != EOF) { + + call ic_fit (ic, cv, Memr[x], Memr[indata], Memr[wts], + nx, new, YES, new, new) + new = NO + + # Be careful because the indata and outdata buffers may be the same. + switch (ntype) { + case 1: + call cvvector (cv, Memr[x], Memr[outdata], nx) + case 2: + do i = 0, nx-1 + Memr[outdata+i] = Memr[indata+i] - cveval (cv, Memr[x+i]) + case 3: + do i = 0, nx-1 { + div = cveval (cv, Memr[x+i]) + if (abs (div) < 1E-20) + div = 1 + Memr[outdata+i] = Memr[indata+i] / div + } + } + } + + call cvfree (cv) + call mfree (wts, TY_REAL) + call sfree (sp) +end + + +# F1D_IMMAP -- Map images for fit1d. + +procedure f1d_immap (input, output, bpm, ntype, in, out, bp, same) + +char input[ARB] # Input image +char output[ARB] # Output image +char bpm[ARB] # Bad pixel mask +int ntype # Type of fit1d output +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +pointer bp # Mask IMIO pointer +bool same # Same image? + +int i +pointer sp, iroot, isect, oroot, osect, line, data + +bool streq() +int imaccess(), impnlr() +pointer immap(), yt_pmmap() +errchk immap, yt_pmmap + +begin + # Get the root name and section of the input image. + + call smark (sp) + call salloc (iroot, SZ_FNAME, TY_CHAR) + call salloc (isect, SZ_FNAME, TY_CHAR) + call salloc (oroot, SZ_FNAME, TY_CHAR) + call salloc (osect, SZ_FNAME, TY_CHAR) + + call imgimage (input, Memc[iroot], SZ_FNAME) + call imgsection (input, Memc[isect], SZ_FNAME) + call imgimage (output, Memc[oroot], SZ_FNAME) + call imgsection (output, Memc[osect], SZ_FNAME) + same = streq (Memc[iroot], Memc[oroot]) + + # If the output image is not accessible then create it as a new copy + # of the full input image and initialize according to ntype. + + if (imaccess (output, READ_WRITE) == NO) { + in = immap (Memc[iroot], READ_ONLY, 0) + out = immap (Memc[oroot], NEW_COPY, in) + IM_PIXTYPE(out) = TY_REAL + + call salloc (line, IM_MAXDIM, TY_LONG) + call amovkl (long (1), Meml[line], IM_MAXDIM) + + switch (ntype) { + case 1, 2: + while (impnlr (out, data, Meml[line]) != EOF) + call aclrr (Memr[data], IM_LEN(out, 1)) + case 3: + while (impnlr (out, data, Meml[line]) != EOF) + call amovkr (1., Memr[data], IM_LEN(out, 1)) + } + + call imunmap (in) + call imunmap (out) + } + + # Map the images. If the output image has a section + # then use it. If the input image has a section and the output image + # does not then add the image section to the output image. Finally + # check the input and output images have the same size. + + in = immap (input, READ_ONLY, 0) + + if (Memc[isect] != EOS && Memc[osect] == EOS) { + call sprintf (Memc[osect], SZ_FNAME, "%s%s") + call pargstr (Memc[oroot]) + call pargstr (Memc[isect]) + } else + call strcpy (output, Memc[osect], SZ_FNAME) + + if (streq (input, Memc[osect])) { + call imunmap (in) + in = immap (input, READ_WRITE, 0) + out = in + } else + out = immap (Memc[osect], READ_WRITE, 0) + + do i = 1, IM_NDIM(in) + if (IM_LEN(in, i) != IM_LEN(out, i)) { + call imunmap (in) + if (!same) + call imunmap (out) + call sfree (sp) + call error (0, "Input and output images have different sizes") + } + + bp = yt_pmmap (bpm, in, Memc[iroot], SZ_FNAME) + + call sfree (sp) +end + + +# F1D_GETDATA -- Get a line of image data. + +int procedure f1d_getdata (in, out, bp, axis, maxbuf, indata, outdata, wts) + +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer +pointer bp # Bad pixel mask IMIO pointer +int axis # Image axis +int maxbuf # Maximum buffer size for column axis +pointer indata # Input data pointer +pointer outdata # Output data pointer +pointer wts # Weights pointer + +int i, index, last_index, col1, col2, nc, ncols, nlines, ncols_block +pointer inbuf, outbuf, wtbuf, ptr +pointer imgl1r(), imgl1s(), impl1r(), imgl2r(), imgl2s(), impl2r() +pointer imgs2r(), imgs2s(), imps2r() +data index/0/ + +begin + # Increment to the next image vector. + + index = index + 1 + + # Initialize for the first vector. + + if (index == 1) { + ncols = IM_LEN (in, 1) + if (IM_NDIM (in) == 1) + nlines = 1 + else + nlines = IM_LEN (in, 2) + + switch (axis) { + case 1: + last_index = nlines + + call malloc (wts, ncols, TY_REAL) + case 2: + last_index = ncols + ncols_block = max (1, min (ncols, maxbuf / nlines)) + col2 = 0 + + call malloc (indata, nlines, TY_REAL) + call malloc (outdata, nlines, TY_REAL) + call malloc (wts, nlines, TY_REAL) + } + } + + # Finish up if the last vector has been done. + + if (index > last_index) { + switch (axis) { + case 1: + call mfree (wts, TY_REAL) + case 2: + ptr = outbuf + index - 1 - col1 + do i = 1, nlines { + Memr[ptr] = Memr[outdata+i-1] + ptr = ptr + nc + } + + call mfree (indata, TY_REAL) + call mfree (outdata, TY_REAL) + call mfree (wts, TY_REAL) + } + + index = 0 + return (EOF) + } + + # Get the next image vector. + + switch (axis) { + case 1: + ncols = IM_LEN(in,1) + if (IM_NDIM (in) == 1) { + indata = imgl1r (in) + outdata = impl1r (out) + if (bp == NULL) + call amovkr (1., Memr[wts], ncols) + else { + wtbuf = imgl1s (bp) + do i = 0, ncols-1 { + if (Mems[wtbuf+i] == 0) + Memr[wts+i] = 1. + else + Memr[wts+i] = 0. + } + } + } else { + indata = imgl2r (in, index) + outdata = impl2r (out, index) + if (bp == NULL) + call amovkr (1., Memr[wts], ncols) + else { + wtbuf = imgl2s (bp, index) + do i = 0, ncols-1 { + if (Mems[wtbuf+i] == 0) + Memr[wts+i] = 1. + else + Memr[wts+i] = 0. + } + } + } + case 2: + if (index > 1) { + ptr = outbuf + index - 1 - col1 + do i = 1, nlines { + Memr[ptr] = Memr[outdata+i-1] + ptr = ptr + nc + } + } + + if (index > col2) { + col1 = col2 + 1 + col2 = min (ncols, col1 + ncols_block - 1) + inbuf = imgs2r (in, col1, col2, 1, nlines) + outbuf = imps2r (out, col1, col2, 1, nlines) + if (bp != NULL) + wtbuf = imgs2s (bp, col1, col2, 1, nlines) + nc = col2 - col1 + 1 + } + + ptr = inbuf + index - col1 + do i = 0, nlines-1 { + Memr[indata+i] = Memr[ptr] + ptr = ptr + nc + } + + if (bp == NULL) + call amovkr (1., Memr[wts], nlines) + else { + ptr = wtbuf + index - col1 + do i = 0, nlines-1 { + if (Mems[ptr] == 0) + Memr[wts+i] = 1. + else + Memr[wts+i] = 0. + ptr = ptr + nc + } + } + } + + return (index) +end + + +# F1D_GETLINE -- Get image data to be fit interactively. Return EOF +# when the user enters EOF or CR. Default is 1 and the out of bounds +# requests are silently limited to the nearest in edge. + +int procedure f1d_getline (ic, gt, im, bp, axis, title, data, wts) + +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer +pointer im # IMIO pointer input image +pointer bp # IMIO pointer for bad pixel mask +int axis # Image axis +char title[ARB] # Title +pointer data # Image data +pointer wts # Weights + +pointer x, wtbuf +char line[SZ_LINE] +int i, j, stat, imlen +int getline(), nscan() +pointer imgl1r(), imgl1s() +data stat/EOF/ + +begin + # If the image is one dimensional do not prompt. + + if (IM_NDIM (im) == 1) { + if (stat == EOF) { + call sprintf (title, SZ_LINE, "%s\n%s") + call pargstr (title) + call pargstr (IM_TITLE(im)) + call gt_sets (gt, GTTITLE, title) + call mfree (data, TY_REAL) + imlen = IM_LEN(im,1) + call malloc (data, imlen, TY_REAL) + call amovr (Memr[imgl1r(im)], Memr[data], imlen) + call malloc (wts, imlen, TY_REAL) + if (bp == NULL) + call amovkr (1., Memr[wts], imlen) + else { + wtbuf = imgl1s (bp) + do i = 0, imlen-1 { + if (Mems[wtbuf+i] == 0) + Memr[wts+i] = 1. + else + Memr[wts+i] = 0. + } + } + stat = OK + } else + stat = EOF + + return (stat) + } + + # If the image is two dimensional prompt for the line or column. + + switch (axis) { + case 1: + imlen = IM_LEN (im, 2) + call sprintf (title, SZ_LINE, "%s: Fit line =") + call pargstr (title) + case 2: + imlen = IM_LEN (im, 1) + call sprintf (title, SZ_LINE, "%s: Fit column =") + call pargstr (title) + } + + call printf ("%s ") + call pargstr (title) + call flush (STDOUT) + + if (getline(STDIN, line) == EOF) + return (EOF) + + call sscan (line) + call gargi (i) + call gargi (j) + + switch (nscan()) { + case 0: + stat = EOF + return (stat) + case 1: + i = max (1, min (imlen, i)) + j = i + case 2: + i = max (1, min (imlen, i)) + j = max (1, min (imlen, j)) + } + + call sprintf (title, SZ_LINE, "%s %d - %d\n%s") + call pargstr (title) + call pargi (i) + call pargi (j) + call pargstr (IM_TITLE(im)) + + call gt_sets (gt, GTTITLE, title) + + call malloc (wts, imlen, TY_REAL) + switch (axis) { + case 1: + call ic_pstr (ic, "xlabel", "Column") + call xt_21imavg (im, axis, 1, IM_LEN(im,1), i, j, x, data, imlen) + if (bp != NULL) + call xt_21imsum (bp, axis, 1, IM_LEN(im,1), i, j, x, wts, imlen) + case 2: + call ic_pstr (ic, "xlabel", "Line") + call xt_21imavg (im, axis, i, j, 1, IM_LEN(im,2), x, data, imlen) + if (bp != NULL) + call xt_21imsum (bp, axis, i, j, 1, IM_LEN(im,2), x, wts, imlen) + } + if (bp == NULL) { + call mfree (wts, TY_REAL) + call malloc (wts, imlen, TY_REAL) + call amovkr (1., Memr[wts], imlen) + } else { + do i = 0, imlen-1 { + if (Memr[wts+i] == 0.) + Memr[wts+i] = 1. + else + Memr[wts+i] = 0. + } + } + call mfree (x, TY_REAL) + + stat = OK + return (stat) +end diff --git a/pkg/images/imfit/src/imsurfit.h b/pkg/images/imfit/src/imsurfit.h new file mode 100644 index 00000000..84c077ec --- /dev/null +++ b/pkg/images/imfit/src/imsurfit.h @@ -0,0 +1,40 @@ +# Header file for IMSURFIT + +define LEN_IMSFSTRUCT 20 + +# surface parameters +define SURFACE_TYPE Memi[$1] +define XORDER Memi[$1+1] +define YORDER Memi[$1+2] +define CROSS_TERMS Memi[$1+3] +define TYPE_OUTPUT Memi[$1+4] + +# median processing parameters +define MEDIAN Memi[$1+5] +define XMEDIAN Memi[$1+6] +define YMEDIAN Memi[$1+7] +define MEDIAN_PERCENT Memr[P2R($1+8)] + +# pixel rejection parameters +define REJECT Memi[$1+9] +define NGROW Memi[$1+10] +define NITER Memi[$1+11] +define LOWER Memr[P2R($1+12)] +define UPPER Memr[P2R($1+13)] + +define DIV_MIN Memr[P2R($1+14)] + +# definitions for type_output +define FIT 1 +define CLEAN 2 +define RESID 3 +define RESP 4 + +# definitions for good regions +define ALL 1 +define COLUMNS 2 +define ROWS 3 +define BORDER 4 +define SECTIONS 5 +define CIRCLE 6 +define INVCIRCLE 7 diff --git a/pkg/images/imfit/src/imsurfit.x b/pkg/images/imfit/src/imsurfit.x new file mode 100644 index 00000000..9f655f52 --- /dev/null +++ b/pkg/images/imfit/src/imsurfit.x @@ -0,0 +1,1172 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <imhdr.h> +include <math/surfit.h> +include "imsurfit.h" + +# IMSURFIT -- Procedure to fit a surface to a single image including +# optional pixel rejection. + +procedure imsurfit (imin, imout, imfit, gl) + +pointer imin # pointer to the input image +pointer imout # pointer to the output image +pointer imfit # pointer to the imsurfit parameters +pointer gl # pointer to the good regions list + +pointer sf, rl +errchk isfree, prl_free +errchk all_pixels, good_pixels, good_median, all_medians, do_reject +errchk set_outimage + +begin + sf = NULL + rl = NULL + + # Accumulate and solve the surface. + if (gl == NULL) { + if (MEDIAN(imfit) == NO) + call all_pixels (imin, imfit, sf) + else + call all_medians (imin, imfit, sf) + } else { + if (MEDIAN(imfit) == NO) + call good_pixels (imin, imfit, gl, sf) + else + call good_medians (imin, imfit, gl, sf) + } + + # Perform the reject cycle. + if (REJECT(imfit) == YES || TYPE_OUTPUT(imfit) == CLEAN) + call do_reject (imin, imfit, gl, sf, rl) + + # Evaluate surface for appropriate output type. + call set_outimage (imin, imout, imfit, sf, rl) + + # Cleanup. + call prl_free (rl) + call isfree (sf) + + rl = NULL + sf = NULL +end + + +# ALL_PIXELS -- Accumulate surface when there are no bad regions +# and no median processing. + +procedure all_pixels (im, imfit, sf) + +pointer im # pointer to the input image +pointer imfit # pointer to the imsurfit structure +pointer sf # pointer to the surface descriptor + +int i, lp, ncols, nlines, ier +long v[IM_MAXDIM] +pointer sp, cols, lines, wgt, lbuf +int imgnlr() +errchk smark, salloc, sfree, imgnlr +errchk isinit, islfit, islrefit, issolve + +begin + ncols = IM_LEN(im, 1) + nlines = IM_LEN(im,2) + + # Initialize the surface fit. + call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit), + CROSS_TERMS(imfit), ncols, nlines) + + # Allocate working space for fitting. + call smark (sp) + call salloc (cols, ncols, TY_INT) + call salloc (lines, nlines, TY_INT) + call salloc (wgt, ncols, TY_REAL) + + # Initialize the x and weight buffers. + do i = 1, ncols + Memi[cols - 1 + i] = i + call amovkr (1.0, Memr[wgt], ncols) + + # Loop over image lines. + lp = 0 + call amovkl (long (1), v, IM_MAXDIM) + do i = 1, nlines { + + # Read in the image line. + if (imgnlr (im, lbuf, v) == EOF) + call error (0, "Error reading image") + + # Fit each image line. + if (i == 1) + call islfit (sf, Memi[cols], i, Memr[lbuf], Memr[wgt], + ncols, SF_USER, ier) + else + call islrefit (sf, Memi[cols], i, Memr[lbuf], Memr[wgt]) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call eprintf ("Warning: Too few columns to fit line: %d\n") + call pargi (i) + case SINGULAR: + call eprintf ("Warning: Solution singular for line: %d\n") + call pargi (i) + Memi[lines + lp] = i + lp = lp + 1 + default: + Memi[lines + lp] = i + lp = lp + 1 + } + + } + + # Solve the surface. + call issolve (sf, Memi[lines], lp, ier) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call error (0, "ALL_PIXELS: Cannot fit surface.") + case SINGULAR: + call eprintf ("Warning: Solution singular for surface.\n") + default: + # everything OK + } + + # Free space. + call sfree (sp) +end + + +# GOOD_PIXELS -- Get surface when good regions are defined and median +# processing is off. + +procedure good_pixels (im, imfit, gl, sf) + +pointer im # input image +pointer imfit # pointer to imsurfit header structure +pointer gl # pointer to good region list +pointer sf # pointer to the surface descriptor + +int lp, lineno, prevlineno, ncols, nlines, npts, nranges, ier, ijunk +int max_nranges +pointer sp, colsfit, lines, buf, fbuf, wgt, ranges +int prl_nextlineno(), prl_eqlines(), prl_get_ranges(), is_expand_ranges() +int is_choose_rangesr() +pointer imgl2r() + +errchk smark, salloc, sfree, imgl2r +errchk isinit, islfit, islrefit, issolve +errchk prl_nextlineno, prl_eqlines, prl_get_ranges +errchk is_choose_rangesr + +begin + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + max_nranges = ncols + + # Initialize the surface fit. + call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit), + CROSS_TERMS(imfit), ncols, nlines) + + # Allocate temporary space for fitting. + call smark (sp) + call salloc (colsfit, ncols, TY_INT) + call salloc (lines, nlines, TY_INT) + call salloc (fbuf, ncols, TY_REAL) + call salloc (wgt, ncols, TY_REAL) + call salloc (ranges, 3 * max_nranges + 1, TY_INT) + call amovkr (1., Memr[wgt], ncols) + + # Intialize counters and pointers. + lp = 0 + lineno = 0 + prevlineno = 0 + + # Loop over those lines to be fit. + while (prl_nextlineno (gl, lineno) != EOF) { + + # Read in the image line. + buf = imgl2r (im, lineno) + if (buf == EOF) + call error (0, "GOOD_PIXELS: Error reading image.") + + # Get the ranges for that image line. + nranges = prl_get_ranges (gl, lineno, Memi[ranges], max_nranges) + if (nranges == 0) + next + + # If ranges are not equal to previous line fit else refit. + if (lp == 0 || prl_eqlines (gl, lineno, prevlineno) == NO) { + npts = is_expand_ranges (Memi[ranges], Memi[colsfit], ncols) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf], + npts, 1, ncols) + call islfit (sf, Memi[colsfit], lineno, Memr[fbuf], Memr[wgt], + npts, SF_USER, ier) + } else { + ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf], + npts, 1, ncols) + call islrefit (sf, Memi[colsfit], lineno, Memr[fbuf], Memr[wgt]) + } + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call eprintf ("Warning: Too few columns to fit line: %d\n") + call pargi (lineno) + case SINGULAR: + call eprintf ("Warning: Solution singular for line: %d\n") + call pargi (lineno) + Memi[lines+lp] = lineno + lp = lp + 1 + default: + Memi[lines+lp] = lineno + lp = lp + 1 + } + + prevlineno = lineno + } + + # Solve the surface. + call issolve (sf, Memi[lines], lp, ier) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call error (0, "GOOD_PIXELS: Cannot fit surface.") + case SINGULAR: + call eprintf ("Warning: Solution singular for surface.\n") + default: + # everything OK + } + + # Free space. + call sfree (sp) +end + + +# ALL_MEDIANS -- Get surface when median processor on and all +# pixels good. + +procedure all_medians (im, imfit, sf) + +pointer im # input image +pointer imfit # pointer to the imsurfit header structure +pointer sf # pointer to the surface descriptor + +int i, lp, cp, op, lineno, x1, x2, y1, y2, ier +int nimcols, nimlines, ncols, nlines, npts +pointer sp, cols, lines, wgt, z, med, sbuf, lbuf, buf + +pointer imgs2r() +real asokr() +errchk salloc, sfree, smark +errchk isinit, islfit, islrefit, issolve + +begin + # Determine the number of lines and columns for a median processed + # image. + nimcols = IM_LEN(im,1) + if (mod (int (IM_LEN(im,1)), XMEDIAN(imfit)) != 0) + ncols = IM_LEN(im,1) / XMEDIAN(imfit) + 1 + else + ncols = IM_LEN(im,1) / XMEDIAN(imfit) + nimlines = IM_LEN(im,2) + if (mod (int (IM_LEN(im,2)), YMEDIAN(imfit)) != 0) + nlines = IM_LEN(im,2) / YMEDIAN(imfit) + 1 + else + nlines = IM_LEN(im,2) / YMEDIAN(imfit) + + # Initialize the surface fitting. + call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit), + CROSS_TERMS(imfit), ncols, nlines) + + # Allocate workin memory. + call smark (sp) + call salloc (cols, ncols, TY_INT) + call salloc (wgt, ncols, TY_REAL) + call salloc (lines, nlines, TY_INT) + call salloc (z, ncols, TY_REAL) + call salloc (med, XMEDIAN(imfit) * YMEDIAN(imfit), TY_REAL) + + # Intialize the x and weight arrays. + do i = 1, ncols + Memi[cols - 1 + i] = i + call amovkr (1.0, Memr[wgt], ncols) + + # Loop over image sections. + lp = 0 + lineno = 1 + for (y1 = 1; y1 <= nimlines; y1 = y1 + YMEDIAN(imfit)) { + + # Get image section. + y2 = min (y1 + YMEDIAN(imfit) - 1, nimlines) + sbuf = imgs2r (im, 1, nimcols, y1, y2) + if (sbuf == EOF) + call error (0, "Error reading image section.") + + # Loop over median boxes. + cp = 0 + for (x1 = 1; x1 <= nimcols; x1 = x1 + XMEDIAN(imfit)) { + + x2 = min (x1 + XMEDIAN(imfit) - 1, nimcols) + npts = x2 - x1 + 1 + lbuf = sbuf - 1 + x1 + + # Loop over lines in the median box. + op = 0 + buf = lbuf + for (i = 1; i <= y2 - y1 + 1; i = i + 1) { + call amovr (Memr[buf], Memr[med+op], npts) + op = op + npts + buf = buf + nimcols + } + + # Calculate the median. + Memr[z+cp] = asokr (Memr[med], op, (op + 1) / 2) + cp = cp + 1 + + } + + # Fit each image "line". + if (y1 == 1) + call islfit (sf, Memi[cols], lineno, Memr[z], Memr[wgt], + ncols, SF_USER, ier) + else + call islrefit (sf, Memi[cols], lineno, Memr[z], Memr[wgt]) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call eprintf ("Warning: Too few columns to fit line: %d\n") + call pargi (lineno) + case SINGULAR: + call eprintf ("Warning: Solution singular for line: %d\n") + call pargi (lineno) + Memi[lines + lp] = lineno + lp = lp + 1 + default: + Memi[lines + lp] = lineno + lp = lp + 1 + } + + lineno = lineno + 1 + } + + # Solve th surface. + call issolve (sf, Memi[lines], lp, ier) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call error (0, "ALL_MEDIANS: Cannot fit surface.") + case SINGULAR: + call eprintf ("Warning: Solution singular for surface.\n") + default: + # everything OK + } + + # Free space + call sfree (sp) +end + + +# GOOD_MEDIANS -- Procedure to fetch medians when the good regions +# list is defined. + +procedure good_medians (im, imfit, gl, sf) + +pointer im # input image +pointer imfit # pointer to surface descriptor structure +pointer gl # pointer to good regions list +pointer sf # pointer the surface descriptor + +int i, cp, lp, x1, x2, y1, y2, ier, ntemp +int nimcols, nimlines, ncols, nlines, nranges, nbox, nxpts +int lineno, current_line, lines_per_box, max_nranges +pointer sp, colsfit, cols, lines, wgt, npts, lbuf, med, mbuf, z, ranges + +int prl_get_ranges(), prl_nextlineno(), is_expand_ranges() +int is_choose_rangesr() +pointer imgl2r() +real asokr() +errchk smark, salloc, sfree, imgl2r +errchk isinit, islfit, issolve +errchk prl_get_ranges, prl_nextlineno, is_choose_rangesr() + +begin + # Determine the number of lines and columns for a median processed + # image. + nimcols = IM_LEN(im,1) + if (mod (int (IM_LEN(im,1)), XMEDIAN(imfit)) != 0) + ncols = IM_LEN(im,1) / XMEDIAN(imfit) + 1 + else + ncols = IM_LEN(im,1) / XMEDIAN(imfit) + nimlines = IM_LEN(im,2) + if (mod (int (IM_LEN(im,2)), YMEDIAN(imfit)) != 0) + nlines = IM_LEN(im,2) / YMEDIAN(imfit) + 1 + else + nlines = IM_LEN(im,2) / YMEDIAN(imfit) + nbox = XMEDIAN(imfit) * YMEDIAN(imfit) + max_nranges = nimcols + + # Initialize the surface fitting. + call isinit (sf, SURFACE_TYPE(imfit), XORDER(imfit), YORDER(imfit), + CROSS_TERMS(imfit), ncols, nlines) + + # Allocate working memory. + call smark (sp) + call salloc (colsfit, nimcols, TY_INT) + call salloc (cols, ncols, TY_INT) + call salloc (npts, ncols, TY_INT) + call salloc (lines, nlines, TY_INT) + call salloc (wgt, ncols, TY_REAL) + call salloc (med, nbox * ncols, TY_REAL) + call salloc (z, ncols, TY_REAL) + call salloc (ranges, 3 * max_nranges + 1, TY_INT) + call amovkr (1., Memr[wgt], ncols) + + # Loop over median boxes in y. + lp = 0 + lineno = 0 + for (y1 = 1; y1 <= nimlines; y1 = y1 + YMEDIAN(imfit)) { + + lineno = lineno + 1 + current_line = y1 - 1 + y2 = min (y1 + YMEDIAN(imfit) - 1, nimlines) + + # If lines not in range, next image section. + lines_per_box = 0 + while (prl_nextlineno (gl, current_line) != EOF) { + if (current_line > y2) + break + lines_per_box = lines_per_box + 1 + } + if (lines_per_box < (YMEDIAN(imfit) * (MEDIAN_PERCENT(imfit)/100.))) + next + + # Loop over the image lines. + call aclri (Memi[npts], ncols) + do i = y1, y2 { + + # Get image line, and check the good regions list. + lbuf = imgl2r (im, i) + nranges = prl_get_ranges (gl, i, Memi[ranges], max_nranges) + if (nranges == 0) + next + nxpts = is_expand_ranges (Memi[ranges], Memi[colsfit], nimcols) + + # Loop over the median boxes in x. + cp= 0 + mbuf = med + for (x1 = 1; x1 <= nimcols; x1 = x1 + XMEDIAN(imfit)) { + x2 = min (x1 + XMEDIAN(imfit) - 1, nimcols) + ntemp = is_choose_rangesr (Memi[colsfit], Memr[lbuf], + Memr[mbuf+Memi[npts+cp]], nxpts, x1, x2) + Memi[npts+cp] = Memi[npts+cp] + ntemp + mbuf = mbuf + nbox + cp = cp + 1 + } + } + + # Calculate the medians. + nxpts = 0 + mbuf = med + do i = 1, ncols { + if (Memi[npts+i-1] > ((MEDIAN_PERCENT(imfit) / 100.) * nbox)) { + Memr[z+nxpts] = asokr (Memr[mbuf], Memi[npts+i-1], + (Memi[npts+i-1] + 1) / 2) + Memi[cols+nxpts] = i + nxpts = nxpts + 1 + } + mbuf = mbuf + nbox + } + + # Fit the line. + call islfit (sf, Memi[cols], lineno, Memr[z], Memr[wgt], nxpts, + SF_USER, ier) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call eprintf ("Warning: Too few columns to fit line: %d\n") + call pargi (lineno) + case SINGULAR: + call eprintf ("Warning: Solution singular for line: %d\n") + call pargi (lineno) + Memi[lines+lp] = lineno + lp = lp + 1 + default: + Memi[lines+lp] = lineno + lp = lp + 1 + } + } + + # Solve the surface. + call issolve (sf, Memi[lines], lp, ier) + + # Handle fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call error (0, "GOOD_MEDIANS: Cannot fit surface.") + case SINGULAR: + call eprintf ("Warning: Solution singular for surface.") + default: + # everyting OK + } + + # Free space. + call sfree (sp) +end + + +# SET_OUTIMAGE -- Procedure to write an output image of the desired type. + +procedure set_outimage (imin, imout, imfit, sf, rl) + +pointer imin # input image +pointer imout # output image +pointer imfit # pointer to the imsurfut header structure +pointer sf # pointer to the surface descriptor +pointer rl # pointer to the rejected pixel list regions list + +int i, k, ncols, nlines, max_nranges +long u[IM_MAXDIM], v[IM_MAXDIM] +real b1x, b2x, b1y, b2y +pointer sp, x, y, inbuf, outbuf, ranges + +int impnlr(), imgnlr() +real ims_divzero() +extern ims_divzero +errchk malloc, mfree, imgnlr, impnlr + +begin + ncols = IM_LEN(imin,1) + nlines = IM_LEN(imin,2) + max_nranges = ncols + + # Calculate transformation constants from real coordinates to + # median coordinates if median processing specified. + if (MEDIAN(imfit) == YES) { + b1x = (1. + XMEDIAN(imfit)) / (2. * XMEDIAN(imfit)) + b2x = (2. * ncols + XMEDIAN(imfit) - 1.) / (2. * XMEDIAN(imfit)) + b1y = (1. + YMEDIAN(imfit)) / (2. * YMEDIAN(imfit)) + b2y = (2. * nlines + YMEDIAN(imfit) - 1.) / (2. * YMEDIAN(imfit)) + } + + # Allocate space for x coordinates, initialize to image coordinates + # and transform to median coordinates. + call smark (sp) + call salloc (x, ncols, TY_REAL) + call salloc (y, ncols, TY_REAL) + call salloc (ranges, 3 * max_nranges + 1, TY_INT) + + # Intialize the x array. + do i = 1, ncols + Memr[x - 1 + i] = i + if (MEDIAN(imfit) == YES) + call amapr (Memr[x], Memr[x], ncols, 1., real (ncols), b1x, b2x) + + # loop over the images lines + call amovkl (long (1), v, IM_MAXDIM) + call amovkl (long (1), u, IM_MAXDIM) + do i = 1, nlines { + + # Get input and output image buffers. + if (TYPE_OUTPUT(imfit) != FIT) { + if (imgnlr (imin, inbuf, v) == EOF) + call error (0, "Error reading input image.") + } + if (impnlr (imout, outbuf, u) == EOF) + call error (0, "Error writing output image.") + + # Intialize y coordinates to image coordinates, and + # transform to median coordinates. + if (MEDIAN(imfit) == YES) { + Memr[y] = real (i) + call amapr (Memr[y], Memr[y], 1, 1., real (nlines), + b1y, b2y) + call amovkr (Memr[y], Memr[y+1], (ncols-1)) + } else + call amovkr (real (i), Memr[y], ncols) + + # Write output image. + switch (TYPE_OUTPUT(imfit)) { + case FIT: + call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols) + case CLEAN: + call clean_line (Memr[x], Memr[y], Memr[inbuf], ncols, nlines, + rl, sf, i, NGROW(imfit)) + call amovr (Memr[inbuf], Memr[outbuf], ncols) + case RESID: + call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols) + call asubr (Memr[inbuf], Memr[outbuf], Memr[outbuf], ncols) + case RESP: + call isvector (sf, Memr[x], Memr[y], Memr[outbuf], ncols) + if (IS_INDEF(DIV_MIN(imfit))) { + iferr (call adivr (Memr[inbuf], Memr[outbuf], Memr[outbuf], + ncols)) + call advzr (Memr[inbuf], Memr[outbuf], Memr[outbuf], + ncols, ims_divzero) + } else { + do k = 1, ncols { + if (Memr[outbuf-1+k] < DIV_MIN(imfit)) + Memr[outbuf-1+k] = 1. + else + Memr[outbuf-1+k] = Memr[inbuf-1+k] / + Memr[outbuf-1+k] + } + } + default: + call error (0, "SET_OUTIMAGE: Unknown output type.") + } + } + + # Free space + call sfree (sp) +end + + +# CLEAN_LINE -- Procedure to set weights of rejected points to zero + +procedure clean_line (x, y, z, ncols, nlines, rl, sf, line, ngrow) + +real x[ARB] # array of weights set to 1 +real y # y value of line +real z[ARB] # line of data +int ncols # number of image columns +int nlines # number of image lines +pointer rl # pointer to reject pixel list +pointer sf # surface fitting +int line # line number +int ngrow # radius for region growing + +int cp, j, k, nranges, dist, yreg_min, yreg_max, xreg_min, xreg_max +pointer sp, branges +real r2 +int prl_get_ranges(), is_next_number() +real iseval() + +begin + call smark (sp) + call salloc (branges, 3 * ncols + 1, TY_INT) + + r2 = ngrow ** 2 + yreg_min = max (1, line - ngrow) + yreg_max = min (nlines, line + ngrow) + + do j = yreg_min, yreg_max { + nranges = prl_get_ranges (rl, j, Memi[branges], ncols) + if (nranges == 0) + next + dist = int (sqrt (r2 - (j - line) ** 2)) + cp = 0 + while (is_next_number (Memi[branges], cp) != EOF) { + xreg_min = max (1, cp - dist) + xreg_max = min (ncols, cp + dist) + do k = xreg_min, xreg_max + z[k] = iseval (sf, x[k], y) + cp = xreg_max + } + } + + call sfree (sp) +end + + +# DO_REJECT -- Procedure to detect rejected pixels in an image. + +procedure do_reject (im, imfit, gl, sf, rl) + +pointer im # pointer to in put image +pointer imfit # pointer to image fitting structure +pointer gl # pointer to good regions list +pointer sf # pointer to surface descriptor +pointer rl # pointer to rejected pixel list + +int niter, nrejects +real sigma +int detect_rejects() +real get_sigma() +errchk prl_init, detect_rejects, get_sigma, refit_surface + +begin + # Initialize rejected pixel list. + call prl_init (rl, int(IM_LEN(im,1)), int(IM_LEN(im,2))) + + # Do an iterative rejection cycle on the image. + niter = 0 + repeat { + + # Get the sigma of the fit. + sigma = get_sigma (im, gl, sf, rl) + + # Detect rejected pixels. + nrejects = detect_rejects (im, imfit, gl, sf, rl, sigma) + + # If no rejected pixels quit, else refit surface. + if (nrejects == 0 || NITER(imfit) == 0) + break + call refit_surface (im, imfit, gl, sf, rl) + + niter = niter + 1 + + } until (niter == NITER(imfit)) +end + + +# REFIT_SURFACE -- Procedure tp refit the surface. + +procedure refit_surface (im, imfit, gl, sf, rl) + +pointer im # pointer to image +pointer imfit # pointer to surface fitting structure +pointer gl # pointer to good regions list +pointer sf # pointer to surface descriptor +pointer rl # pointer to rejected pixels list + +int i, ijunk, lp, ier, max_nranges +int ncols, nlines, npts, nfree, nrejects, nranges, ncoeff +pointer sp, cols, colsfit, lines, buf, fbuf, wgt, granges + +int prl_get_ranges(), grow_regions(), is_expand_ranges() +int is_choose_rangesr() +pointer imgl2r() +errchk smark, salloc, sfree, imgl2r +errchk iscoeff, islfit, issolve +errchk prl_get_ranges, grow_regions +errchk is_choose_rangesr + +begin + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + max_nranges = ncols + + # Allocate up temporary storage. + call smark (sp) + call salloc (cols, ncols, TY_INT) + call salloc (colsfit, ncols, TY_INT) + call salloc (lines, nlines, TY_INT) + call salloc (fbuf, ncols, TY_REAL) + call salloc (wgt, ncols, TY_REAL) + call salloc (granges, 3 * max_nranges + 1, TY_INT) + + # Initialize columns. + do i = 1, ncols + Memi[cols+i-1] = i + call amovi (Memi[cols], Memi[colsfit], ncols) + + # Get number of coefficients. + switch (SURFACE_TYPE(imfit)) { + case SF_LEGENDRE, SF_CHEBYSHEV: + ncoeff = XORDER(imfit) + case SF_SPLINE3: + ncoeff = XORDER(imfit) + 3 + case SF_SPLINE1: + ncoeff = XORDER(imfit) + 1 + } + + # Refit affected lines and solve for surface. + lp = 0 + do i = 1, nlines { + + # Determine whether image line is good. + if (gl != NULL) { + nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges) + if (nranges == 0) + next + } + + # Define rejected points with region growing. + call amovkr (1., Memr[wgt], ncols) + nrejects = grow_regions (Memr[wgt], ncols, nlines, rl, i, + NGROW(imfit)) + + # Get number of data points. + if (gl == NULL) + npts = ncols + else + npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols) + nfree = npts - nrejects + + # If no rejected pixels skip to next line. + if (nrejects == 0) { + if (nfree >= ncoeff ) { + Memi[lines+lp] = i + lp = lp + 1 + } + next + } + + # Read in image line. + buf = imgl2r (im, i) + if (buf == EOF) + call error (0, "REFIT_SURFACE: Error reading image.") + + # Select the data. + if (gl == NULL) { + npts = ncols + if (nfree >= ncoeff) + call islfit (sf, Memi[colsfit], i, Memr[buf], Memr[wgt], + npts, SF_USER, ier) + else + ier = NO_DEG_FREEDOM + } else { + ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], + Memr[fbuf], npts, 1, ncols) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt], Memr[wgt], + npts, 1, ncols) + if (nfree >= ncoeff) + call islfit (sf, Memi[colsfit], i, Memr[fbuf], Memr[wgt], + npts, SF_USER, ier) + else + ier = NO_DEG_FREEDOM + } + + # Evaluate fitting errors. + switch (ier) { + case NO_DEG_FREEDOM: + call eprintf ("REFIT_SURFACE: Too few points to fit line: %d\n") + call pargi (i) + case SINGULAR: + call eprintf ("REFIT_SURFACE: Solution singular for line: %d\n") + call pargi (i) + Memi[lines+lp] = i + lp = lp + 1 + default: + Memi[lines+lp] = i + lp = lp + 1 + } + } + + # Resolve surface. + call issolve (sf, Memi[lines], lp, ier) + + # Evaluate fitting errors for surface. + switch (ier) { + case NO_DEG_FREEDOM: + call error (0, "REFIT_SURFACE: Too few points to fit surface\n") + case SINGULAR: + call eprintf ("REFIT_SURFACE: Solution singular for surface\n") + default: + # everything OK + } + + call sfree (sp) + +end + + +# GROW_REGIONS -- Procedure to set weights of rejected points to zero. + +int procedure grow_regions (wgt, ncols, nlines, rl, line, ngrow) + +real wgt[ARB] # array of weights set to 1 +int ncols # number of image columnspoints +int nlines # number of images lines +pointer rl # pointer to reject pixel list +int line # line number +int ngrow # radius for region growing + +int cp, j, k, nrejects, nranges, max_nranges +int dist, yreg_min, yreg_max, xreg_min, xreg_max +pointer sp, branges +real r2 +int prl_get_ranges(), is_next_number() +errchk smark, salloc, sfree +errchk prl_get_ranges, is_next_number + +begin + max_nranges = ncols + + call smark (sp) + call salloc (branges, 3 * max_nranges + 1, TY_INT) + + r2 = ngrow ** 2 + nrejects = 0 + yreg_min = max (1, line - ngrow) + yreg_max = min (nlines, line + ngrow) + + do j = yreg_min, yreg_max { + nranges = prl_get_ranges (rl, j, Memi[branges], max_nranges) + if (nranges == 0) + next + dist = int (sqrt (r2 - (j - line) ** 2)) + cp = 0 + while (is_next_number (Memi[branges], cp) != EOF) { + xreg_min = max (1, cp - dist) + xreg_max = min (ncols, cp + dist) + do k = xreg_min, xreg_max { + if (wgt[k] > 0.) { + wgt[k] = 0. + nrejects = nrejects + 1 + } + } + cp = xreg_max + } + } + + call sfree (sp) + return (nrejects) +end + + +# GET_SIGMA -- Procedure to calculate the sigma of the surface fit + +real procedure get_sigma (im, gl, sf, rl) + +pointer im # pointer to image +pointer gl # pointer to good pixel list +pointer sf # pointer to surface deascriptor +pointer rl # pointer to rejected pixel list + +int i, ijunk, cp, nranges, npts, ntpts, ncols, nlines, max_nranges +pointer sp, colsfit, x, xfit, y, zfit, buf, fbuf, wgt, granges, branges +real sum, sigma +int prl_get_ranges(), is_next_number(), is_expand_ranges() +int is_choose_rangesr() +pointer imgl2r() +real asumr(), awssqr() +errchk smark, salloc, sfree, imgl2r + +begin + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + max_nranges = ncols + + # Allocate working space. + call smark (sp) + call salloc (colsfit, ncols, TY_INT) + call salloc (x, ncols, TY_REAL) + call salloc (xfit, ncols, TY_REAL) + call salloc (y, ncols, TY_REAL) + call salloc (fbuf, ncols, TY_REAL) + call salloc (zfit, ncols, TY_REAL) + call salloc (wgt, ncols, TY_REAL) + call salloc (granges, 3 * max_nranges + 1, TY_INT) + call salloc (branges, 3 * max_nranges + 1, TY_INT) + + # Intialize the x array. + do i = 1, ncols + Memr[x+i-1] = i + call amovr (Memr[x], Memr[xfit], ncols) + + sum = 0. + sigma = 0. + ntpts = 0 + + # Loop over the image. + do i = 1, nlines { + + # Check that line is in range. + if (gl != NULL) { + nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges) + if (nranges == 0) + next + npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols) + } + + # Read in image. + buf = imgl2r (im, i) + if (buf == EOF) + call error (0, "GET_SIGMA: Error reading image.") + + # Select appropriate data and fit. + call amovkr (real (i), Memr[y], ncols) + if (gl == NULL) { + npts = ncols + call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts) + call asubr (Memr[buf], Memr[zfit], Memr[zfit], npts) + } else { + ijunk = is_choose_rangesr (Memi[colsfit], Memr[x], Memr[xfit], + npts, 1, ncols) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf], + npts, 1, ncols) + call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts) + call asubr (Memr[fbuf], Memr[zfit], Memr[zfit], npts) + } + + # Get ranges of rejected pixels for the line and set weights. + call amovkr (1., Memr[wgt], ncols) + nranges = prl_get_ranges (rl, i, Memi[branges], max_nranges) + if (nranges > 0) { + cp = 0 + while (is_next_number (Memi[branges], cp) != EOF) + Memr[wgt+cp-1] = 0. + if (gl != NULL) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt], + Memr[wgt], npts, 1, ncols) + } + + # Calculate sigma. + sigma = sigma + awssqr (Memr[zfit], Memr[wgt], npts) + ntpts = ntpts + asumr (Memr[wgt], npts) + } + + call sfree (sp) + + return (sqrt (sigma / (ntpts - 1))) +end + + +# DETECT_REJECTS - Procedure to detect rejected pixels. + +int procedure detect_rejects (im, imfit, gl, sf, rl, sigma) + +pointer im # pointer to image +pointer imfit # pointer to surface fitting structure +pointer gl # pointer to good pixel list +pointer sf # pointer to surface descriptor +pointer rl # pointer to rejected pixel list +real sigma # standard deviation of fit + +int i, j, ijunk, cp, ncols, nlines, npts, nranges, nlrejects, ntrejects +int norejects, max_nranges +pointer sp, granges, branges, x, xfit, cols, colsfit, y, zfit, buf, fbuf +pointer wgt, list +real upper, lower + +int prl_get_ranges(), is_next_number(), is_make_ranges(), is_expand_ranges() +int is_choose_rangesr() +pointer imgl2r() + +begin + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + max_nranges = ncols + + # Allocate temporary space. + call smark (sp) + call salloc (x, ncols, TY_REAL) + call salloc (xfit, ncols, TY_REAL) + call salloc (cols, ncols, TY_INT) + call salloc (colsfit, ncols, TY_INT) + call salloc (y, ncols, TY_REAL) + call salloc (fbuf, ncols, TY_REAL) + call salloc (zfit, ncols, TY_REAL) + call salloc (wgt, ncols, TY_REAL) + call salloc (granges, 3 * max_nranges + 1, TY_INT) + call salloc (branges, 3 * max_nranges + 1, TY_INT) + call salloc (list, ncols, TY_INT) + + # Intialize x and column values. + do i = 1, ncols { + Memi[cols+i-1] = i + Memr[x+i-1] = i + } + call amovr (Memr[x], Memr[xfit], ncols) + call amovi (Memi[cols], Memi[colsfit], ncols) + + ntrejects = 0 + if (LOWER(imfit) <= 0.0) + lower = -MAX_REAL + else + lower = -sigma * LOWER(imfit) + if (UPPER(imfit) <= 0.0) + upper = MAX_REAL + else + upper = sigma * UPPER(imfit) + + # Loop over the image. + do i = 1, nlines { + + # Get ranges if appropriate. + if (gl != NULL) { + nranges = prl_get_ranges (gl, i, Memi[granges], max_nranges) + if (nranges == 0) + next + npts = is_expand_ranges (Memi[granges], Memi[colsfit], ncols) + } + + # Read in image. + buf = imgl2r (im, i) + if (buf == EOF) + call error (0, "GET_SIGMA: Error reading image.") + + # Select appropriate data and fit. + call amovkr (real (i), Memr[y], ncols) + if (gl == NULL) { + npts = ncols + call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts) + call asubr (Memr[buf], Memr[zfit], Memr[zfit], npts) + } else { + ijunk = is_choose_rangesr (Memi[colsfit], Memr[x], Memr[xfit], + npts, 1, ncols) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[buf], Memr[fbuf], + npts, 1, ncols) + call isvector (sf, Memr[xfit], Memr[y], Memr[zfit], npts) + call asubr (Memr[fbuf], Memr[zfit], Memr[zfit], npts) + } + + # Get ranges of rejected pixels for the line and set weights. + call amovkr (1., Memr[wgt], ncols) + nranges = prl_get_ranges (rl, i, Memi[branges], max_nranges) + norejects = 0 + if (nranges > 0) { + cp = 0 + while (is_next_number (Memi[branges], cp) != EOF) { + Memi[list+norejects] = cp + norejects = norejects + 1 + Memr[wgt+cp-1] = 0. + } + if (gl != NULL) + ijunk = is_choose_rangesr (Memi[colsfit], Memr[wgt], + Memr[wgt], npts, 1, ncols) + } + + # Detect deviant pixels. + nlrejects = 0 + do j = 1, npts { + if ((Memr[zfit+j-1] < lower || Memr[zfit+j-1] > upper) && + Memr[wgt+j-1] != 0.) { + Memi[list+norejects+nlrejects] = Memi[colsfit+j-1] + nlrejects = nlrejects + 1 + } + } + + # Add to rejected pixel list. + if (nlrejects > 0) { + call asrti (Memi[list], Memi[list], norejects + nlrejects) + nranges = is_make_ranges (Memi[list], norejects + nlrejects, + Memi[granges], max_nranges) + call prl_put_ranges (rl, i, i, Memi[granges]) + } + + ntrejects = ntrejects + nlrejects + } + + call sfree (sp) + return (ntrejects) +end + + +# AWSSQR -- Procedure to calculate the weighted sum of the squares + +real procedure awssqr (a, w, npts) + +real a[npts] # array of data +real w[npts] # array of points +int npts # number of data points + +int i +real sum + +begin + sum = 0. + do i = 1, npts + sum = sum + w[i] * a[i] ** 2 + + return (sum) +end + + +# IMS_DIVZER0 -- Return 1. on a divide by zero + +real procedure ims_divzero (x) + +real x + +begin + return (1.) +end diff --git a/pkg/images/imfit/src/mkpkg b/pkg/images/imfit/src/mkpkg new file mode 100644 index 00000000..3bc27b8f --- /dev/null +++ b/pkg/images/imfit/src/mkpkg @@ -0,0 +1,15 @@ +# Library for the IMAGES IMFIT Subpackage Tasks + +$checkout libpkg.a ../../ +$update libpkg.a +$checkin libpkg.a ../../ +$exit + +libpkg.a: + fit1d.x <imhdr.h> <pkg/gtools.h> <error.h> + imsurfit.x imsurfit.h <math/surfit.h> <imhdr.h> <mach.h> + pixlist.x pixlist.h + ranges.x <ctype.h> <mach.h> + t_imsurfit.x imsurfit.h <error.h> <imhdr.h> + t_lineclean.x <imhdr.h> <pkg/gtools.h> + ; diff --git a/pkg/images/imfit/src/pixlist.h b/pkg/images/imfit/src/pixlist.h new file mode 100644 index 00000000..39c875a9 --- /dev/null +++ b/pkg/images/imfit/src/pixlist.h @@ -0,0 +1,11 @@ +# PIXEL LIST descriptor structure + +define LEN_PLSTRUCT 10 + +define PRL_NCOLS Memi[$1] # number of columns +define PRL_NLINES Memi[$1+1] # number of lines +define PRL_LINES Memi[$1+2] # pointer to the line offsets +define PRL_LIST Memi[$1+3] # pointer to list of ranges +define PRL_SZLIST Memi[$1+4] # size of list in INTS +define PRL_LP Memi[$1+5] # offset to next space in list + diff --git a/pkg/images/imfit/src/pixlist.x b/pkg/images/imfit/src/pixlist.x new file mode 100644 index 00000000..066637fd --- /dev/null +++ b/pkg/images/imfit/src/pixlist.x @@ -0,0 +1,369 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "pixlist.h" + +.help pixels xtools "Pixel List Handling Tools" +.nf ________________________________________________________________________ +.fi +.ih +PURPOSE +These routines provide simple pixel list handling facilities and are +intended as a temporary facility pending full scale completion of +image masking. The list is stored in the form of ranges as a function +of line number. Each image line has a offset which may be NULL for +no entry or an offset into the list itself. The actual list is a set of +ranges with the ranges for each line delimited by a NULL. Routines +exist to fetch the ranges for a given line, add or append ranges to a +given line, fetch the next or previous line number with a non-NULL +range and specify whether two lines have the same ranges. At present +the list can grow indefinitely, with additional memory being added as +necessary. No attempt is made to clean up redundant entries though +such a faclity could easily be added. The ranges arguments conform +with the design of the ranges routinesr, with each range consisting +of and intitial and final entry and a step size. A list of ranges +is terminated with a NULL +.ih +PROCEDURE +.nf +prl_init (pl, ncols, nlines) + + pointer pl # pointer to list descriptor + int ncols # number of image columns + int nlines # number of image lines + +nranges = prl_get_ranges (pl, lineno, ranges, max_nranges) + + pointer pl # pointer to list descriptor + int lineno # line number of ranges to be fetched + int ranges[ARB] # ranges to be output + int max_nranges # the maximum number of ranges to be output + +prl_put_ranges (pl, linemin, linemax, ranges) + + pointer pl # pointer to list descriptor + int linemin # minimum line number + int linemax # maximum line number + int ranges[ARB] # ranges to be added to list + +prl_append_ranges (pl, linemin, linemax, ranges) + + pointer pl # pointer to list descriptor + int linemin # minimum line number + int linemax # maximum line number + int ranges[ARB] # ranges to be added to list + +next_lineno/EOF = prl_nextlineno (pl, current_lineno) + + pointer pl # pointer to list descriptor + int current_lineno # current line number + +prev_lineno/EOF = prl_prevlineno (pl, current_lineno) + + pointer pl # pointer to the list descriptor + int current_lineno # current line number + +YES/NO = prl_eqlines (pl, line1, line2) + + pointer pl # pointer to the list descriptor + int line1 # first line number + int line2 # second line number + +prl_free (pl) + + pointer pl # pointer to list descriptor +.fi +.endhelp ________________________________________________________________ + + +# PRL_ADD_RANGES -- Procedure to add the ranges for a given range of +# line numbers to the pixel list. The new ranges will be appended to any +# previously existing ranges for the specified line numbers. + +procedure prl_add_ranges (pl, linemin, linemax, ranges) + +pointer pl # pointer to the list descriptor +int linemin # minimum line number +int linemax # maximum line number +int ranges[ARB] # ranges + +int i, j, lc +int olp, lp, lnull, lold +int nr, nnewr, noldr + +begin + # check conditions + if ((linemin < 1) || (linemax > PRL_NLINES(pl)) || linemin > linemax) + return + + # calculate the length of the range to be appended minus the null + nr = 0 + while (ranges[nr+1] != NULL) + nr = nr + 1 + + + lc = 1 + olp = -1 + do i = linemin, linemax { + + # get offset for line i + lp = Memi[PRL_LINES(pl)+i-1] + + # if line pointer is undefined + if (lp == NULL) { + + if (lc == 1) { + + # set line pointer and store + Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl) + lnull = PRL_LP(pl) + + # check the size of the list + if (PRL_SZLIST(pl) < (nr + PRL_LP(pl))) { + PRL_SZLIST(pl) = PRL_SZLIST(pl) + nr + 1 + call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT) + } + + # move ranges and reset pointers + call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], nr) + PRL_LP(pl) = PRL_LP(pl) + nr + 1 + Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL + lc = lc + 1 + + } else + + # set line pointer + Memi[PRL_LINES(pl)+i-1] = lnull + + } else { + + if (lp != olp) { + + # set line pointer and store + Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl) + lold = PRL_LP(pl) + + # find length of previously defined range and calculate + # length of new ranges + for (j = lp; Memi[PRL_LIST(pl)+j-1] != NULL; j = j + 1) + ; + noldr = j - lp + nnewr = noldr + nr + + # check size of list + if (PRL_SZLIST(pl) < (nnewr + PRL_LP(pl))) { + PRL_SZLIST(pl) = PRL_SZLIST(pl) + nnewr + 1 + call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT) + } + + # add ranges to list and update pointers + call amovi (Memi[PRL_LIST(pl)+lp-1], + Memi[PRL_LIST(pl)+PRL_LP(pl)-1], noldr) + PRL_LP(pl) = PRL_LP(pl) + noldr + call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], nr) + PRL_LP(pl) = PRL_LP(pl) + nr + 1 + Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL + + } else + + # set line pointers + Memi[PRL_LINES(pl)+i-1] = lold + + olp = lp + } + } + +end + +# PRL_EQLINES -- Routine to test whether two lines have equal ranges. +# The routine returns YES or NO. + +int procedure prl_eqlines (pl, line1, line2) + +pointer pl # pointer to the list +int line1 # line numbers +int line2 + +begin + if (Memi[PRL_LINES(pl)+line1-1] == Memi[PRL_LINES(pl)+line2-1]) + return (YES) + else + return (NO) +end + +# PRL_GET_RANGES -- Procedure to fetch the ranges for the specified lineno. +# Zero is returned if there are no ranges otherwise the number of ranges +# are returned. The ranges are stored in an integer array. Three positive +# numbers are used to define a range a minimum, maximum and a step size. +# The ranges are delimited by a NULL. + +int procedure prl_get_ranges (pl, lineno, ranges, max_nranges) + +pointer pl # pointer to the pixel list descriptor +int lineno # line number +int ranges[ARB] # array of ranges +int max_nranges # the maximum number of ranges + +int lp, ip +int nranges + +begin + # check for existence of ranges + if (Memi[PRL_LINES(pl)+lineno-1] == NULL) { + ranges[1] = NULL + return (0) + } + + # set pointer to the first element in list for line lineno + lp = PRL_LIST(pl) + Memi[PRL_LINES(pl)+lineno-1] - 1 + + # get ranges + nranges = 0 + ip = 1 + while (Memi[lp+ip-1] != NULL && nranges <= 3 * max_nranges) { + ranges[ip] = Memi[lp+ip-1] + ip = ip + 1 + nranges = nranges + 1 + } + ranges[ip] = NULL + + # return nranges + if (nranges == 0) + return (nranges) + else + return (nranges / 3) +end + +# PRL_NEXTLINENO -- Procedure to fetch the next line number with a set of +# defined ranges given the current line number. Note that the current +# line number is updated. + +int procedure prl_nextlineno (pl, current_lineno) + +pointer pl # pointer to the pixel list descriptor +int current_lineno # current line number + +int findex, lp + +begin + findex = max (1, current_lineno + 1) + do lp = findex, PRL_NLINES(pl) { + if (Memi[PRL_LINES(pl)+lp-1] != NULL) { + current_lineno = lp + return (lp) + } + } + + return (EOF) +end + +# PRL_PREVLINENO -- Procedure to fetch the first previous line number +# with a set of defined ranges given the current line number. +# Note that the current line number is updated. + +int procedure prl_prevlineno (pl, current_lineno) + +pointer pl # pointer to the pixel list descriptor +int current_lineno # current line number + +int findex, lp + +begin + findex = min (current_lineno - 1, PRL_NLINES(pl)) + do lp = findex, 1, -1 { + if (Memi[PRL_LINES(pl)+lp-1] != NULL) { + current_lineno = lp + return (lp) + } + } + + return (EOF) +end + +# PRL_PUT_RANGES -- Procedure to add the ranges for a given range of +# lines to the pixel list. Note that any previously defined ranges are +# lost. + +procedure prl_put_ranges (pl, linemin, linemax, ranges) + +pointer pl # pointer to the list +int linemin # minimum line +int linemax # maximum line +int ranges[ARB] # list of ranges + +int i +int len_range + +begin + # check boundary conditions + if ((linemin < 1) || (linemax > PRL_NLINES(pl)) || (linemin > linemax)) + return + + # determine length of range string minus the NULL + len_range = 0 + while (ranges[len_range+1] != NULL) + len_range = len_range + 1 + + # check space allocation + if (PRL_SZLIST(pl) < (len_range + PRL_LP(pl))) { + PRL_SZLIST(pl) = PRL_SZLIST(pl) + len_range + 1 + call realloc (PRL_LIST(pl), PRL_SZLIST(pl), TY_INT) + } + + # set the line pointers + do i = linemin, linemax + Memi[PRL_LINES(pl)+i-1] = PRL_LP(pl) + + # add ranges + call amovi (ranges, Memi[PRL_LIST(pl)+PRL_LP(pl)-1], len_range) + PRL_LP(pl) = PRL_LP(pl) + len_range + 1 + Memi[PRL_LIST(pl)+PRL_LP(pl)-2] = NULL +end + + +# PLR_FREE -- Procedure to free the pixel list descriptor + +procedure prl_free (pl) + +pointer pl # pointer to pixel list descriptor + +begin + if (pl == NULL) + return + + if (PRL_LIST(pl) != NULL) + call mfree (PRL_LIST(pl), TY_INT) + if (PRL_LINES(pl) != NULL) + call mfree (PRL_LINES(pl), TY_INT) + + call mfree (pl, TY_STRUCT) +end + +# PRL_INIT -- Procedure to initialize the pixel list. Ncols and nlines are +# the number of columns and lines respectively in the associated IRAF +# image. + +procedure prl_init (pl, ncols, nlines) + +pointer pl # pixel list descriptor +int ncols # number of image columns +int nlines # number of image lines + +begin + # allocate space for a pixel list descriptor + call malloc (pl, LEN_PLSTRUCT, TY_STRUCT) + + # initialize + PRL_NCOLS(pl) = ncols + PRL_NLINES(pl) = nlines + + # allocate space for the line pointers + call malloc (PRL_LINES(pl), PRL_NLINES(pl), TY_INT) + call amovki (NULL, Memi[PRL_LINES(pl)], PRL_NLINES(pl)) + + # set pointer to next free element + PRL_LP(pl) = 1 + + # allocate space for the actual list + call malloc (PRL_LIST(pl), PRL_NLINES(pl), TY_INT) + PRL_SZLIST(pl) = PRL_NLINES(pl) +end diff --git a/pkg/images/imfit/src/ranges.x b/pkg/images/imfit/src/ranges.x new file mode 100644 index 00000000..19dc5c0e --- /dev/null +++ b/pkg/images/imfit/src/ranges.x @@ -0,0 +1,524 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <ctype.h> + +.help ranges xtools "Range Parsing Tools" +.ih +PURPOSE + +These tools +parse a string using a syntax to represent integer values, ranges, and +steps. The parsed string is used to generate a list of integers for various +purposes such as specifying lines or columns in an image or tape file numbers. +.ih +SYNTAX + +The syntax for the range string consists of positive integers, '-' (minus), +'x', ',' (comma), and whitespace. The commas and whitespace are ignored +and may be freely used for clarity. The remainder of the string consists +of sequences of five fields. The first field is the beginning of a range, +the second is a '-', the third is the end of the range, the fourth is +a 'x', and the fifth is a step size. Any of the five fields may be +missing causing various default actions. The defaults are illustrated in +the following table. + +.nf +-3x1 A missing starting value defaults to 1. +2-x1 A missing ending value defaults to MAX_INT. +2x1 A missing ending value defaults to MAX_INT. +2-4 A missing step defaults to 1. +4 A missing ending value and step defaults to an ending + value equal to the starting value and a step of 1. +x2 Missing starting and ending values defaults to + the range 1 to MAX_INT with the specified step. +"" The null string is equivalent to "1 - MAX_INT x 1", + i.e all positive integers. +.fi + +The specification of several ranges yields the union of the ranges. +.ih +EXAMPLES + +The following examples further illustrate the range syntax. + +.nf +- All positive integers. +1,5,9 A list of integers equivalent to 1-1x1,5-5x1,9-9x1. +x2 Every second positive integer starting with 1. +2x3 Every third positive integer starting with 2. +-10 All integers between 1 and 10. +5- All integers greater than or equal to 5. +9-3x1 The integers 3,6,9. +.fi +.ih +PROCEDURES + +.ls 4 is_decode_ranges + +.nf +int procedure is_decode_ranges (range_string, ranges, max_ranges, minimum, + maximum, nvalues) + +char range_string[ARB] # Range string to be decoded +int ranges[3, max_ranges] # Range array +int max_ranges # Maximum number of ranges +int minimum, maximum # Minimum and maximum range values allowed +int nvalues # The number of values in the ranges +.fi + +The range string is decoded into an integer array of maximum dimension +3 * max_ranges. Each range consists of three consecutive integers +corresponding to the starting and ending points of the range and the +step size. The number of integers covered by the ranges is returned +as nvalue. The end of the set of ranges is marked by a NULL. +The returned status is either ERR or OK. +.le +.ls 4 is_next_number, is_previous_number + +.nf +int procedure is_next_number (ranges, number) +int procedure is_previous_number (ranges, number) + +int ranges[ARB] # Range array +int number # Both input and output parameter +.fi + +Given a value for number the procedures find the next (previous) number in +increasing (decreasing) +value within the set of ranges. The next (previous) number is returned in +the number argument. A returned status is either OK or EOF. +EOF indicates that there are no greater values. The usual usage would +be in a loop of the form: + +.nf + number = 0 + while (is_next_number (ranges, number) != EOF) { + <Statements using number> + } +.fi +.le +.ls 4 is_in_rangelist + +.nf +bool procedure is_in_rangelist (ranges, number) + +int ranges[ARB] # Ranges array +int number # Number to check againts ranges +.fi + +A boolean value is returned indicating whether number is covered by +the ranges. + +.endhelp + + +# IS_DECODE_RANGES -- Parse a string containing a list of integer numbers or +# ranges, delimited by either spaces or commas. Return as output a list +# of ranges defining a list of numbers, and the count of list numbers. +# Range limits must be positive nonnegative integers. ERR is returned as +# the function value if a conversion error occurs. The list of ranges is +# delimited by a single NULL. + + +int procedure is_decode_ranges (range_string, ranges, max_ranges, minimum, + maximum, nvalues) + +char range_string[ARB] # Range string to be decoded +int ranges[3, max_ranges] # Range array +int max_ranges # Maximum number of ranges +int minimum, maximum # Minimum and maximum range values allowed +int nvalues # The number of values in the ranges + +int ip, nrange, out_of_range, a, b, first, last, step, ctoi() + +begin + ip = 1 + nrange = 1 + nvalues = 0 + out_of_range = 0 + + while (nrange < max_ranges) { + # Default values + a = minimum + b = maximum + step = 1 + + # Skip delimiters + while (IS_WHITE(range_string[ip]) || range_string[ip] == ',') + ip = ip + 1 + + # Get first limit. + # Must be a number, '*', '-', 'x', or EOS. If not return ERR. + if (range_string[ip] == EOS) { # end of list + if (nrange == 1) { + if (out_of_range == 0) { + # Null string defaults + ranges[1, 1] = a + ranges[2, 1] = b + ranges[3, 1] = step + ranges[1, 2] = NULL + nvalues = (b - a) / step + 1 + return (OK) + } else { + # Only out of range data + return (ERR) + } + } else { + ranges[1, nrange] = NULL + return (OK) + } + } else if (range_string[ip] == '-') + ; + else if (range_string[ip] == '*') + ; + else if (range_string[ip] == 'x') + ; + else if (IS_DIGIT(range_string[ip])) { # ,n.. + if (ctoi (range_string, ip, a) == 0) + return (ERR) + } else + return (ERR) + + # Skip delimiters + while (IS_WHITE(range_string[ip]) || range_string[ip] == ',') + ip = ip + 1 + + # Get last limit + # Must be '-', '*', or 'x' otherwise b = a. + if (range_string[ip] == 'x') + ; + else if ((range_string[ip] == '-') || (range_string[ip] == '*')) { + ip = ip + 1 + while (IS_WHITE(range_string[ip]) || range_string[ip] == ',') + ip = ip + 1 + if (range_string[ip] == EOS) + ; + else if (IS_DIGIT(range_string[ip])) { + if (ctoi (range_string, ip, b) == 0) + return (ERR) + } else if (range_string[ip] == 'x') + ; + else + return (ERR) + } else + b = a + + # Skip delimiters + while (IS_WHITE(range_string[ip]) || range_string[ip] == ',') + ip = ip + 1 + + # Get step. + # Must be 'x' or assume default step. + if (range_string[ip] == 'x') { + ip = ip + 1 + while (IS_WHITE(range_string[ip]) || range_string[ip] == ',') + ip = ip + 1 + if (range_string[ip] == EOS) + ; + else if (IS_DIGIT(range_string[ip])) { + if (ctoi (range_string, ip, step) == 0) + ; + } else if (range_string[ip] == '-') + ; + else if (range_string[ip] == '*') + ; + else + return (ERR) + } + + # Output the range triple. + first = min (a, b) + last = max (a, b) + if (first < minimum) + first = minimum + mod (step - mod (minimum - first, step), step) + if (last > maximum) + last = maximum - mod (last - maximum, step) + if (first <= last) { + ranges[1, nrange] = first + ranges[2, nrange] = last + ranges[3, nrange] = step + nvalues = nvalues + (last - first) / step + 1 + nrange = nrange + 1 + } else + out_of_range = out_of_range + 1 + } + + return (ERR) # ran out of space +end + + +# IS_NEXT_NUMBER -- Given a list of ranges and the current file number, +# find and return the next file number. Selection is done in such a way +# that list numbers are always returned in monotonically increasing order, +# regardless of the order in which the ranges are given. Duplicate entries +# are ignored. EOF is returned at the end of the list. + +int procedure is_next_number (ranges, number) + +int ranges[ARB] # Range array +int number # Both input and output parameter + +int ip, first, last, step, next_number, remainder + +begin + # If number+1 is anywhere in the list, that is the next number, + # otherwise the next number is the smallest number in the list which + # is greater than number+1. + + number = number + 1 + next_number = MAX_INT + + for (ip=1; ranges[ip] != NULL; ip=ip+3) { + first = ranges[ip] + last = ranges[ip+1] + step = ranges[ip+2] + if (number >= first && number <= last) { + remainder = mod (number - first, step) + if (remainder == 0) + return (number) + if (number - remainder + step <= last) + next_number = number - remainder + step + } else if (first > number) + next_number = min (next_number, first) + } + + if (next_number == MAX_INT) + return (EOF) + else { + number = next_number + return (number) + } +end + + +# IS_PREVIOUS_NUMBER -- Given a list of ranges and the current file number, +# find and return the previous file number. Selection is done in such a way +# that list numbers are always returned in monotonically decreasing order, +# regardless of the order in which the ranges are given. Duplicate entries +# are ignored. EOF is returned at the end of the list. + +int procedure is_previous_number (ranges, number) + +int ranges[ARB] # Range array +int number # Both input and output parameter + +int ip, first, last, step, next_number, remainder + +begin + # If number-1 is anywhere in the list, that is the previous number, + # otherwise the previous number is the largest number in the list which + # is less than number-1. + + number = number - 1 + next_number = 0 + + for (ip=1; ranges[ip] != NULL; ip=ip+3) { + first = ranges[ip] + last = ranges[ip+1] + step = ranges[ip+2] + if (number >= first && number <= last) { + remainder = mod (number - first, step) + if (remainder == 0) + return (number) + if (number - remainder >= first) + next_number = number - remainder + } else if (last < number) { + remainder = mod (last - first, step) + if (remainder == 0) + next_number = max (next_number, last) + else if (last - remainder >= first) + next_number = max (next_number, last - remainder) + } + } + + if (next_number == 0) + return (EOF) + else { + number = next_number + return (number) + } +end + + +# IS_IN_RANGELLIST -- Test number to see if it is in range. + +bool procedure is_in_rangelist (ranges, number) + +int ranges[ARB] # Range array +int number # Number to be tested against ranges + +int ip, first, last, step + +begin + for (ip=1; ranges[ip] != NULL; ip=ip+3) { + first = ranges[ip] + last = ranges[ip+1] + step = ranges[ip+2] + if (number >= first && number <= last) + if (mod (number - first, step) == 0) + return (TRUE) + } + + return (FALSE) +end + + +# IS_EXPAND_RANGES -- Expand a range string into a array of values. + +int procedure is_expand_ranges (ranges, array, max_nvalues) + +int ranges[ARB] # Range array +int array[max_nvalues] # Array of values +int max_nvalues # Maximum number of values + +int n, value + +int is_next_number() + +begin + n = 0 + value = 0 + while ((n < max_nvalues) && (is_next_number (ranges, value) != EOF)) { + n = n + 1 + array[n] = value + } + + return (n) +end + + +# IS_SELECT_RANGES -- Select array values in the ranges. +# The input and output arrays may be the same. + +procedure is_select_ranges (a, b, ranges) + +real a[ARB] # Input array +real b[ARB] # Output array +int ranges[3, ARB] # Ranges + +int i, j, npts, nmove + +begin + npts = 0 + for (i = 1; ranges[1, i] != NULL; i = i + 1) { + if (ranges[3, i] == 1) { + nmove = ranges[2, i] - ranges[1, i] + 1 + call amovr (a[ranges[1, i]], b[npts + 1], nmove) + npts = npts + nmove + } else { + do j = ranges[1, i], ranges[2, i], ranges[3, i] { + npts = npts + 1 + b[npts] = a[j] + } + } + } +end + + +# IS_CHOOSE_RANGESI -- Copy the selected values from array a to b. + +int procedure is_choose_rangesi (indices, a, b, npts, ifirst, ilast) + +int indices[ARB] # array of indices +int a[ARB] # input array +int b[ARB] # output array +int npts # number of points +int ifirst # first index +int ilast # last index + +int i, element + +begin + element = 1 + do i = 1, npts { + if (indices[i] < ifirst || indices[i] > ilast) + next + b[element] = a[indices[i]] + element = element + 1 + } + return (element - 1) +end + + +# IS_CHOOSE_RANGESR -- Copy the selected values from array a to b. + +int procedure is_choose_rangesr (indices, a, b, npts, ifirst, ilast) + +int indices[ARB] # array of indices +real a[ARB] # input array +real b[ARB] # output array +int npts # number of points +int ifirst # first element to be extracted +int ilast # last element to be extracted + +int i, element + +begin + element = 1 + do i = 1, npts { + if (indices[i] < ifirst || indices[i] > ilast) + next + b[element] = a[indices[i]] + element = element + 1 + } + return (element - 1) +end + + +# IS_MAKE_RANGES -- Procedure to make a set of ranges from an ordered list +# of column numbers. Only a step size of 1 is checked for. + +int procedure is_make_ranges (list, npts, ranges, max_nranges) + +int list[ARB] # list of column numbers in increasing order +int npts # number of list elements +int ranges[ARB] # output ranges +int max_nranges # the maximum number of ranges + +bool next_range +int ip, op, nranges + +begin + # If zero list elements return + if (npts == 0) { + ranges[1] = NULL + return (0) + } + + # Initialize + nranges = 0 + ranges[1] = list[1] + op = 2 + next_range = false + + # Loop over column list + for (ip = 2; ip <= npts && nranges < max_nranges; ip = ip + 1) { + if ((list[ip] != (list[ip-1] + 1))) { + ranges[op] = list[ip-1] + op = op + 1 + ranges[op] = 1 + op = op + 1 + nranges = nranges + 1 + ranges[op] = list[ip] + op = op + 1 + } + } + + # finish off + if (npts == 1) { + ranges[op] = list[npts] + ranges[op+1] = 1 + ranges[op+2] = NULL + nranges = 1 + } else if (nranges == max_nranges) { + ranges[op-1] = NULL + } else { + ranges[op] = list[npts] + ranges[op+1] = 1 + ranges[op+2] = NULL + nranges = nranges + 1 + } + + return (nranges) +end diff --git a/pkg/images/imfit/src/t_imsurfit.x b/pkg/images/imfit/src/t_imsurfit.x new file mode 100644 index 00000000..2a93b8b2 --- /dev/null +++ b/pkg/images/imfit/src/t_imsurfit.x @@ -0,0 +1,400 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imsurfit.h" + +# T_IMSURFIT -- Fit a surface function to an image +# +# 1. A user selected function is fit to each surface. +# 2. Only the selected regions of the image are fit. +# 3. Deviant pixels may be rejected from the fit. +# 4. The user selects the type of output image. The choices are: +# a. the fitted image. +# b. the input image with deviant pixels replaced by +# the fitted values +# c. the input image minus the fitted image. +# d. the ratio of the input image and the fit where +# pixels less than div_min are set to a ratio of 1. + + +procedure t_imsurfit () + +char imtlist1[SZ_LINE] # Input image list +char imtlist2[SZ_LINE] # Output image list +char image1[SZ_FNAME] # Input image +char image2[SZ_FNAME] # Output image + +char str[SZ_LINE], region_str[SZ_LINE] +int list1, list2, region_type +pointer im1, im2, imfit, gl, sp + +bool clgetb() +int imtopen(), imtgetim(), imtlen(), btoi(), clgeti(), clgwrd() +pointer immap() +real clgetr() + +begin + # Allocate space for imfit structure. + call smark (sp) + call salloc (imfit, LEN_IMSFSTRUCT, TY_STRUCT) + + # Get task parameters. + call clgstr ("input", imtlist1, SZ_FNAME) + call clgstr ("output", imtlist2, SZ_FNAME) + TYPE_OUTPUT(imfit) = clgwrd ("type_output", str, SZ_LINE, + ",fit,clean,residual,response,") + DIV_MIN(imfit) = clgetr ("div_min") + + # Get surface ftting parameters. + SURFACE_TYPE(imfit) = clgwrd ("function", str, SZ_LINE, + ",legendre,chebyshev,spline3,spline1,") + XORDER(imfit) = clgeti ("xorder") + YORDER(imfit) = clgeti ("yorder") + CROSS_TERMS(imfit) = btoi (clgetb ("cross_terms")) + + # Get median processing parameters. + XMEDIAN(imfit) = clgeti ("xmedian") + YMEDIAN(imfit) = clgeti ("ymedian") + MEDIAN_PERCENT(imfit) = clgetr ("median_percent") + if (XMEDIAN(imfit) > 1 || YMEDIAN(imfit) > 1) + MEDIAN(imfit) = YES + else + MEDIAN(imfit) = NO + + # Get rejection cycle parameters. + NITER(imfit) = clgeti ("niter") + LOWER(imfit) = clgetr ("lower") + UPPER(imfit) = clgetr ("upper") + NGROW(imfit) = clgeti ("ngrow") + + if (MEDIAN(IMFIT) == YES) { + REJECT(imfit) = NO + NITER(imfit) = 0 + } else if (NITER(imfit) > 0 && (LOWER(imfit) > 0. || UPPER(imfit) > 0.)) + REJECT(imfit) = YES + else { + REJECT(imfit) = NO + NITER(imfit) = 0 + } + + # Checking sigmas for cleaning. + if (TYPE_OUTPUT(imfit) == CLEAN && MEDIAN(imfit) == YES) + call error (0, + "T_IMSURFIT: Clean option and median processing are exclusive.") + if (TYPE_OUTPUT(imfit) == CLEAN && NITER(imfit) <= 0) + call error (0, "T_IMSURFIT: Clean option requires non-zero niter.") + if (TYPE_OUTPUT(imfit) == CLEAN && LOWER(imfit) <= 0. && + UPPER(imfit) <= 0.) + call error (0, "T_IMSURFIT: Clean option requires non-zero sigma.") + + # Get regions to be fit. + gl = NULL + region_type = clgwrd ("regions", str, SZ_LINE, + ",all,columns,rows,border,sections,circle,invcircle,") + switch (region_type) { + case ALL: + ; + case BORDER: + call clgstr ("border", region_str, SZ_LINE) + case SECTIONS: + call clgstr ("sections", region_str, SZ_LINE) + case COLUMNS: + call clgstr ("columns", region_str, SZ_LINE) + case ROWS: + call clgstr ("rows", region_str, SZ_LINE) + case CIRCLE, INVCIRCLE: + call clgstr ("circle", region_str, SZ_LINE) + } + + # Expand the input and output image lists. + list1 = imtopen (imtlist1) + list2 = imtopen (imtlist2) + if (imtlen (list1) != imtlen (list2)) { + call imtclose (list1) + call imtclose (list2) + call error (0, "Number of input and output images not the same.") + } + + # Do each set of input and output images. + while ((imtgetim (list1, image1, SZ_FNAME) != EOF) && + (imtgetim (list2, image2, SZ_FNAME) != EOF)) { + + im1 = immap (image1, READ_ONLY, 0) + im2 = immap (image2, NEW_COPY, im1) + + iferr { + if (region_type != ALL) + call make_good_list (im1, gl, region_type, region_str) + call imsurfit (im1, im2, imfit, gl) + } then + call erract (EA_WARN) + + call imunmap (im1) + call imunmap (im2) + call prl_free (gl) + } + + # Cleanup. + call sfree (sp) + call imtclose (list1) + call imtclose (list2) +end + + +# MAKE_GOOD_LIST -- Procedure to make a list of good regions. The program +# returns an error message if no good regions are defined. The good +# list parameter is set to NULL if the whole image is to be fit. This routine +# uses both the ranges and pixlist package which will be replaced by image +# masking. + +procedure make_good_list (im, gl, region_type, region_string) + +pointer im # pointer to the image +pointer gl # good pixel list descriptor +int region_type # type of good region +char region_string[ARB] # region parameters + +int i, ip, zero, nvals, range_min, r2, xdist, max_nranges +int x1, x2, y1, y2, temp, border, xcenter, ycenter, radius +int columns[7] +pointer sp, ranges, list + +bool is_in_rangelist() +int is_next_number(), is_decode_ranges(), open(), fscan(), nscan(), ctoi() +errchk open, close + +begin + # Determine the maximum number of images. + max_nranges = IM_LEN(im,1) + + # Allocate working space. + call smark (sp) + call salloc (ranges, 3 * max_nranges + 1, TY_INT) + + # Compute the good pixel list. + switch (region_type) { + case ROWS: + + # Decode the row ranges. + if (is_decode_ranges (region_string, Memi[ranges], max_nranges, 1, + int (IM_LEN(im,2)), nvals) == ERR) + call error (0, "MAKE_GOOD_LIST: Error decoding row string.") + if (nvals == 0) + call error (0, "MAKE_GOOD_LIST: no good rows.") + if (nvals == IM_LEN(im,2)) { + call sfree (sp) + return + } + + # Intialize the good pixel list. + call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2))) + + # Set column limits using the ranges format. + columns[1] = 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + + # Set column limits for the specied lines. + zero = 0 + range_min = is_next_number (Memi[ranges], zero) + while (range_min != EOF) { + for (i = range_min; i <= IM_LEN(im,2) + 1; i = i + 1) { + if (!is_in_rangelist (Memi[ranges], i) || + i == IM_LEN(im,2)+1) { + call prl_put_ranges (gl, range_min, i-1, columns) + break + } + } + range_min = is_next_number (Memi[ranges], i) + } + + case COLUMNS: + + # Set the specified columns. + if (is_decode_ranges (region_string, Memi[ranges], max_nranges, 1, + int (IM_LEN(im,1)), nvals) == ERR) + call error (0, "MAKE_GOOD_LIST: Error decoding column string.") + if (nvals == 0) + call error (0, "MAKE_GOOD_LIST: No good columns.") + if (nvals == IM_LEN(im,1)) { + call sfree (sp) + return + } + + # Make the good pixel list. + call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2))) + call prl_add_ranges (gl, 1, int (IM_LEN(im,2)), Memi[ranges]) + + case CIRCLE, INVCIRCLE: + + # Get the parameters of the circle. + ip = 1 + if (ctoi (region_string, ip, xcenter) <= 0) + call error (0, "MAKE_GOOD_LIST: Error decoding xcenter.") + if (ctoi (region_string, ip, ycenter) <= 0) + call error (0, "MAKE_GOOD_LIST: Error decoding ycenter.") + if (ctoi (region_string, ip, radius) <= 0) + call error (0, "MAKE_GOOD_LIST: Error decoding radius.") + + y1 = max (1, ycenter - radius) + y2 = min (int (IM_LEN(im,2)), ycenter + radius) + x1 = max (1, xcenter - radius) + x2 = min (int (IM_LEN(im,1)), xcenter + radius) + if (region_type == CIRCLE) { + if (y1 > IM_LEN(im,2) || y2 < 1 || x1 > IM_LEN(im,1) || x2 < 1) + call error (0, "MAKE_GOOD_LIST: No good regions.") + } + + # Create the good pixel list. + call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2))) + + r2 = radius ** 2 + if (region_type == CIRCLE) { + do i = y1, y2 { + xdist = sqrt (real (r2 - (ycenter - i) ** 2)) + x1 = max (1, xcenter - xdist) + x2 = min (IM_LEN(im,1), xcenter + xdist) + columns[1] = x1 + columns[2] = x2 + columns[3] = 1 + columns[4] = NULL + call prl_put_ranges (gl, i, i, columns) + } + } else if (region_type == INVCIRCLE) { + do i = 1, y1 - 1 { + columns[1] = 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + call prl_put_ranges (gl, i, i, columns) + } + do i = y2 + 1, IM_LEN(im,2) { + columns[1] = 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + call prl_put_ranges (gl, i, i, columns) + } + do i = y1, y2 { + xdist = sqrt (real (r2 - (ycenter - i) ** 2)) + x1 = max (1, xcenter - xdist) + x2 = min (IM_LEN(im,1), xcenter + xdist) + if (x1 > 1) { + columns[1] = 1 + columns[2] = x1 - 1 + columns[3] = 1 + if (x2 < IM_LEN(im,1)) { + columns[4] = x2 + 1 + columns[5] = IM_LEN(im,1) + columns[6] = 1 + columns[7] = NULL + } else + columns[4] = NULL + } else if (x2 < IM_LEN(im,1)) { + columns[1] = x2 + 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + } else + columns[1] = NULL + call prl_put_ranges (gl, i, i, columns) + } + } + + + case SECTIONS: + + # Open file of sections. + list = open (region_string, READ_ONLY, TEXT_FILE) + call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2))) + + # Scan the list. + while (fscan (list) != EOF) { + + # Fetch parameters from list. + call gargi (x1) + call gargi (x2) + call gargi (y1) + call gargi (y2) + if (nscan() != 4) + next + + # Check and correct for out of bounds limits. + x1 = max (1, min (IM_LEN(im,1), x1)) + x2 = min (IM_LEN(im,1), max (1, x2)) + y1 = max (1, min (IM_LEN(im,2), y1)) + y2 = min (IM_LEN(im,2), max (1, y2)) + + # Check the order. + if (x2 < x1) { + temp = x1 + x1 = x2 + x2 = temp + } + if (y2 < y1) { + temp = y1 + y1 = y2 + y2 = temp + } + + # If entire image return. + if ((x1 == 1) && (x2 == IM_LEN(im,1)) && (y1 == 1) && + (y2 == IM_LEN(im,2))) { + call prl_free (gl) + gl = NULL + break + } + + # Set ranges. + columns[1] = x1 + columns[2] = x2 + columns[3] = 1 + columns[4] = NULL + call prl_add_ranges (gl, y1, y2, columns) + } + + call close (list) + + case BORDER: + + # Decode border parameter. + ip = 1 + if (ctoi (region_string, ip, border) == ERR) + call error (0, "MAKE_GOOD_LIST: Error decoding border string.") + if (border < 1) + call error (0, "MAKE_GOOD_LIST: No border.") + if ((border > IM_LEN(im,1)/2) && (border > IM_LEN(im,2)/2)) { + call sfree (sp) + return + } + + # Intialize list. + call prl_init (gl, int (IM_LEN(im,1)), int (IM_LEN(im,2))) + y1 = 1 + border - 1 + y2 = IM_LEN(im,2) - border + 1 + columns[1] = 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + + # Set ranges for top and bottom edges of image. + call prl_put_ranges (gl, 1, y1, columns) + call prl_put_ranges (gl, y2, int (IM_LEN(im,2)), columns) + + columns[1] = 1 + columns[2] = y1 + columns[3] = 1 + columns[4] = NULL + call prl_put_ranges (gl, y1 + 1, y2 - 1, columns) + + columns[1] = IM_LEN(im,1) - border + 1 + columns[2] = IM_LEN(im,1) + columns[3] = 1 + columns[4] = NULL + call prl_add_ranges (gl, y1 + 1, y2 - 1, columns) + } + + call sfree (sp) +end diff --git a/pkg/images/imfit/src/t_lineclean.x b/pkg/images/imfit/src/t_lineclean.x new file mode 100644 index 00000000..4acb9752 --- /dev/null +++ b/pkg/images/imfit/src/t_lineclean.x @@ -0,0 +1,270 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <pkg/gtools.h> + +# LINECLEAN -- Fit a function to the image lines and output an image +# with rejected points replaced by the fit. The fitting parameters may be +# set interactively using the icfit package. + +procedure t_lineclean () + +int listin # Input image list +int listout # Output image list +char sample[SZ_LINE] # Sample ranges +int naverage # Sample averaging size +char function[SZ_LINE] # Curve fitting function +int order # Order of curve fitting function +real low_reject, high_reject # Rejection threshold +int niterate # Number of rejection iterations +real grow # Rejection growing radius +bool interactive # Interactive? + +char input[SZ_LINE] # Input image +char output[SZ_FNAME] # Output image +pointer in, out # IMIO pointers +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer + +int imtopen(), imtgetim(), imtlen(), gt_init() +int clgeti() +real clgetr() +bool clgetb() + +begin + # Get input and output lists and check that the number of images + # are the same. + + call clgstr ("input", input, SZ_LINE) + listin = imtopen (input) + call clgstr ("output", input, SZ_LINE) + listout = imtopen (input) + if (imtlen (listin) != imtlen (listout)) { + call imtclose (listin) + call imtclose (listout) + call error (0, "Input and output image lists do not match") + } + + # Get task parameters. + + call clgstr ("sample", sample, SZ_LINE) + naverage = clgeti ("naverage") + call clgstr ("function", function, SZ_LINE) + order = clgeti ("order") + low_reject = clgetr ("low_reject") + high_reject = clgetr ("high_reject") + niterate = clgeti ("niterate") + grow = clgetr ("grow") + interactive = clgetb ("interactive") + + + # Set the ICFIT pointer structure. + call ic_open (ic) + call ic_pstr (ic, "sample", sample) + call ic_puti (ic, "naverage", naverage) + call ic_pstr (ic, "function", function) + call ic_puti (ic, "order", order) + call ic_putr (ic, "low", low_reject) + call ic_putr (ic, "high", high_reject) + call ic_puti (ic, "niterate", niterate) + call ic_putr (ic, "grow", grow) + call ic_pstr (ic, "ylabel", "") + + gt = gt_init() + call gt_sets (gt, GTTYPE, "line") + + # Clean the lines in each input image. + + while ((imtgetim (listin, input, SZ_FNAME) != EOF) && + (imtgetim (listout, output, SZ_FNAME) != EOF)) { + + call lc_immap (input, output, in, out) + call lineclean (in, out, ic, gt, input, interactive) + call imunmap (in) + call imunmap (out) + } + + call ic_closer (ic) + call gt_free (gt) + call imtclose (listin) + call imtclose (listout) +end + + +# LINECLEAN -- Given the image descriptor determine the fitting function +# for each line and create an output image. If the interactive flag +# is set then set the fitting parameters interactively. + +procedure lineclean (in, out, ic, gt, title, interactive) + +pointer in # IMIO pointer for input image +pointer out # IMIO pointer for output image +pointer ic # ICFIT pointer +pointer gt # GTIO pointer +char title[ARB] # Title +bool interactive # Interactive? + +char graphics[SZ_FNAME] +int i, nx, new +long inline[IM_MAXDIM], outline[IM_MAXDIM] +pointer cv, gp, sp, x, wts, indata, outdata + +int lf_getline(), imgnlr(), impnlr(), strlen() +pointer gopen() + +begin + # Allocate memory for curve fitting. + + nx = IM_LEN(in, 1) + + call smark (sp) + call salloc (x, nx, TY_REAL) + call salloc (wts, nx, TY_REAL) + + do i = 1, nx + Memr[x+i-1] = i + call amovkr (1., Memr[wts], nx) + + call ic_putr (ic, "xmin", Memr[x]) + call ic_putr (ic, "xmax", Memr[x+nx-1]) + + # If the interactive flag is set then use icg_fit to set the + # fitting parameters. Get_fitline returns EOF when the user + # is done. The weights are reset since the user may delete + # points. + + if (interactive) { + call clgstr ("graphics", graphics, SZ_FNAME) + gp = gopen ("stdgraph", NEW_FILE, STDGRAPH) + + i = strlen (title) + while (lf_getline (ic, gt, in, indata, inline, title)!=EOF) { + title[i + 1] = EOS + call icg_fit (ic, gp, "cursor", gt, cv, Memr[x], Memr[indata], + Memr[wts], nx) + call amovkr (1., Memr[wts], nx) + } + call gclose (gp) + } + + # Loop through each input image line and create an output image line. + + new = YES + call amovkl (long(1), inline, IM_MAXDIM) + call amovkl (long(1), outline, IM_MAXDIM) + + while (imgnlr (in, indata, inline) != EOF) { + if (impnlr (out, outdata, outline) == EOF) + call error (0, "Error writing output image") + + call ic_fit (ic, cv, Memr[x], Memr[indata], Memr[wts], + nx, new, YES, new, new) + new = NO + + call ic_clean (ic, cv, Memr[x], Memr[indata], Memr[wts], nx) + + call amovr (Memr[indata], Memr[outdata], nx) + } + + call cvfree (cv) + call sfree (sp) +end + + +# LC_IMMAP -- Map images for lineclean. + +procedure lc_immap (input, output, in, out) + +char input[ARB] # Input image +char output[ARB] # Output image +pointer in # Input IMIO pointer +pointer out # Output IMIO pointer + +pointer sp, root, sect +int imaccess() +pointer immap() + +begin + # Get the root name and section of the input image. + + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (sect, SZ_FNAME, TY_CHAR) + + call get_root (input, Memc[root], SZ_FNAME) + call get_section (input, Memc[sect], SZ_FNAME) + + # If the output image is not accessible then create it as a new copy + # of the full input image. + + if (imaccess (output, 0) == NO) + call img_imcopy (Memc[root], output, false) + + # Map the input and output images. + + in = immap (input, READ_ONLY, 0) + + call sprintf (Memc[root], SZ_FNAME, "%s%s") + call pargstr (output) + call pargstr (Memc[sect]) + out = immap (Memc[root], READ_WRITE, 0) + + call sfree (sp) +end + + +# LF_GETLINE -- Get an image line to be fit interactively. Return EOF +# when the user enters EOF or CR. Default line is the first line and +# the out of bounds lines are silently limited to the nearest in bounds line. + +int procedure lf_getline (ic, gt, im, data, v, title) + +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer +pointer im # IMIO pointer +pointer data # Image data +long v[ARB] # Image line vector +char title[ARB] # Title + +int i +char line[SZ_LINE] +int getline(), nscan(), imgnlr() + +begin + call sprintf (title, SZ_LINE, "%s: Fit line =") + call pargstr (title) + + call printf ("%s ") + call pargstr (title) + call flush (STDOUT) + + if (getline(STDIN, line) == EOF) + return (EOF) + call sscan (line) + + call amovkl (long (1), v, IM_MAXDIM) + do i = 2, max (2, IM_NDIM(im)) { + call gargl (v[i]) + if (nscan() == 0) + return (EOF) + else if (nscan() != i - 1) + break + + if (IM_NDIM(im) == 1) + v[i] = 1 + else + v[i] = max (1, min (IM_LEN(im, i), v[i])) + + call sprintf (title, SZ_LINE, "%s %d") + call pargstr (title) + call pargl (v[i]) + } + + call sprintf (title, SZ_LINE, "%s\n%s") + call pargstr (title) + call pargstr (IM_TITLE(im)) + call ic_pstr (ic, "xlabel", "Column") + call gt_sets (gt, GTTITLE, title) + + return (imgnlr (im, data, v)) +end |