From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/proto/masks/rsstats.x | 492 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 492 insertions(+) create mode 100644 pkg/proto/masks/rsstats.x (limited to 'pkg/proto/masks/rsstats.x') diff --git a/pkg/proto/masks/rsstats.x b/pkg/proto/masks/rsstats.x new file mode 100644 index 00000000..9c7a1b32 --- /dev/null +++ b/pkg/proto/masks/rsstats.x @@ -0,0 +1,492 @@ +include +include +include +include +include "mimstat.h" +include "rskysub.h" + + +# RS_STATS -- Compute the input image scaling factors. + +procedure rs_stats (inlist, imsklist, omsklist, sclist, rs, msk_invert, + cache, verbose) + +int inlist #I the input image list +int imsklist #I the input mask list +int omsklist #I the output mask list +int sclist #I the input scale factors list +pointer rs #I the sky subtraction descriptor +bool msk_invert #I invert the pixel masks ? +bool cache #I cache the image i/o buffers ? +bool verbose #I print image statistics ? + +real fscale +pointer sp, image, imaskname, omaskname, masktemp, str +pointer im, ims, pmim, pmout +int ip, old_size +real imgetr() +pointer immap(), im_pmmap(), mp_open() +int imtgetim(), imtlen(), imtrgetim(), ctor(), ctowrd(), btoi() +int fntgfnb(), imstati(), imaccess() +bool strne(), streq() +errchk immap() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (imaskname, SZ_FNAME, TY_CHAR) + call salloc (omaskname, SZ_FNAME, TY_CHAR) + call salloc (masktemp, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Loop over the input images and compute the scale factors. + # At some point we might combine this with the later running + # mean / median loop for more efficient operation especially in an + # observing environment but that can easily be rearranged later. + + while (imtgetim (inlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the input image. This image is opened READ_WRITE + # so some header information can be added ... + iferr (im = immap (Memc[image], READ_WRITE, 0)) { + call printf ("Error opening image %s ...\n") + call pargstr (Memc[image]) + next + } + + # Check for a statistics section. If the image image already + # includes a section strip it off, append the statistics + # section to the input image name, and open the statistics + # section image. + + if (streq (RS_ISCALES(rs), "median") && RS_STATSEC(rs) != EOS) { + call imgimage (Memc[image], Memc[str], SZ_FNAME) + call strcat (RS_STATSEC(rs), Memc[str], SZ_FNAME) + iferr (ims = immap (Memc[str], READ_ONLY, 0)) { + call imunmap (im) + call printf ("Error opening image %s ...\n") + call pargstr (Memc[image]) + next + } + } else + ims = NULL + + # Open input the mask if any. The input and output mask + # lists are ignored if the scaling factor is not median + # or if the list lengths are both zero. + if (strne (RS_ISCALES(rs), "median")) { + pmim = NULL + pmout = NULL + } else if (imtlen (omsklist) == 0 && imtlen (imsklist) == 0) { + pmim = NULL + pmout = NULL + } else { + + # Get the input mask which defaults to the empty mask if + # there is none. + if (imtgetim (imsklist, Memc[str+1], SZ_FNAME) != EOF) { + if (msk_invert) { + Memc[str] = '^' + pmim = mp_open (Memc[str], im, Memc[imaskname], + SZ_FNAME) + } else + pmim = mp_open (Memc[str+1], im, Memc[imaskname], + SZ_FNAME) + } else if (imtrgetim (imsklist, 1, Memc[str], + SZ_FNAME) != EOF) { + pmim = mp_open (Memc[str], im, Memc[imaskname], SZ_FNAME) + } else { + pmim = mp_open ("", im, Memc[imaskname], SZ_FNAME) + } + if (pmim == NULL) { + call printf ("Error reading mask for image %s ...\n") + call pargstr (Memc[image]) + call imunmap (im) + next + } + + # Get the output mask name if any. + if (imtlen (omsklist) > 0) { + if (imtgetim (omsklist, Memc[omaskname], SZ_FNAME) == EOF) { + call imunmap (pmim) + call imunmap (im) + next + } else { + if (Memc[imaskname] == '^') + call xt_mkimtemp (Memc[imaskname+1], + Memc[omaskname], Memc[masktemp], SZ_FNAME) + else + call xt_mkimtemp (Memc[imaskname], Memc[omaskname], + Memc[masktemp], SZ_FNAME) + pmout = im_pmmap (Memc[omaskname], NEW_IMAGE, 0) + call mp_mpcopy (im, pmim, pmout) + } + } else { + pmout = NULL + } + } + + + # Print title. + if (verbose) { + if (pmim == NULL) { + call printf ("Computing scale factor for image %s\n") + call pargstr (Memc[image]) + } else { + call printf ( + "Computing scale factor for image %s using mask %s\n") + call pargstr (Memc[image]) + call pargstr (Memc[imaskname]) + } + call flush (STDOUT) + } + + # Check for existence of scaling keyword. If the keyword is + # present and the rescaling flag is turned off then proceed + # to the next image, otherwise compute the new scale factor. + + if (RS_RESCALE(rs) == NO) { + ifnoerr (fscale = imgetr (im, RS_KYFSCALE(rs))) { + if (verbose) { + call printf (" Using precomputed value %g\n") + call pargr (fscale) + } + call imunmap (pmim) + if (ims != NULL) + call imunmap (ims) + call imunmap (im) + next + } + } + + # Compute the scaling factor. The scaling factor defaults + # to, 1 if the scaling method is "none", the value of the image + # header keyowrd if the scaling factor is !KEYWORD, 1 / median + # if the the scaling methid is "median", or the value in the + # scaling factors file if the scaling factor is "@file". If an + # error occurs the scaling factor is set to 1.0. + + if (streq (RS_ISCALES(rs), "none")) { + fscale = 1.0 + } else if (RS_ISCALES(rs) == '!') { + ip = 2 + if (ctowrd (RS_ISCALES(rs), ip, Memc[str], SZ_FNAME) <= 0) + Memc[str] = EOS + iferr (fscale = imgetr (im, Memc[str])) + fscale = 1.0 + } else if (streq (RS_ISCALES(rs), "median")) { + if (ims != NULL) + call rs_cache1 (btoi(cache), ims, old_size) + else + call rs_cache1 (btoi(cache), im, old_size) + if (pmim == NULL) { + if (ims != NULL) + call rs_med (ims, rs, fscale) + else + call rs_med (im, rs, fscale) + } else { + if (ims != NULL) + call rs_mmed (im, ims, pmim, pmout, rs, fscale) + else + call rs_mmed (im, im, pmim, pmout, rs, fscale) + } + if (IS_INDEFR(fscale)) + fscale = 1.0 + else + fscale = 1.0 / fscale + call fixmem (old_size) + } else if (fntgfnb (sclist, Memc[str], SZ_FNAME) != EOF) { + ip = 1 + if (ctor (Memc[str], ip, fscale) <= 0) + fscale = 1.0 + } else { + fscale = 1.0 + } + + # Print the computed scaling factor. + if (verbose) { + call printf (" New scale factor is 1 / %g\n") + call pargr (1.0 / fscale) + if (pmout != NULL) { + call printf (" Writing new image mask %s\n") + call pargstr (Memc[masktemp]) + } + call flush (STDOUT) + } + + # Store the new scaling factor in the input image header. + call imaddr (im, RS_KYFSCALE(rs), fscale) + + # Close the input image and mask. + if (pmout != NULL) { + if (imaccess (Memc[omaskname], YES) == YES) + call imdelete (Memc[omaskname]) + call pm_savef (imstati (pmout, IM_PMDES), Memc[omaskname], + "", 0) + call imunmap (pmout) + if (pmim != NULL) + call imunmap (pmim) + call xt_delimtemp (Memc[omaskname], Memc[masktemp]) + } else { + if (pmim != NULL) + call imunmap (pmim) + } + if (ims != NULL) + call imunmap (ims) + call imunmap (im) + } + + call sfree (sp) +end + + +# RS_MMED -- Estimate the image median using iterative rejection and +# a pixel mask. The input image and input statistics image descriptors may +# be the same. + +procedure rs_mmed (im, ims, pmim, pmout, rs, fscale) + +pointer im #I the input image descriptor +pointer ims #I the input image statistics descriptor +pointer pmim #I the input mask image descriptor +pointer pmout #I the output mask image descriptor +pointer rs #I the sky subtraction pointer +real fscale #O the scaling factor + +real low, up, hmin, hmax, hwidth +pointer sp, vs, ve, mst, pm, mp, buf, hgm, smsk +int i, mval, npts, npix, nbins, nbad + +pointer mp_miopen() +int imstati(), mio_glsegr(), mst_ihist(), rs_umask() + +begin + call smark (sp) + call salloc (vs, IM_MAXDIM, TY_LONG) + call salloc (ve, IM_MAXDIM, TY_LONG) + + # Allocate space for statistics structure. + call mst_allocate (mst) + + # Get the selected fields. + #nfields = mst_fields ("midpt,stddev" Memi[fields], MIS_NFIELDS) + + # Set the processing switches + #call mst_switches (mst, Memi[fields], nfields, RS_MAXITER(rs)) + + # Set up the region masking parameters. + mp = mp_miopen (ims, pmim) + + # Compute the image statistics. + low = RS_LOWER(rs) + up = RS_UPPER(rs) + do i = 0 , RS_MAXITER(rs) { + + # Set up the mask i/o boundaries. + call amovkl (long(1), Meml[vs], IM_NDIM(ims)) + call amovl (IM_LEN(ims,1), Meml[ve], IM_NDIM(ims)) + call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(ims)) + + # Initialize the statistics computation. + call mst_initialize (mst, low, up) + + # Accumulate the statistics. + while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) + call mst_accumulate2 (mst, Memr[buf], npts, low, up, YES) + + # Compute the 2nd order central moment statistics. + call mst_stats (mst) + + # Compute new limits and iterate. + if (i < RS_MAXITER(rs)) { + if (IS_INDEFR(RS_LNSIGREJ(rs))) + low = -MAX_REAL + else if (RS_LNSIGREJ(rs) > 0.0 || IS_INDEFR(MIS_MEAN(mst)) || + IS_INDEFR(MIS_STDDEV(mst))) + low = MIS_MEAN(mst) - RS_LNSIGREJ(rs) * MIS_STDDEV(mst) + else + low = -MAX_REAL + if (IS_INDEFR(RS_UNSIGREJ(rs))) + up = MAX_REAL + else if (RS_UNSIGREJ(rs) > 0.0 || IS_INDEFR(MIS_MEAN(mst)) || + IS_INDEFR(MIS_STDDEV(mst))) + up = MIS_MEAN(mst) + RS_UNSIGREJ(rs) * MIS_STDDEV(mst) + else + up = MAX_REAL + if (i > 0) { + if (MIS_NPIX(mst) == npix) + break + } + npix = MIS_NPIX(mst) + } + + } + + # Estimate the median and the mode by accumulating the histogram. + hgm = NULL + if (mst_ihist (mst, RS_BINWIDTH(rs), hgm, nbins, hwidth, hmin, + hmax) == YES) { + call aclri (Memi[hgm], nbins) + call amovkl (long(1), Meml[vs], IM_NDIM(ims)) + call amovl (IM_LEN(ims,1), Meml[ve], IM_NDIM(ims)) + call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(ims)) + while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) + call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax) + call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin, hmax) + } + if (hgm != NULL) + call mfree (hgm, TY_INT) + + # Set the statistic + fscale = MIS_MEDIAN(mst) + + if (pmout != NULL) { + call malloc (smsk, IM_LEN(im,1), TY_SHORT) + call amovkl (long(1), Meml[vs], IM_NDIM(im)) + call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im)) + call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im)) + pm = imstati (pmout, IM_PMDES) + while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) { + nbad = rs_umask (Memr[buf], Mems[smsk], npts, low, up) + if (nbad > 0) + call pm_plps (pm, Meml[vs], Mems[smsk], 1, npts, PIX_SRC) + } + call mp_invert (pm) + call imseti (pmout, IM_PMDES, pm) + call mfree (smsk, TY_SHORT) + } + + # Close the maskio descriptor. + call mio_close (mp) + + call mst_free (mst) + call sfree (sp) +end + + +# RS_MED -- Estimate the image median using iterative rejection and +# no pixel mask. + +procedure rs_med (im, rs, fscale) + +pointer im #I the input image descriptor +pointer rs #I the sky subtraction descriptor +real fscale #I the computed scaling factor + + +real low, up, hmin, hmax, hwidth +pointer sp, v, mst, buf, hgm +int i, npts, npix, nbins +int imgnlr(), mst_ihist() + +begin + call smark (sp) + call salloc (v, IM_MAXDIM, TY_LONG) + + # Allocate space for statistics structure. + call mst_allocate (mst) + + # Get the selected fields. + #nfields = mst_fields ("midpt,stddev" Memi[fields], MIS_NFIELDS) + + # Set the processing switches + #call mst_switches (mst, Memi[fields], nfields, RS_MAXITER(rs)) + + # Compute the image statistics. + low = RS_LOWER(rs) + up = RS_UPPER(rs) + do i = 0 , RS_MAXITER(rs) { + + # Initialize the statistics computation. + call mst_initialize (mst, low, up) + + # Accumulate the statistics. + npts = IM_LEN(im,1) + call amovkl (long(1), Meml[v], IM_NDIM(im)) + while (imgnlr (im, buf, Meml[v]) != EOF) + call mst_accumulate2 (mst, Memr[buf], npts, low, up, YES) + + # Compute the 2nd order central moment statistics. + call mst_stats (mst) + + # Compute new limits and iterate. + if (i < RS_MAXITER(rs)) { + if (IS_INDEFR(RS_LNSIGREJ(rs))) + low = -MAX_REAL + else if (RS_LNSIGREJ(rs) > 0.0) + low = MIS_MEAN(mst) - RS_LNSIGREJ(rs) * MIS_STDDEV(mst) + else + low = -MAX_REAL + if (IS_INDEFR(RS_UNSIGREJ(rs))) + up = MAX_REAL + else if (RS_UNSIGREJ(rs) > 0.0) + up = MIS_MEAN(mst) + RS_UNSIGREJ(rs) * MIS_STDDEV(mst) + else + up = MAX_REAL + if (i > 0) { + if (MIS_NPIX(mst) == npix) + break + } + npix = MIS_NPIX(mst) + } + + } + + # Estimate the median and the mode by accumulating the histogram. + hgm = NULL + if (mst_ihist (mst, RS_BINWIDTH(rs), hgm, nbins, hwidth, hmin, + hmax) == YES) { + call aclri (Memi[hgm], nbins) + call amovkl (long(1), Meml[v], IM_NDIM(im)) + while (imgnlr (im, buf, Meml[v]) != EOF) + call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax) + call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin, hmax) + } + if (hgm != NULL) + call mfree (hgm, TY_INT) + + # Set the statistic + fscale = MIS_MEDIAN(mst) + + call mst_free (mst) + call sfree (sp) +end + + +# RS_UMASK -- Update the mask. + +int procedure rs_umask (pix, msk, npts, lower, upper) + +real pix[ARB] #I array of image pixels +short msk[ARB] #O array of mask pixels, set to 1 and 0 +int npts #I the number of pixels +real lower #I the lower good data limit +real upper #I the upper good data limit + +real lo, up +int i, nbad + +begin + if (IS_INDEFR(lower) && IS_INDEFR(upper)) + return (0) + + if (IS_INDEFR(lower)) + lo = -MAX_REAL + else + lo = lower + if (IS_INDEFR(upper)) + up = MAX_REAL + else + up = upper + + nbad = 0 + do i = 1, npts { + if (pix[i] < lo || pix[i] > up) { + msk[i] = 0 + nbad = nbad + 1 + } else + msk[i] = 1 + } + + return (nbad) +end -- cgit