Index: trunk/psLib/src/imageops/Makefile.am =================================================================== --- trunk/psLib/src/imageops/Makefile.am (revision 21206) +++ trunk/psLib/src/imageops/Makefile.am (revision 21207) @@ -7,4 +7,5 @@ psImageBackground.c \ psImageConvolve.c \ + psImageCovariance.c \ psImageGeomManip.c \ psImageInterpolate.c \ @@ -25,4 +26,5 @@ psImageBackground.h \ psImageConvolve.h \ + psImageCovariance.h \ psImageGeomManip.h \ psImageInterpolate.h \ Index: trunk/psLib/src/imageops/psImageCovariance.c =================================================================== --- trunk/psLib/src/imageops/psImageCovariance.c (revision 21207) +++ trunk/psLib/src/imageops/psImageCovariance.c (revision 21207) @@ -0,0 +1,112 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include + +#include "psMemory.h" +#include "psConstants.h" +#include "psImageConvolve.h" +#include "psTrace.h" + +#include "psImageCovariance.h" + + +psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance) +{ + PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); + + // See http://en.wikipedia.org/wiki/Error_propagation + // + // If + // f_k = sum_i A_ik x_i + // is a set of functions, then the covariance matrix for f is given by: + // M^f_ij = sum_k sum_l A_ik M^x_kl A_lj + // 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 + if (covariance) { + covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const + } else { + covar = psKernelAlloc(0, 0, 0, 0); + covar->kernel[0][0] = 1.0; + } + + // The above (double) sum for the covariance matrix means that, for each point in the output covariance + // matrix, we need to work out all combinations of getting to the central point via a kernel, input + // covariance matrix and another kernel. This means that the resultant covariance matrix has twice the + // size of the kernel plus the size of the input covariance matrix. + + // Covariance matrix for output + psKernel *out = psKernelAlloc(2 * kernel->xMin + covar->xMin, 2 * kernel->xMax + covar->xMax, + 2 * kernel->yMin + covar->yMin, 2 * kernel->yMax + covar->yMax); + + // Need to go: + // (x,y) --kernel--> (u,v) --covar--> (p,q) --kernel--> (0,0) + // All coordinates (x,y), (u,v), and (p,q) are "absolute", meaning that they are in x,y space, which is + // the space of the output covariance matrix. + // + // So, the ranges in the coordinates are: + // (x,y): -outSize:outSize + // (u,v): (x,y)-kernelSize:(x,y)+kernelSize AND -kernelSize-inSize:kernelSize+inSize + // (p,q): -kernelSize:kernelSize AND (u,v)-inSize:(u,v)+inSize + // Here outSize is the size of the output covariance matrix, kernelSize is the size of the convolution + // kernel, and inSize is the size of the input covariance matrix. + // Since (u,v) and (p,q) have two ways of specifying the range (one from the target coordinate and one + // from the source coordinate), we take the smallest possible (because everything else is zero outside). + + double total = 0.0; // Total covariance + for (int y = out->yMin; y <= out->yMax; y++) { + // Range for v + int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin); + int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax); + for (int x = out->xMin; x <= out->xMax; x++) { + // Range for u + int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin); + int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax); + + double sum = 0.0; // Sum for value of covariance matrix at (x,y) + for (int v = vMin; v <= vMax; v++) { + // Range for q + int qMin = PS_MAX(v + covar->yMin, kernel->yMin); + int qMax = PS_MIN(v + covar->yMax, kernel->yMax); + for (int u = uMin; u <= uMax; u++) { + // Range for p + int pMin = PS_MAX(u + covar->xMin, kernel->xMin); + int pMax = PS_MIN(u + covar->xMax, kernel->xMax); + + double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v) + + double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0) + for (int q = qMin; q <= qMax; q++) { + for (int p = pMin; p <= pMax; p++) { + uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p]; + } + } + sum += xyuvValue * uvpqValue; + } + } + out->kernel[y][x] = sum; + total += sum; + } + } + psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]); + + psFree(covar); + return out; +} + +float psImageCovarianceFactor(const psKernel *covariance) +{ + return covariance ? covariance->kernel[0][0] : 1.0; +} + Index: trunk/psLib/src/imageops/psImageCovariance.h =================================================================== --- trunk/psLib/src/imageops/psImageCovariance.h (revision 21207) +++ trunk/psLib/src/imageops/psImageCovariance.h (revision 21207) @@ -0,0 +1,33 @@ +/* @file psImageCovariance.h + * + * @brief Calculations involving image covariance + * + * @author Paul Price, IfA + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-01-28 22:16:33 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PS_IMAGE_COVARIANCE_H +#define PS_IMAGE_COVARIANCE_H + +/// @addtogroup ImageOps Image Operations +/// @{ + +#include + +/// Calculate the covariance pseudo-matrix for a convolution kernel +psKernel *psImageCovarianceCalculate( + const psKernel *kernel, ///< Convolution kernel + const psKernel *covariance ///< Current covariance pseudo-matrix + ); + +/// Return the pixel-to-pixel covariance factor +float psImageCovarianceFactor( + const psKernel *covariance ///< Covariance pseudo-matrix + ); + + +/// @} +#endif // #ifndef PS_IMAGE_COVARIANCE_H