Index: /trunk/psLib/src/imageops/psImageCovariance.c =================================================================== --- /trunk/psLib/src/imageops/psImageCovariance.c (revision 21279) +++ /trunk/psLib/src/imageops/psImageCovariance.c (revision 21280) @@ -11,4 +11,11 @@ #include "psImageCovariance.h" + +psKernel *psImageCovarianceNone(void) +{ + psKernel *covar = psKernelAlloc(0, 0, 0, 0); // Covariance pseudo-matrix + covar->kernel[0][0] = 1.0; + return covar; +} @@ -25,12 +32,4 @@ // where M^x is the covariance matrix for x. // Note that the errors in f are correlated (covariance) even if the errors in x are not. - // - // We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which - // makes storage difficult; and if that's not enough, the time to do the calculation is definitely - // impractical). Since there are (generally) lots of zeros in the covariance matrix, and the same basic - // pattern repeats (for background pixels), we can just carry that pattern. We carry this in a psKernel, - // since the values are the covariance between the pixel of consideration (at 0,0 in the kernel) and - // neighbouring pixels. Note that this may not be strictly correct near sources, but this is the best we - // can do (and much better than most currently do). psKernel *covar; // Covariance matrix to use @@ -108,5 +107,44 @@ float psImageCovarianceFactor(const psKernel *covariance) { - return covariance ? covariance->kernel[0][0] : 1.0; + return covariance ? covariance->kernel[0][0] : NAN; } +psKernel *psImageCovarianceAverage(const psArray *array) +{ + PS_ASSERT_ARRAY_NON_NULL(array, NULL); + PS_ASSERT_ARRAY_NON_EMPTY(array, NULL); + + int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance + int num = 0; // Number of good matrices to average + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + xMin = PS_MIN(xMin, covar->xMin); + xMax = PS_MAX(xMax, covar->xMax); + yMin = PS_MIN(yMin, covar->yMin); + yMax = PS_MIN(yMax, covar->yMax); + num++; + } + if (num == 0) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging."); + return NULL; + } + + psKernel *average = psKernelAlloc(xMin, xMax, yMin, yMax); // Average covariance + for (int i = 0; i < array->n; i++) { + psKernel *covar = array->data[i]; // Covariance matrix + if (!covar) { + continue; + } + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + average->kernel[y][x] += covar->kernel[y][x]; + } + } + } + psBinaryOp(average->image, average->image, "/", psScalarAlloc(num, PS_TYPE_F32)); + + return average; +} Index: /trunk/psLib/src/imageops/psImageCovariance.h =================================================================== --- /trunk/psLib/src/imageops/psImageCovariance.h (revision 21279) +++ /trunk/psLib/src/imageops/psImageCovariance.h (revision 21280) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-28 22:16:33 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:55:27 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ @@ -17,4 +17,15 @@ #include + +// We don't carry the entire covariance matrix for an image (the size goes as N^2, for N pixels, which makes +// storage difficult; and if that's not enough, the time to do the calculation is definitely impractical). +// Since there are (generally) lots of zeros in the covariance matrix, and the same basic pattern repeats (for +// background pixels), we can just carry that pattern. We carry this in a psKernel, since the values are the +// covariance between the pixel of consideration (at 0,0 in the kernel) and neighbouring pixels. Note that +// this may not be strictly correct near sources, but this is the best we can do (and much better than most +// currently do). + +/// Allocate a covariance pseudo-matrix with no covariance +psKernel *psImageCovarianceNone(void); /// Calculate the covariance pseudo-matrix for a convolution kernel @@ -29,4 +40,9 @@ ); +/// Average many covariance pseudo-matrices +psKernel *psImageCovarianceAverage( + const psArray *array ///< Array of covariance pseudo-matrices + ); + /// @} Index: /trunk/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.c (revision 21279) +++ /trunk/psLib/src/imageops/psImageInterpolate.c (revision 21280) @@ -7,6 +7,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:37 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:55:27 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -866,6 +866,7 @@ -psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, psImageMaskType *maskValue, - float x, float y, const psImageInterpolation *interp) +psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) { PS_ASSERT_PTR_NON_NULL(imageValue, PS_INTERPOLATE_STATUS_ERROR); @@ -963,4 +964,29 @@ +psKernel *psImageInterpolationKernel(float x, float y, psImageInterpolateMode mode) +{ + int size = kernelSizes[mode]; // Size of kernel + + // Kernel basics + INTERPOLATE_SETUP(x, y); + xExact = yExact = false; + + psF32 xKernel[size], yKernel[size]; // Interpolation kernels + interpolationKernel(xKernel, xFrac, mode); + interpolationKernel(yKernel, yFrac, mode); + + int min = -size/2, max = (size - 1) / 2; // Range for kernel + psKernel *kernel = psKernelAlloc(min, max, min, max); // Kernel to return + + for (int y = 0; y < size; y++) { + for (int x = 0; x < size; x++) { + kernel->kernel[y][x] = yKernel[y] * xKernel[x]; + } + } + + return kernel; +} + + psImageInterpolateMode psImageInterpolateModeFromString(const char *name) { Index: /trunk/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /trunk/psLib/src/imageops/psImageInterpolate.h (revision 21279) +++ /trunk/psLib/src/imageops/psImageInterpolate.h (revision 21280) @@ -7,6 +7,6 @@ * @author Paul Price, Institute for Astronomy * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:37 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-04 02:55:27 $ * Copyright 2004-2007 Institute for Astronomy, University of Hawaii */ @@ -50,9 +50,9 @@ const psImage *variance; ///< Variance image for interpolation const psImage *mask; ///< Mask image for interpolation - psImageMaskType maskVal; ///< Value to mask + psImageMaskType maskVal; ///< Value to mask double badImage; ///< Image value if x,y location is not good double badVariance; ///< Variance value if x,y location is not good - psImageMaskType badMask; ///< Mask value to give bad pixels - psImageMaskType poorMask; ///< Mask value to give poor pixels + psImageMaskType badMask; ///< Mask value to give bad pixels + psImageMaskType poorMask; ///< Mask value to give poor pixels float poorFrac; ///< Fraction of flux in bad pixels before output is marked bad bool shifting; ///< Shifting images? Don't interpolate if the shift is exact. @@ -101,3 +101,8 @@ ); +/// Generate the appropriate interpolation kernel +psKernel *psImageInterpolationKernel(float x, float y, ///< Position of interest + psImageInterpolateMode mode ///< Interpolation mode + ); + #endif