diff options
Diffstat (limited to 'noao/twodspec/multispec/fitsmooth.x')
| -rw-r--r-- | noao/twodspec/multispec/fitsmooth.x | 168 |
1 files changed, 168 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/fitsmooth.x b/noao/twodspec/multispec/fitsmooth.x new file mode 100644 index 00000000..413f3c46 --- /dev/null +++ b/noao/twodspec/multispec/fitsmooth.x @@ -0,0 +1,168 @@ + +# FIT_SMOOTH -- Least-squares fit of smoothed profiles to data profiles with +# cleaning of deviant pixels. +# +# The profile fitting and cleaning are combined in order to minimize +# the calculations in re-evaluating the least-squares fit after rejecting +# deviant pixels. +# +# The sigma used for rejection is calculated from the sigma of the fit +# before rejecting any pixels. Pixels whose residuals exceed +# +/- sigma_cut * sigma are rejected. The maximum number of pixels to be +# replaced in each spectrum is max_replace. If max_replace is zero then +# only the model fitting is performed. +# +# The output of this routine are the cleaned data profiles and the least-square +# fitted model profiles. The number of pixels replaced is returned. + + +procedure fit_smooth (line, data, model, profiles, len_prof, nspectra, nlines) + +int line # Image line of data +real data[len_prof, nspectra] # Data profiles +real model[len_prof, nspectra] # Model profiles +real profiles[len_prof, nspectra, ARB] # Work array for SMOOTH profiles +int len_prof # Length of profile +int nspectra # Number of spectra +int nlines # Number of lines profiles + +int max_replace # Maximum number of bad pixels +real sigma_cut # Sigma cutoff on the residuals + +int i, spectrum +int nmax, ntotal, nreplace, nreject, nindef, nsigma +real sum1, sum2, scale, sigma +real lower, upper, residual, resid_min, resid_max +pointer sp, a, b, c + +begin + # Allocate working memory. + call smark (sp) + call salloc (a, nspectra, TY_REAL) + call salloc (b, nspectra, TY_REAL) + call salloc (c, nspectra, TY_INT) + + # Fit each spectrum and compute sigma of fit. + sigma = 0. + nsigma = 0 + do spectrum = 1, nspectra { + # Accumulate least squares sums. + sum1 = 0. + sum2 = 0. + nindef = 0 + do i = 1, len_prof { + if (IS_INDEFR (data[i, spectrum])) + nindef = nindef + 1 + else if (model[i, spectrum] > 0.) { + sum1 = sum1 + data[i, spectrum] * model[i, spectrum] + sum2 = sum2 + model[i, spectrum] * model[i, spectrum] + } + } + + # Compute sigma if cleanning is desired. + if (nmax != 0) { + scale = sum1 / sum2 + do i = 1, len_prof { + if (!IS_INDEFR (data[i, spectrum]) && + (model[i, spectrum] > 0.)) { + sigma = sigma + + (data[i,spectrum] - scale * model[i,spectrum]) ** 2 + nsigma = nsigma + 1 + } + } + } + + Memr[a + spectrum - 1] = sum1 + Memr[b + spectrum - 1] = sum2 + Memi[c + spectrum - 1] = nindef + } + sigma = sqrt (sigma / nsigma) + + # Reject deviant pixels from the fit, scale the model to data, + # and replace rejected and INDEFR pixels with model values. + ntotal = 0 + do spectrum = 1, nspectra { + sum1 = Memr[a + spectrum - 1] + sum2 = Memr[b + spectrum - 1] + nindef = Memi[c + spectrum - 1] + + # If there are no model data points go to the next spectrum. + if (sum2 == 0.) + next + + # Reject pixels if desired. + nreplace = 0 + if (nmax != 0) { + # Compare each pixel in the profile against the model and set + # deviant pixels to INDEFR. If the number of pixels to be + # replaced is equal to the maximum allowed or the number of + # pixels rejected equals the entire profile or the number of + # deviant pixels is zero in an iteration stop cleaning and + # exit the loop. Ignore INDEFR pixels. + + repeat { + nreject = 0 + scale = sum1 / sum2 + resid_min = -lower * sigma + resid_max = upper * sigma + do i = 1, len_prof { + if (IS_INDEFR (data[i, spectrum])) + next + + # Compute the residual and remove point if it exceeds + # the residual limits. + + residual = data[i,spectrum] - scale * model[i,spectrum] + if ((residual < resid_min) || (residual > resid_max)) { + # Remove point from the least squares fit + # and flag the deviant pixel with INDEFR. + sum1 = sum1 - data[i,spectrum] * model[i,spectrum] + sum2 = sum2 - model[i,spectrum] * model[i,spectrum] + data[i,spectrum] = INDEFR + nreplace = nreplace + 1 + nreject = nreject + 1 + } + if (nreplace == nmax) + break + } + } until ((nreplace == nmax) || (nreject == 0) || (sum2 == 0.)) + } + + # If there are good pixels remaining scale the model to the + # data profile. + if (sum2 > 0.) + call amulkr (model[1, spectrum], sum1 / sum2, + model[1, spectrum], len_prof) + + # Replace bad pixels by the model values. + if ((nindef > 0) || (nreplace > 0)) { + do i = 1, len_prof { + if (IS_INDEFR (data[i, spectrum])) + data[i, spectrum] = model[i, spectrum] + } + ntotal = ntotal + nreplace + } + } + + # Print the number of pixel replaced. + call ex_prnt3 (ntotal) + + # Replace the cleaned data profiles in future SMOOTH profiles. + if (ntotal > 0) + call update_smooth (line, data, profiles, len_prof, nspectra, + nlines) + + # Free allocated memory. + call sfree (sp) + + return + +# SET_FIT_SMOOTH -- Set the cleaning parameters. + +entry set_fit_smooth (max_replace, sigma_cut) + + nmax = max_replace + lower = sigma_cut + upper = sigma_cut + return +end |