include include include include include "wpdef.h" ################################################################################# # SCALE -- Compute scale factors for the images. This procedure does CLIO # # to determine the type of scaling desired. The output header # # parameters for exposure time and NCOMBINE are set. The scaling # # and weighting factors are logged. The logging is done here # # because some of the information is only available here. This # # routine is taken from the `images.imcombine' package. # # # # Development version: 11/90 by RAShaw # bool procedure scale (in, out, nimages, log, low, high) include "wpcom.h" # Calling arguments: pointer in[nimages] # Input images pointer out # Output image int nimages # Number of images to combine int log # Log file descriptor real low, high # Rejection limits # Local variables: bool doscale # Perform input image scaling? int i, nout # real mean # real exposure # bool expscale # pointer exptime # Exposure time for input images pointer expname # FITS keyword for exposure time pointer fname # bool modescale # Scale image [sections] by the mode? bool modeoffset # pointer modes # pointer ncombine # Number of images previously combined pointer sp # Pointer to stack memory char str[SZ_FNAME] # Name of selected option for log pointer time # bool weight # Are images to be weighted? pointer x1, x2, xs # # Functions used: real asumr() # real asumi() # bool clgetb() # Get Boolean parameter from CL long clktime() # real imc_moder() # short imc_modes() # int imgeti() # Get integer value from image header real imgetr() # Get real value from image header begin errchk imc_moder, imc_modes, imgeti, imgetr call smark (sp) call salloc (ncombine, nimages, TY_INT) call salloc (exptime, nimages, TY_REAL) call salloc (modes, nimages, TY_REAL) call salloc (expname, SZ_FNAME, TY_CHAR) call salloc (time, SZ_DATE, TY_CHAR) call salloc (fname, SZ_FNAME, TY_CHAR) call salloc (x1, IM_MAXDIM, TY_INT) call salloc (x2, IM_MAXDIM, TY_INT) call salloc (xs, IM_MAXDIM, TY_INT) # Determine type of scaling desired and exposure time keyword. expscale = clgetb ("exposure") modescale = clgetb ("scale") modeoffset = clgetb ("offset") weight = clgetb ("weight") call clgstr ("modesec", Memc[fname], SZ_FNAME) call clgstr ("expname", Memc[expname], SZ_FNAME) call clgstr ("option", str, SZ_FNAME) # Get the number of images previously combined and the exposure times. # The default combine number is 1 and the default exposure is 0. do i = 1, nimages { iferr (Memi[ncombine+i-1] = imgeti (in[i], "ncombine")) Memi[ncombine+i-1] = 1 if (Memc[expname] != EOS) { iferr (Memr[exptime+i-1] = imgetr (in[i], Memc[expname])) Memr[exptime+i-1] = 0. } else Memr[exptime+i-1] = INDEF } # Set the default scaling factors. call amovkr (INDEF, Memr[modes], nimages) call amovkr (1., SCALES, nimages) call amovkr (0., ZEROS, nimages) call amovkr (1., WTS, nimages) # Set scaling factors. Mode scaling overrides exposure scaline and # offset scaling. if (modescale) { call amovki (1, Memi[x1], IM_NDIM(in[1])) call amovi (IM_LEN(in[1],1), Memi[x2], IM_NDIM(in[1])) call amovki (1, Memi[xs], IM_NDIM(in[1])) call subsectn (Memc[fname], Memi[x1], Memi[x2], Memi[xs], IM_NDIM(in[1])) do i = 1, nimages { switch (IM_PIXTYPE(in[i])) { case TY_SHORT: Memr[modes+i-1] = imc_modes (in[i], Memi[x1], Memi[x2], Memi[xs]) default: Memr[modes+i-1] = imc_moder (in[i], Memi[x1], Memi[x2], Memi[xs]) } SCALES[i] = Memr[modes+i-1] if (SCALES[i] <= 0.) call error (0, "Mode must be positive for scaling") if (weight) WTS[i] = sqrt (Memi[ncombine+i-1] * SCALES[i]) } } else { if (expscale) do i = 1, nimages { SCALES[i] = max (0.001, Memr[exptime+i-1]) if (weight) WTS[i] = sqrt (Memi[ncombine+i-1] * SCALES[i]) } if (modeoffset) { call amovki (1, Memi[x1], IM_NDIM(in[1])) call amovi (IM_LEN(in[1],1), Memi[x2], IM_NDIM(in[1])) call amovki (1, Memi[xs], IM_NDIM(in[1])) call subsectn (Memc[fname], Memi[x1], Memi[x2], Memi[xs], IM_NDIM(in[1])) do i = 1, nimages { switch (IM_PIXTYPE(in[i])) { case TY_SHORT: Memr[modes+i-1] = imc_modes (in[i], Memi[x1], Memi[x2], Memi[xs]) default: Memr[modes+i-1] = imc_moder (in[i], Memi[x1], Memi[x2], Memi[xs]) } ZEROS[i] = Memr[modes+i-1] / SCALES[i] if (weight) { if (ZEROS[i] <= 0.) call error (0, "Mode must be positive for weighting") WTS[i] = sqrt (Memi[ncombine+i-1]*SCALES[i]/ZEROS[i]) } } } } # Change to relative scaling factors. mean = asumr (ZEROS, nimages) / nimages call asubkr (ZEROS, mean, ZEROS, nimages) mean = asumr (SCALES, nimages) / nimages call adivkr (SCALES, mean, SCALES, nimages) call amulkr (ZEROS, mean, ZEROS, nimages) mean = asumr (WTS, nimages) call adivkr (WTS, mean, WTS, nimages) # Because of finite arithmetic it is possible for the offsets to be nonzero even # when they are all equal. Just for the sake of a nice log set the offsets # in this case. for (i=2; (i<=nimages) && (ZEROS[i] == ZEROS[1]); i=i+1) ; if (i > nimages) call aclrr (ZEROS, nimages) # If all scale factors, offsets, and weights are equal then don't actually # scale. for (i=2; (i<=nimages) && (SCALES[i] == SCALES[1]); i=i+1) if ((ZEROS[i] != ZEROS[1]) || (WTS[i] != WTS[1])) break if (i > nimages) doscale = false else doscale = true # Set the output header parameters. nout = asumi (Memi[ncombine], nimages) call imaddi (out, "ncombine", nout) if (Memc[expname] != EOS) { exposure = 0. do i = 1, nimages exposure = exposure + WTS[i] * Memr[exptime+i-1] / SCALES[i] call imaddr (out, Memc[expname], exposure) } else exposure = INDEF # Append to the log if not null. call cnvdate (clktime(0), Memc[time], SZ_DATE) if (log != NULL) { if ((low > 0.) || (high > 0.)) { call fprintf (log, "%s combine: %s, lowreject=%g, highreject=%g\n") call pargstr (Memc[time]) call pargstr (str) call pargr (low) call pargr (high) } else { call fprintf (log, "%s combine: %s\n") call pargstr (Memc[time]) call pargstr (str) } if (doscale) { call fprintf (log, " %20s %6s %6s %6s %6s %6s %6s\n") call pargstr ("Images") call pargstr ("N") call pargstr ("Exp") call pargstr ("Mode") call pargstr ("Scale") call pargstr ("Offset") call pargstr ("Weight") do i = 1, nimages { call imstats (in[i], IM_IMAGENAME, Memc[fname], SZ_FNAME) call fprintf (log, " %20s %6d %6.1f %6g %6.3f %6g %6.3f\n") call pargstr (Memc[fname]) call pargi (Memi[ncombine+i-1]) call pargr (Memr[exptime+i-1]) call pargr (Memr[modes+i-1]) call pargr (1./SCALES[i]) call pargr (-ZEROS[i]) call pargr (WTS[i]) } call fprintf (log, " -------------------- ------ ------\n") call imstats (out, IM_IMAGENAME, Memc[fname], SZ_FNAME) call fprintf (log, " %20s %6d %6.1f\n") call pargstr (Memc[fname]) call pargi (nout) call pargr (exposure) } else { call fprintf (log, " %20s %6s\n") call pargstr ("Images") call pargstr ("N") do i = 1, nimages { call imstats (in[i], IM_IMAGENAME, Memc[fname], SZ_FNAME) call fprintf (log, " %20s %6d\n") call pargstr (Memc[fname]) call pargi (Memi[ncombine+i-1]) } call fprintf (log, " -------------------- ------\n") call imstats (out, IM_IMAGENAME, Memc[fname], SZ_FNAME) call fprintf (log, " %20s %6d\n") call pargstr (Memc[fname]) call pargi (nout) } call flush (log) } call sfree (sp) return (doscale) end