Index: /branches/2dbias/psLib/src/math/psVectorSmooth.c =================================================================== --- /branches/2dbias/psLib/src/math/psVectorSmooth.c (revision 42714) +++ /branches/2dbias/psLib/src/math/psVectorSmooth.c (revision 42715) @@ -1,4 +1,4 @@ #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif @@ -15,4 +15,5 @@ #include "psAssert.h" #include "psConstants.h" +#include "psStats.h" #include "psVectorSmooth.h" @@ -21,6 +22,5 @@ const psVector *input, double sigma, - double Nsigma - ) + double Nsigma) { PS_ASSERT_VECTOR_NON_NULL(input, NULL); @@ -28,5 +28,6 @@ PS_ASSERT_FLOAT_LARGER_THAN(Nsigma, 0.0, NULL); - if (output == input) { + if (output == input) + { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot smooth vector in-place."); return NULL; @@ -34,143 +35,175 @@ // relevant terms - long Nrange = sigma*Nsigma + 0.5; // Extent of smoothing - long Npixel = 2*Nrange + 1; // Total size of smoothing kernel - long Nbin = input->n; // Number of elements - double factor = -0.5/(sigma*sigma); // Factor for Gaussian - - if (Nbin < Npixel) { - // cannot smooth narrow vector - return NULL; - } - - #define VECTOR_SMOOTH_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - output = psVectorRecycle(output, Nbin, PS_TYPE_##TYPE); \ - /* generate normalized gaussian */ \ - psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \ - double sum = 0.0; \ - for (long i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] = exp(factor*i*i); \ - sum += gaussnorm->data.TYPE[i+Nrange]; \ - } \ - for (long i = -Nrange; i < Nrange + 1; i++) { \ - gaussnorm->data.TYPE[i+Nrange] /= sum; \ - } \ - ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \ - \ - /* smooth vector */ \ - psVector *temp = psVectorAlloc(Nbin, PS_TYPE_##TYPE); \ - ps##TYPE *vi = input->data.TYPE; \ - ps##TYPE *vo = temp->data.TYPE; \ - /* smooth first Nrange pixels, with renorm */ \ + long Nrange = sigma * Nsigma + 0.5; // Extent of smoothing + long Npixel = 2 * Nrange + 1; // Total size of smoothing kernel + long Nbin = input->n; // Number of elements + double factor = -0.5 / (sigma * sigma); // Factor for Gaussian + + if (Nbin < Npixel) + { + // cannot smooth narrow vector + return NULL; + } + +#define VECTOR_SMOOTH_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + { \ + output = psVectorRecycle(output, Nbin, PS_TYPE_##TYPE); \ + /* generate normalized gaussian */ \ + psVector *gaussnorm = psVectorAlloc(Npixel, PS_TYPE_##TYPE); \ + double sum = 0.0; \ + for (long i = -Nrange; i < Nrange + 1; i++) \ + { \ + gaussnorm->data.TYPE[i + Nrange] = exp(factor * i * i); \ + sum += gaussnorm->data.TYPE[i + Nrange]; \ + } \ + for (long i = -Nrange; i < Nrange + 1; i++) \ + { \ + gaussnorm->data.TYPE[i + Nrange] /= sum; \ + } \ + ps##TYPE *gauss = &gaussnorm->data.TYPE[Nrange]; \ + \ + /* smooth vector */ \ + psVector *temp = psVectorAlloc(Nbin, PS_TYPE_##TYPE); \ + ps##TYPE *vi = input->data.TYPE; \ + ps##TYPE *vo = temp->data.TYPE; \ + /* smooth first Nrange pixels, with renorm */ \ /* XXX need to check that this does not run over end for narrow vectors */ \ - for (long i = 0; i < Nrange; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - i; \ - ps##TYPE *vg = gauss - i; \ - double g = 0; \ - double s = 0; \ - for (int n = -i; n < Nrange + 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - g += *vg; \ - } \ - *vo = s / g; \ - } \ - /* smooth middle pixels */ \ - for (long i = Nrange; i < Nbin - Nrange; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - Nrange; \ - ps##TYPE *vg = gauss - Nrange; \ - double s = 0; \ - for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - } \ - *vo = s; \ - } \ - /* smooth last Nrange pixels, with renorm */ \ - /* XXX does this miss the last column? */ \ - for (long i = Nbin - Nrange; i < Nbin; i++, vi++, vo++) { \ - ps##TYPE *vr = vi - Nrange; \ - ps##TYPE *vg = gauss - Nrange; \ - double g = 0; \ - double s = 0; \ - for (int n = -Nrange; n < Nbin - i - 1; n++, vr++, vg++) { \ - s += *vg * *vr; \ - g += *vg; \ - } \ - *vo = s / g; \ - } \ - memcpy(output->data.TYPE, temp->data.TYPE, Nbin*sizeof(ps##TYPE)); \ - psFree(temp); \ - psFree(gaussnorm); \ - break; \ - } - - switch (input->type.type) { + for (long i = 0; i < Nrange; i++, vi++, vo++) \ + { \ + ps##TYPE *vr = vi - i; \ + ps##TYPE *vg = gauss - i; \ + double g = 0; \ + double s = 0; \ + for (int n = -i; n < Nrange + 1; n++, vr++, vg++) \ + { \ + s += *vg * *vr; \ + g += *vg; \ + } \ + *vo = s / g; \ + } \ + /* smooth middle pixels */ \ + for (long i = Nrange; i < Nbin - Nrange; i++, vi++, vo++) \ + { \ + ps##TYPE *vr = vi - Nrange; \ + ps##TYPE *vg = gauss - Nrange; \ + double s = 0; \ + for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) \ + { \ + s += *vg * *vr; \ + } \ + *vo = s; \ + } \ + /* smooth last Nrange pixels, with renorm */ \ + /* XXX does this miss the last column? */ \ + for (long i = Nbin - Nrange; i < Nbin; i++, vi++, vo++) \ + { \ + ps##TYPE *vr = vi - Nrange; \ + ps##TYPE *vg = gauss - Nrange; \ + double g = 0; \ + double s = 0; \ + for (int n = -Nrange; n < Nbin - i - 1; n++, vr++, vg++) \ + { \ + s += *vg * *vr; \ + g += *vg; \ + } \ + *vo = s / g; \ + } \ + memcpy(output->data.TYPE, temp->data.TYPE, Nbin * sizeof(ps##TYPE)); \ + psFree(temp); \ + psFree(gaussnorm); \ + break; \ + } + + switch (input->type.type) + { VECTOR_SMOOTH_CASE(F32); VECTOR_SMOOTH_CASE(F64); - default: { - char* typeStr; - PS_TYPE_NAME(typeStr, input->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type %s is not valid.", typeStr); - return NULL; - } + default: + { + char *typeStr; + PS_TYPE_NAME(typeStr, input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type %s is not valid.", typeStr); + return NULL; + } } return output; } - - psVector *psVectorBoxcar(psVector *output, const psVector *input, - int size - ) + int size, + bool robust) { PS_ASSERT_VECTOR_NON_NULL(input, NULL); PS_ASSERT_INT_POSITIVE(size, NULL); - if (output == input) { + if (output == input) + { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot smooth vector in-place."); return false; } - long num = input->n; // Number of elements + long num = input->n; // Number of elements output = psVectorRecycle(output, num, input->type.type); - psVector *nums = psVectorAlloc(num, PS_TYPE_U32); // Number of elements in each bin - psU32 *numsData = nums->data.U32; // Dereferenced version - psVectorInit(output, 0.0); - psVectorInit(nums, 0); - - - #define VECTOR_BOXCAR_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - /* Dereference data */ \ - ps##TYPE *outputData = output->data.TYPE; \ - ps##TYPE *inputData = input->data.TYPE; \ - /* Smooth the vector */ \ - for (long i = 0; i < num; i++) { \ - for (long j = PS_MAX(0, i - size); j < PS_MIN(num, i + size + 1); j++) { \ - outputData[j] += inputData[i]; \ - numsData[j]++; \ - } \ - } \ - /* Normalisation */ \ - for (long i = 0; i < num; i++) { \ - outputData[i] /= numsData[i]; \ - } \ - break; \ - } - - switch (input->type.type) { + +#define VECTOR_BOXCAR_CASE(TYPE) \ + case PS_TYPE_##TYPE: \ + { \ + /* Dereference data */ \ + ps##TYPE *outputData = output->data.TYPE; \ + ps##TYPE *inputData = input->data.TYPE; \ + /* Smooth the vector */ \ + for (long i = 0; i < num; i++) \ + { \ + /* sample within the window size */ \ + int windowSize = PS_MIN(num, i + size + 1) - PS_MAX(0, i - size); \ + psVector *sample = psVectorAlloc(windowSize, PS_TYPE_##TYPE); \ + ps##TYPE *sampleData = sample->data.TYPE; \ + /* extract sample from inputData */ \ + for (long j = PS_MAX(0, i - size); j < PS_MIN(num, i + size + 1); j++) \ + { \ + sampleData[j - PS_MAX(0, i - size)] = inputData[j]; \ + } \ + /* find the desired statistic from the sample */ \ + psStatsOptions statistic; \ + psStats *stats; \ + if (robust) \ + { \ + statistic = PS_STAT_CLIPPED_MEAN; \ + stats = psStatsAlloc(statistic); \ + } \ + else \ + { \ + statistic = PS_STAT_SAMPLE_MEDIAN; \ + stats = psStatsAlloc(statistic); \ + } \ + if (!psVectorStats(stats, sample, NULL, NULL, 0)) \ + { \ + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); \ + psFree(sample); \ + psFree(stats); \ + return NULL; \ + } \ + outputData[i] = psStatsGetValue(stats, statistic); \ + psFree(sample); \ + psFree(stats); \ + } \ + break; \ + } + + switch (input->type.type) + { VECTOR_BOXCAR_CASE(F32); VECTOR_BOXCAR_CASE(F64); - default: { - char* typeStr; - PS_TYPE_NAME(typeStr, input->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type %s is not valid.", typeStr); - psFree(nums); - return NULL; - } - } - psFree(nums); + default: + { + char *typeStr; + PS_TYPE_NAME(typeStr, input->type.type); + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type %s is not valid.", typeStr); + return NULL; + } + } return output; } Index: /branches/2dbias/psLib/src/math/psVectorSmooth.h =================================================================== --- /branches/2dbias/psLib/src/math/psVectorSmooth.h (revision 42714) +++ /branches/2dbias/psLib/src/math/psVectorSmooth.h (revision 42715) @@ -14,15 +14,16 @@ /// Smooth a vector with a Gaussian -psVector *psVectorSmooth(psVector *output, ///< Output vector, or NULL +psVector *psVectorSmooth(psVector *output, ///< Output vector, or NULL const psVector *input, ///< Input vector (F32 or F64 only) - double sigma, ///< Gausian width (standard deviations) - double Nsigma ///< Number of standard deviations for Gaussian to extend - ); + double sigma, ///< Gausian width (standard deviations) + double Nsigma ///< Number of standard deviations for Gaussian to extend +); /// Smooth a vector with a boxcar -psVector *psVectorBoxcar(psVector *output, ///< Output vector, or NULL +psVector *psVectorBoxcar(psVector *output, ///< Output vector, or NULL const psVector *input, ///< Input vector (F32 or F64 only) - int size ///< Boxcar size (one-sided size) - ); + int size, ///< Boxcar size (one-sided size) + bool robust ///< compute robust median or normal median +); /// @} Index: /branches/2dbias/psModules/src/detrend/pmOverscan.c =================================================================== --- /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42714) +++ /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42715) @@ -209,5 +209,5 @@ // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels, false); psFree(pixels); if (!reduced) @@ -293,5 +293,5 @@ // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels); + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels, false); psFree(pixels); if (!reduced) @@ -415,5 +415,5 @@ // Reduce the overscans // XXX need to save 2 different chi-square values - psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels); + psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels, false); psFree(yscanPixels); if (!yReduced) @@ -430,5 +430,5 @@ // Reduce the overscans // XXX need to save 2 different chi-square values - psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels); + psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels, true); psFree(xscanPixels); if (!xReduced) @@ -560,6 +560,6 @@ psVector *pmOverscanVector(float *chi2, // chi^2 from fit pmOverscanStatOptions *overscanOpts, // Overscan statistic options - const psArray *pixels // Array of vectors containing the pixel values -) + const psArray *pixels, // Array of vectors containing the pixel values + bool robust) // Use robust statistics for boxcar smoothing { assert(overscanOpts); @@ -611,5 +611,5 @@ if (overscanOpts->boxcar > 0) { - psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector + psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar, robust); // Smoothed vector psFree(reduced); reduced = smoothed; Index: /branches/2dbias/psModules/src/detrend/pmOverscan.h =================================================================== --- /branches/2dbias/psModules/src/detrend/pmOverscan.h (revision 42714) +++ /branches/2dbias/psModules/src/detrend/pmOverscan.h (revision 42715) @@ -18,9 +18,10 @@ /// Type of fit to perform -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLY_ORD, ///< Fit ordinary polynomial - PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial - PM_FIT_SPLINE ///< Fit cubic splines +typedef enum +{ + PM_FIT_NONE, ///< No fit + PM_FIT_POLY_ORD, ///< Fit ordinary polynomial + PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial + PM_FIT_SPLINE ///< Fit cubic splines } pmFit; @@ -35,14 +36,14 @@ { // Inputs - pmFit fitType; ///< Type of fit to overscan - unsigned int order; ///< Order of polynomial, or number of spline pieces - psStats *stat; ///< Statistic to use when reducing the minor direction - int boxcar; ///< Boxcar smoothing radius - float gauss; ///< Gaussian smoothing sigma + pmFit fitType; ///< Type of fit to overscan + unsigned int order; ///< Order of polynomial, or number of spline pieces + psStats *stat; ///< Statistic to use when reducing the minor direction + int boxcar; ///< Boxcar smoothing radius + float gauss; ///< Gaussian smoothing sigma // These are used to carry the fitted functions, but are only used to generate - // an updated reduced vector - psPolynomial1D *poly; ///< Result of polynomial fit - psSpline1D *spline; ///< Result of spline fit + // an updated reduced vector + psPolynomial1D *poly; ///< Result of polynomial fit + psSpline1D *spline; ///< Result of spline fit } pmOverscanStatOptions; @@ -57,11 +58,11 @@ { // Inputs - bool single; ///< Reduce all overscan regions to a single value? - bool constant; ///< use a supplied constant value (do not measure region) - bool TwoD; ///< use a supplied constant value (do not measure region) - float value; ///< supplied value if needed (per above) - float minValid; ///< if overscan is too low, readout is dead : mask - float maxValid; ///< if overscan is too high, readout is dead : mask - psImageMaskType maskVal; ///< Mask value to give dead readouts + bool single; ///< Reduce all overscan regions to a single value? + bool constant; ///< use a supplied constant value (do not measure region) + bool TwoD; ///< use a supplied constant value (do not measure region) + float value; ///< supplied value if needed (per above) + float minValid; ///< if overscan is too low, readout is dead : mask + float maxValid; ///< if overscan is too high, readout is dead : mask + psImageMaskType maskVal; ///< Mask value to give dead readouts pmOverscanStatOptions *primary; @@ -74,14 +75,14 @@ pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void); -psVector *pmOverscanVector(float *chi2, // chi^2 from fit - pmOverscanStatOptions *overscanOpts, // Overscan options - const psArray *pixels // Array of vectors containing the pixel values - ); +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanStatOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + bool robust // use robust statistics for boxcar smoothing +); // bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2); -bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan); +bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan); /// @} #endif -