Index: trunk/psLib/src/imageops/Makefile.am =================================================================== --- trunk/psLib/src/imageops/Makefile.am (revision 36364) +++ trunk/psLib/src/imageops/Makefile.am (revision 36375) @@ -7,4 +7,5 @@ psImageBackground.c \ psImageConvolve.c \ + psImageConvolve2dCache.c \ psImageCovariance.c \ psImageGeomManip.c \ Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 36364) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 36375) @@ -36,9 +36,6 @@ #define MIN_GAUSS_FRAC 0.25 // Minimum Gaussian fraction to accept when smoothing - - static bool threaded = false; // Run image convolution threaded? static pthread_mutex_t threadMutex = PTHREAD_MUTEX_INITIALIZER; - static void kernelFree(psKernel *kernel) @@ -780,5 +777,5 @@ } -void psImageSmooth_PreAlloc_DataFree (psImageSmooth_PreAlloc_Data *smdata) { +void psImageSmoothCacheDataFree (psImageSmoothCacheData *smdata) { psFree (smdata->resultX); psFree (smdata->resultY); @@ -786,8 +783,11 @@ } -psImageSmooth_PreAlloc_Data *psImageSmooth_PreAlloc_DataAlloc (psImage *image, double sigma, double Nsigma) { - - psImageSmooth_PreAlloc_Data *smdata = psAlloc(sizeof(psImageSmooth_PreAlloc_Data)); - psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth_PreAlloc_DataFree); +// allocate the psImageSmoothCache data structure, but do not define the kernel +psImageSmoothCacheData *psImageSmoothCacheAlloc (psImage *image, double sigma, double Nsigma) { + + psImageSmoothCacheData *smdata = psAlloc(sizeof(psImageSmoothCacheData)); + psMemSetDeallocator(smdata, (psFreeFunc) psImageSmoothCacheDataFree); + + smdata->kernel = NULL; if (!image) { @@ -796,5 +796,4 @@ smdata->Nx = 0; smdata->Ny = 0; - smdata->kernel = NULL; smdata->resultX = NULL; smdata->resultY = NULL; @@ -807,5 +806,6 @@ smdata->Ny = image->numRows; // Number of rows - IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32); + // XXX drop this : we now require a call to a kernel-creation function (like psImageSmoothCacheKernel_Gauss) + // IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32); // use a temp running buffer for X and Y directions. @@ -819,8 +819,18 @@ } +// generate a Gaussian smoothing kernel for supplied sigma. sigma here does not need to match +// that used to allocate the structure, but it is recommended +bool psImageSmoothCacheKernel_Gauss (psImageSmoothCacheData *smdata, float sigma) { + // check for NULL structure elements? + psFree (smdata->kernel); + IMAGE_SMOOTH_GAUSS(smdata->kernel, smdata->Nrange, sigma, F32); + return true; +} + // we can use the same DATA structure on multiple images of the same size -bool psImageSmooth_PreAlloc_F32(psImage *image, psImageSmooth_PreAlloc_Data *smdata) +bool psImageSmoothCache_F32(psImage *image, psImageSmoothCacheData *smdata) { PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_VECTOR_NON_NULL(smdata->kernel, false); // assert on data type @@ -2156,3 +2166,2 @@ return threaded; } - Index: trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.h (revision 36364) +++ trunk/psLib/src/imageops/psImageConvolve.h (revision 36375) @@ -26,5 +26,5 @@ #define PS_TYPE_KERNEL_NAME "psF32" ///< the data type for kernel as a string */ -/// a structure to contain data related to image smoothing with a 1D gauss kernel +/// a structure to contain data related to image smoothing with a pre-cached 1D gauss kernel typedef struct { int Nx; @@ -34,5 +34,12 @@ psF32 *resultY; psVector *kernel; -} psImageSmooth_PreAlloc_Data; +} psImageSmoothCacheData; + +/// a structure to contain data related to image smoothing with a pre-cached 1D gauss kernel +typedef struct { + float Nsigma; + int Ns; // number of pixel radii + float *radflux; // conv kernel in special positions +} psImageSmooth2dCacheData; /// A convolution kernel @@ -306,6 +313,7 @@ ); -psImageSmooth_PreAlloc_Data *psImageSmooth_PreAlloc_DataAlloc (psImage *image, double sigma, double Nsigma); -bool psImageSmooth_PreAlloc_F32(psImage *image, psImageSmooth_PreAlloc_Data *smdata); +psImageSmoothCacheData *psImageSmoothCacheAlloc (psImage *image, double sigma, double Nsigma); +bool psImageSmoothCache_F32(psImage *image, psImageSmoothCacheData *smdata); +bool psImageSmoothCacheKernel_Gauss (psImageSmoothCacheData *smdata, float sigma); /// Control threading for image convolution functions @@ -318,4 +326,9 @@ bool psImageConvolveGetThreads(void); +psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma); +bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa); +bool psImageSmooth2dCacheKernel_Gauss (psImageSmooth2dCacheData *smdata, float sigma); +bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata); + /// @} #endif // #ifndef PS_IMAGE_CONVOLVE_H Index: trunk/psLib/src/imageops/psImageConvolve2dCache.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve2dCache.c (revision 36375) +++ trunk/psLib/src/imageops/psImageConvolve2dCache.c (revision 36375) @@ -0,0 +1,405 @@ +/// @file psImageConvolve2dCache.c -- specialized 2D convolution for psf matching +/// @author Eugene Magnier, IfA +/// +/// @date $Date: 2009-02-05 23:56:14 $ +/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include "psAbort.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psError.h" +#include "psAssert.h" +#include "psScalar.h" +#include "psBinaryOp.h" +#include "psImageFFT.h" +#include "psImageStructManip.h" +#include "psImagePixelManip.h" +#include "psTrace.h" +#include "psThread.h" + +#include "psImageConvolve.h" + +// fixed set of radii and number of entries in circularly symmetric profile +# define NRAD_MAX 87 +# define RAD_MAX 15 +static float radii2[NRAD_MAX] = { 0.0, 1.0, 2.0, 4.0, 5.0, + 8.0, 9.0, 10.0, 13.0, 16.0, + 17.0, 18.0, 20.0, 25.0, 26.0, + + 29.0, 32.0, 34.0, 36.0, 37.0, + 40.0, 41.0, 45.0, 49.0, 50.0, + 52.0, 53.0, 58.0, 61.0, 64.0, + + 65.0, 68.0, 72.0, 73.0, 74.0, + 80.0, 81.0, 82.0, 85.0, 89.0, + 90.0, 97.0, 98.0, 100.0, 101.0, + + 104.0, 106.0, 109.0, 113.0, 116.0, + 117.0, 121.0, 122.0, 125.0, 128.0, + 130.0, 136.0, 137.0, 144.0, 145.0, + + 146.0, 148.0, 149.0, 153.0, 157.0, + 160.0, 162.0, 164.0, 169.0, 170.0, + 173.0, 178.0, 180.0, 181.0, 185.0, + + 193.0, 194.0, 196.0, 197.0, 200.0, + 202.0, 205.0, 208.0, 212.0, 218.0, + 221.0, 225.0 +}; + +static int radiiN[NRAD_MAX] = { 1, 4, 4, 4, 8, + 4, 4, 8, 8, 4, + 8, 4, 8, 12, 8, + + 8, 4, 8, 4, 8, + 8, 8, 8, 4, 12, + 8, 8, 8, 8, 4, + + 16, 8, 4, 8, 8, + 8, 4, 8, 16, 8, + 8, 8, 4, 12, 8, + + 8, 8, 8, 8, 8, + 8, 4, 8, 16, 4, + 16, 8, 8, 8, 16, + + 8, 8, 8, 8, 8, + 8, 4, 8, 12, 16, + 8, 8, 8, 8, 16, + + 8, 8, 4, 8, 12, + 8, 16, 8, 8, 8, + 16, 12 +}; + +# define ADD_AXIS(RAD,DD) s += radflux[RAD]*(vi[iy - DD][ix] + \ + vi[iy + DD][ix] + \ + vi[iy][ix - DD] + \ + vi[iy][ix + DD]); + +# define ADD_DIAG(RAD,DD) s += radflux[RAD]*(vi[iy - DD][ix - DD] + \ + vi[iy - DD][ix + DD] + \ + vi[iy + DD][ix - DD] + \ + vi[iy + DD][ix + DD]); + +# define ADD_RAND(RAD,DX,DY) s += radflux[RAD]*(vi[iy - DY][ix - DX] + \ + vi[iy - DY][ix + DX] + \ + vi[iy + DY][ix - DX] + \ + vi[iy + DY][ix + DX] + \ + vi[iy - DX][ix - DY] + \ + vi[iy - DX][ix + DY] + \ + vi[iy + DX][ix - DY] + \ + vi[iy + DX][ix + DY]); + +# if (0) +# define ADD_AXIS(RAD,DD) s += radflux[RAD]*(vi[p - DD] + vi[p + DD] + vi[p - DD*Nx] + vi[p + DD*Nx]); +# define ADD_DIAG(RAD,DD) s += radflux[RAD]*(vi[p - DD - DD*Nx] + vi[p + DD - DD*Nx] + vi[p - DD + DD*Nx] + vi[p + DD + DD*Nx]); +# define ADD_RAND(RAD,DX,DY) s += radflux[RAD]*(vi[p - DX - DY*Nx] + vi[p + DX - DY*Nx] + vi[p - DX + DY*Nx] + vi[p + DX + DY*Nx] + + vi[p - DY - DX*Nx] + vi[p + DY - DX*Nx] + vi[p - DY + DX*Nx] + vi[p + DY + DX*Nx]); +# endif + +void psImageSmooth2dCacheDataFree (psImageSmooth2dCacheData *smdata) { + + if (smdata->radflux == NULL) return; + psFree (smdata->radflux); +} + +// allocate the psImageSmooth2dCache data structure, but do not define the kernel +psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (float Nsigma) { + + psAssert (isfinite(Nsigma) && (Nsigma > 0.5), "Nsigma is not valid"); + + psImageSmooth2dCacheData *smdata = psAlloc(sizeof(psImageSmooth2dCacheData)); + psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth2dCacheDataFree); + + smdata->radflux = NULL; + smdata->Nsigma = Nsigma; + smdata->Ns = -1; + + return smdata; +} + +// generate a 2D smoothing kernel for supplied sigma & kappa (PS1_V1 profile). +bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa) { + + // check for NULL structure elements? + int Ns = (int)(smdata->Nsigma * sigma); + Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX)); + smdata->Ns = Ns; + + int Ns2 = Ns * Ns; + + // we are going to use a hard-wired set of radial points + smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX); + + double sum = 0.0; + for (int i = 0; i < NRAD_MAX; i++) { + if (radii2[i] > Ns2) { + smdata->radflux[i] = 0.0; + continue; + } + float z = 0.5 * radii2[i] / PS_SQR(sigma); + smdata->radflux[i] = 1.0 / (1.0 + kappa*z + pow(z,1.666)); + sum += radiiN[i] * smdata->radflux[i]; + } + for (int i = 0; i < NRAD_MAX; i++) { + smdata->radflux[i] = smdata->radflux[i] / sum; + } + + return true; +} + +// generate a 2D smoothing kernel for supplied sigma & kappa (PS1_V1 profile). +bool psImageSmooth2dCacheKernel_Gauss (psImageSmooth2dCacheData *smdata, float sigma) { + + // check for NULL structure elements? + int Ns = (int)(smdata->Nsigma * sigma); + Ns = PS_MAX (3, PS_MIN (Ns, RAD_MAX)); + smdata->Ns = Ns; + + int Ns2 = Ns * Ns; + + // we are going to use a hard-wired set of radial points + smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX); + + float sum = 0.0; + for (int i = 0; i < NRAD_MAX; i++) { + if (radii2[i] > Ns2) { + smdata->radflux[i] = 0.0; + continue; + } + float z = 0.5 * radii2[i] / PS_SQR(sigma); + smdata->radflux[i] = exp(-z); + sum += radiiN[i] * smdata->radflux[i]; + } + for (int i = 0; i < NRAD_MAX; i++) { + smdata->radflux[i] = smdata->radflux[i] / sum; + } + + return true; +} + +// we can use the same DATA structure on multiple images of the same size +bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata) +{ + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_PTR_NON_NULL(smdata->radflux, false); + // assert on data type + + // relevant terms + int Ns = smdata->Ns; // Number of pixels either side for convolution kernel + + int Nx = image->numCols; + int Ny = image->numRows; + + int Nxtmp = Nx + 2*Ns; + int Nytmp = Ny + 2*Ns; + + // copy input image into a buffer padded by Ns on either side + float **vi = (float **) psAlloc(sizeof(float *)*Nytmp); + for (int iy = 0; iy < Nytmp; iy ++) { + int Iy = iy - Ns; + vi[iy] = (float *) psAlloc(sizeof(float)*Nxtmp); + memset (vi[iy], 0, sizeof(float)*Nxtmp); + if (Iy < 0) continue; + if (Iy >= image->numRows) continue; + for (int Ix = 0; Ix < image->numCols; Ix ++) { + int ix = Ix + Ns; + vi[iy][ix] = image->data.F32[Iy][Ix]; + } + } + + // set up the output buffer (different from input) + float **vo = (float **) psAlloc(sizeof(float *)*Nytmp); + for (int iy = 0; iy < Nytmp; iy ++) { + vo[iy] = (float *) psAlloc(sizeof(float)*Nxtmp); + memset (vo[iy], 0, sizeof(float)*Nxtmp); + } + + float *radflux = smdata->radflux; + + // smooth in 2D (ix,iy is the coordinate in the input and output images, only need to + // transform the region inside the padding. + for (int iy = Ns; iy < Ny + Ns; iy ++) { + for (int ix = Ns; ix < Nx + Ns; ix ++) { + + float s = radflux[0] * vi[iy][ix]; + + // r <= 1.0 + ADD_AXIS (1, 1); // r^2 = 1 + + // r <= 2.0 + ADD_DIAG (2, 1); // r^2 = 2 + ADD_AXIS (3, 2); // r^2 = 4 + + // r <= 3.0 + ADD_RAND (4, 1, 2); // r^2 = 5 + ADD_DIAG (5, 2); // r^2 = 8 + ADD_AXIS (6, 3); // r^2 = 9 + if (Ns <= 3) goto finish; + + // r <= 4.0 + ADD_RAND (7, 1, 3); // r^2 = 10 + ADD_RAND (8, 2, 3); // r^2 = 13 + ADD_AXIS (9, 4); // r^2 = 16 + if (Ns <= 4) goto finish; + + // r <= 5.0 + ADD_RAND (10, 1, 4); // r^2 = 17 + ADD_DIAG (11, 3); // r^2 = 18 + ADD_RAND (12, 2, 4); // r^2 = 20 + ADD_RAND (13, 3, 4); // r^2 = 25 + ADD_AXIS (13, 5); // r^2 = 25 + if (Ns <= 5) goto finish; + + // r <= 6.0 + ADD_RAND (14, 1, 5); // r^2 = 26 + ADD_RAND (15, 2, 5); // r^2 = 29 + ADD_DIAG (16, 4); // r^2 = 32 + ADD_RAND (17, 3, 5); // r^2 = 34 + ADD_AXIS (18, 6); // r^2 = 36 + if (Ns <= 6) goto finish; + + // r <= 7.0 + ADD_RAND (19, 1, 6); // r^2 = 37 + ADD_RAND (20, 2, 6); // r^2 = 40 + ADD_RAND (21, 4, 5); // r^2 = 41 + ADD_RAND (22, 3, 6); // r^2 = 45 + ADD_AXIS (23, 7); // r^2 = 49 + if (Ns <= 7) goto finish; + + // r <= 8.0 + ADD_RAND (24, 1, 7); // r^2 = 50 + ADD_DIAG (24, 5); // r^2 = 50 + ADD_RAND (25, 4, 6); // r^2 = 52 + ADD_RAND (26, 2, 7); // r^2 = 53 + ADD_RAND (27, 3, 7); // r^2 = 58 + ADD_RAND (28, 5, 6); // r^2 = 61 + ADD_AXIS (29, 8); // r^2 = 64 + if (Ns <= 8) goto finish; + + ADD_RAND (30, 1, 8); // r^2 = 65 * + ADD_RAND (30, 4, 7); // r^2 = 65 * + ADD_RAND (31, 2, 8); // r^2 = 68 + ADD_DIAG (32, 6); // r^2 = 72 + ADD_RAND (33, 3, 8); // r^2 = 73 + ADD_RAND (34, 5, 7); // r^2 = 74 + ADD_RAND (35, 4, 8); // r^2 = 80 + ADD_AXIS (36, 9); // r^2 = 81 + if (Ns <= 9) goto finish; + + ADD_RAND (37, 1, 9); // r^2 = 82 + ADD_RAND (38, 2, 9); // r^2 = 85 * + ADD_RAND (38, 6, 7); // r^2 = 85 * + ADD_RAND (39, 5, 8); // r^2 = 89 + ADD_RAND (40, 3, 9); // r^2 = 90 + ADD_RAND (41, 4, 9); // r^2 = 97 + ADD_DIAG (42, 7); // r^2 = 98 + ADD_RAND (43, 6, 8); // r^2 = 100 * + ADD_AXIS (43, 10); // r^2 = 100 * + if (Ns <= 10) goto finish; + + ADD_RAND ( 44, 1, 10); // r^2 = 101 + ADD_RAND ( 45, 2, 10); // r^2 = 104 + ADD_RAND ( 46, 5, 9); // r^2 = 106 + ADD_RAND ( 47, 3, 10); // r^2 = 109 + ADD_RAND ( 48, 7, 8); // r^2 = 113 + ADD_RAND ( 49, 4, 10); // r^2 = 116 + ADD_RAND ( 50, 6, 9); // r^2 = 117 + ADD_AXIS ( 51, 11) ; // r^2 = 121 + if (Ns <= 11) goto finish; + + ADD_RAND ( 52, 1, 11); // r^2 = 122 + ADD_RAND ( 53, 2, 11); // r^2 = 125 * + ADD_RAND ( 53, 5, 10); // r^2 = 125 * + ADD_DIAG ( 54, 8) ; // r^2 = 128 + ADD_RAND ( 55, 3, 11); // r^2 = 130 * + ADD_RAND ( 55, 7, 9); // r^2 = 130 * + ADD_RAND ( 56, 6, 10); // r^2 = 136 + ADD_RAND ( 57, 4, 11); // r^2 = 137 + ADD_AXIS ( 58, 12) ; // r^2 = 144 + if (Ns <= 12) goto finish; + + ADD_RAND ( 59, 1, 12); // r^2 = 145 * + ADD_RAND ( 59, 8, 9); // r^2 = 145 * + ADD_RAND ( 60, 5, 11); // r^2 = 146 + ADD_RAND ( 61, 2, 12); // r^2 = 148 + ADD_RAND ( 62, 7, 10); // r^2 = 149 + ADD_RAND ( 63, 3, 12); // r^2 = 153 + ADD_RAND ( 64, 6, 11); // r^2 = 157 + ADD_RAND ( 65, 4, 12); // r^2 = 160 + ADD_DIAG ( 66, 9) ; // r^2 = 162 + ADD_RAND ( 67, 8, 10); // r^2 = 164 + ADD_RAND ( 68, 5, 12); // r^2 = 169 * + ADD_AXIS ( 68, 13) ; // r^2 = 169 * + if (Ns <= 13) goto finish; + + ADD_RAND ( 69, 1, 13); // r^2 = 170 * + ADD_RAND ( 69, 7, 11); // r^2 = 170 * + ADD_RAND ( 70, 2, 13); // r^2 = 173 + ADD_RAND ( 71, 3, 13); // r^2 = 178 + ADD_RAND ( 72, 6, 12); // r^2 = 180 + ADD_RAND ( 73, 9, 10); // r^2 = 181 + ADD_RAND ( 74, 4, 13); // r^2 = 185 * + ADD_RAND ( 74, 8, 11); // r^2 = 185 * + ADD_RAND ( 75, 7, 12); // r^2 = 193 + ADD_RAND ( 76, 5, 13); // r^2 = 194 + ADD_AXIS ( 77, 14) ; // r^2 = 196 + if (Ns <= 14) goto finish; + + ADD_RAND ( 78, 1, 14); // r^2 = 197 + ADD_DIAG ( 79, 10) ; // r^2 = 200 * + ADD_RAND ( 79, 2, 14); // r^2 = 200 * + ADD_RAND ( 80, 9, 11); // r^2 = 202 + ADD_RAND ( 81, 3, 14); // r^2 = 205 * + ADD_RAND ( 81, 6, 13); // r^2 = 205 * + ADD_RAND ( 82, 8, 12); // r^2 = 208 + ADD_RAND ( 83, 4, 14); // r^2 = 212 + ADD_RAND ( 84, 7, 13); // r^2 = 218 + ADD_RAND ( 85, 10, 11); // r^2 = 221 * + ADD_RAND ( 85, 5, 14); // r^2 = 221 * + ADD_RAND ( 86, 9, 12); // r^2 = 225 * + ADD_AXIS ( 86, 15); // r^2 = 225 * + if (Ns <= 15) goto finish; + + finish: + vo[iy][ix] = s; + } + } + + for (int iy = 0; iy < Ny; iy ++) { + int Iy = iy + Ns; + for (int ix = 0; ix < Nx; ix ++) { + int Ix = ix + Ns; + image->data.F32[iy][ix] = vo[Iy][Ix]; + } + } + + for (int iy = 0; iy < Nytmp; iy ++) { + psFree (vi[iy]); + psFree (vo[iy]); + } + psFree (vi); + psFree (vo); + + return true; +} + +/* +101, 104, 106, 109, 113, +116, 117, 121, 122, 125, +128, 130, 136, 137, 144, +145, 146, 148, 149, 153, +157, 160, 162, 164, 169, +170, 173, 178, 180, 181, +185, 193, 194, 196, 197, +200, 202, 205, 208, 212, +218, 221, 225 +*/ + Index: trunk/psLib/src/math/psEllipse.c =================================================================== --- trunk/psLib/src/math/psEllipse.c (revision 36364) +++ trunk/psLib/src/math/psEllipse.c (revision 36375) @@ -12,4 +12,7 @@ // f = exp(-z) where z describes the elliptical contour at any flux level: +// NOTE: major, minor are the 1-sigma lengths in the major,minor directions assuming the +// moments represent a Gaussian profile + // sigma shape: z = 0.5((x/sx)^2 + (y/sy)^2 + sxy*x*y) // sigma axes : z = 0.5((x/sa)^2 + (y/sb)^2), x,y rotated by theta @@ -17,5 +20,123 @@ // polarization : e0, e1, e2 -// ellipse rotation (major, minor, theta) -> (x2, y2, xy) +// ellipse rotation (major, minor, theta) -> (Mxx, Mxy, Myy) +psEllipseMoments psEllipseAxesToMoments(psEllipseAxes axes) +{ + psEllipseMoments moments; + + double f1 = PS_SQR(axes.major) + PS_SQR(axes.minor); + double f2 = PS_SQR(axes.major) - PS_SQR(axes.minor); + + moments.x2 = +0.5*f1 + 0.5*f2*cos(2*axes.theta); + moments.y2 = +0.5*f1 - 0.5*f2*cos(2*axes.theta); + moments.xy = +0.5*f2*sin(2*axes.theta); + + assert (isfinite(moments.x2)); + assert (isfinite(moments.y2)); + assert (isfinite(moments.xy)); + + return moments; +} + +// ellipse rotation (Mxx, Mxy, Myy) -> (major, minor, theta). +psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments, double maxAR) +{ + psEllipseAxes axes; + psEllipseAxes badValue = {NAN, NAN, NAN}; + + if (!isfinite(moments.x2)) return badValue; + if (!isfinite(moments.y2)) return badValue; + if (!isfinite(moments.xy)) return badValue; + + if (moments.x2 < 0) return badValue; + if (moments.y2 < 0) return badValue; + + double g1 = moments.x2 + moments.y2; + double g2 = moments.x2 - moments.y2; + double g3 = sqrt(PS_SQR(g2) + 4*PS_SQR(moments.xy)); + + axes.major = sqrt (0.5*(g1 + g3)); + axes.theta = +0.5 * atan2 (+2.0*moments.xy, g2); // theta in radians + + // long, thin objects are likely to have a poorly measured minor axis + // the angle and major axis are likely to be ok. + // restrict the axis ratio + double rAR2 = (g1 - g3) / (g1 + g3); + if (rAR2 < 1.0/PS_SQR(maxAR)) { + axes.minor = axes.major / maxAR; + } else { + axes.minor = sqrt (0.5*(g1 - g3)); + } + + assert (isfinite(axes.major)); + assert (isfinite(axes.minor)); + assert (isfinite(axes.theta)); + + return axes; +} + +// ellipse rotation (major, minor, theta) -> (sx, sy, sxy) +// theta is postive rotation of major axis away from x-axis +psEllipseShape psEllipseAxesToShape(psEllipseAxes axes) +{ + psEllipseShape shape; + psEllipseShape badValue = {NAN, NAN, NAN}; + + if (!isfinite(axes.minor)) return badValue; + if (!isfinite(axes.major)) return badValue; + if (!isfinite(axes.theta)) return badValue; + + if (axes.minor <= 0) return badValue; + if (axes.major <= 0) return badValue; + + double f1 = 1.0 / PS_SQR(axes.minor) + 1.0 / PS_SQR(axes.major); + double f2 = 1.0 / PS_SQR(axes.minor) - 1.0 / PS_SQR(axes.major); + + double sxr = 0.5*f1 - 0.5*f2*cos(2*axes.theta); + double syr = 0.5*f1 + 0.5*f2*cos(2*axes.theta); + + // sxr, syr cannot be < 0 (f1 >= f2) + + shape.sx = +1.0 / sqrt(sxr); + shape.sy = +1.0 / sqrt(syr); + shape.sxy = -0.5*f2*sin(2*axes.theta); + + assert (isfinite(shape.sx)); + assert (isfinite(shape.sy)); + assert (isfinite(shape.sxy)); + + return shape; +} + +// ellipse derotation (sx, sy, sxy) -> (major, minor, theta) +psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape, double maxAR) +{ + psEllipseAxes axes; + + double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx); + double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx); + double f3 = sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); + + axes.minor = sqrt (2.0 / (f1 + f3)); + axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2); + + // long, thin objects are likely to have a poorly measured major axis + // the angle and minor axis are likely to be ok. + // restrict the axis ratio + double rAR2 = (f1 - f3) / (f1 + f3); + if (rAR2 < 1.0/PS_SQR(maxAR)) { + axes.major = axes.minor * maxAR; + } else { + axes.major = sqrt (2.0 / (f1 - f3)); + } + + assert (isfinite(axes.theta)); + assert (isfinite(axes.major)); + assert (isfinite(axes.minor)); + + return axes; +} + +// ellipse rotation (major, minor, theta) -> (e0, e1, e2) psEllipsePol psEllipseAxesToPol(psEllipseAxes axes) { @@ -35,25 +156,5 @@ } -// ellipse rotation (major, minor, theta) -> (x2, y2, xy) -psEllipsePol psEllipseShapeToPol(psEllipseShape shape) -{ - psEllipsePol pol; - - double r = 1.0 / (1.0 - PS_SQR(shape.sx)*PS_SQR(shape.sy)*PS_SQR(shape.sxy)); - - pol.e0 = r*(PS_SQR(shape.sx) + PS_SQR(shape.sy)); - pol.e1 = r*(PS_SQR(shape.sx) - PS_SQR(shape.sy)); - // XXX I do not understand this negative sign - pol.e2 = -r*(2.0*PS_SQR(shape.sx)*PS_SQR(shape.sy)*shape.sxy); - - assert (isfinite(pol.e0)); - assert (isfinite(pol.e1)); - assert (isfinite(pol.e2)); - - return pol; -} - -// ellipse rotation (major, minor, theta) -> (x2, y2, xy) -// XXXX handle case where e0 < LIMIT +// ellipse rotation (e0, e1, e2) -> (major, minor, theta) psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol, const float minMinorAxis) @@ -96,87 +197,41 @@ } -// ellipse rotation (major, minor, theta) -> (x2, y2, xy) -psEllipseMoments psEllipseAxesToMoments(psEllipseAxes axes) -{ - psEllipseMoments moments; - - double f1 = PS_SQR(axes.major) + PS_SQR(axes.minor); - double f2 = PS_SQR(axes.major) - PS_SQR(axes.minor); - - moments.x2 = +0.5*f1 + 0.5*f2*cos(2*axes.theta); - moments.y2 = +0.5*f1 - 0.5*f2*cos(2*axes.theta); - moments.xy = +0.5*f2*sin(2*axes.theta); - - assert (isfinite(moments.x2)); - assert (isfinite(moments.y2)); - assert (isfinite(moments.xy)); - - return moments; -} - -// ellipse rotation (x2, y2, xy) -> (major, minor, theta). NOTE: major, minor are the -// 1-sigma lengths in the major,minor directions assuming the moments represent a Gaussian -// profile -psEllipseAxes psEllipseMomentsToAxes(psEllipseMoments moments, double maxAR) -{ - psEllipseAxes axes; - psEllipseAxes badValue = {NAN, NAN, NAN}; - - if (!isfinite(moments.x2)) return badValue; - if (!isfinite(moments.y2)) return badValue; - if (!isfinite(moments.xy)) return badValue; - - if (moments.x2 < 0) return badValue; - if (moments.y2 < 0) return badValue; - - double g1 = moments.x2 + moments.y2; - double g2 = moments.x2 - moments.y2; - double g3 = sqrt(PS_SQR(g2) + 4*PS_SQR(moments.xy)); - - axes.major = sqrt (0.5*(g1 + g3)); - axes.theta = +0.5 * atan2 (+2.0*moments.xy, g2); // theta in radians - - // long, thin objects are likely to have a poorly measured minor axis - // the angle and major axis are likely to be ok. - // restrict the axis ratio - double rAR2 = (g1 - g3) / (g1 + g3); - if (rAR2 < 1.0/PS_SQR(maxAR)) { - axes.minor = axes.major / maxAR; - } else { - axes.minor = sqrt (0.5*(g1 - g3)); - } - - assert (isfinite(axes.major)); - assert (isfinite(axes.minor)); - assert (isfinite(axes.theta)); - - return axes; -} - -// ellipse rotation (major, minor, theta) -> (sx, sy, sxy) -// theta is postive rotation of major axis away from x-axis -psEllipseShape psEllipseAxesToShape(psEllipseAxes axes) +// ellipse rotation (sx, sy, sxy) -> (e0, e1, e2) +psEllipsePol psEllipseShapeToPol(psEllipseShape shape) +{ + psEllipsePol pol; + + double r = 1.0 / (1.0 - PS_SQR(shape.sx)*PS_SQR(shape.sy)*PS_SQR(shape.sxy)); + + pol.e0 = r*(PS_SQR(shape.sx) + PS_SQR(shape.sy)); + pol.e1 = r*(PS_SQR(shape.sx) - PS_SQR(shape.sy)); + // XXX I do not understand this negative sign + pol.e2 = -r*(2.0*PS_SQR(shape.sx)*PS_SQR(shape.sy)*shape.sxy); + + assert (isfinite(pol.e0)); + assert (isfinite(pol.e1)); + assert (isfinite(pol.e2)); + + return pol; +} + +// ellipse rotation (e0, e1, e2) -> (sx, sy, sxy) +psEllipseShape psEllipsePolToShape(psEllipsePol pol) { psEllipseShape shape; - psEllipseShape badValue = {NAN, NAN, NAN}; - - if (!isfinite(axes.minor)) return badValue; - if (!isfinite(axes.major)) return badValue; - if (!isfinite(axes.theta)) return badValue; - - if (axes.minor <= 0) return badValue; - if (axes.major <= 0) return badValue; - - double f1 = 1.0 / PS_SQR(axes.minor) + 1.0 / PS_SQR(axes.major); - double f2 = 1.0 / PS_SQR(axes.minor) - 1.0 / PS_SQR(axes.major); - - double sxr = 0.5*f1 - 0.5*f2*cos(2*axes.theta); - double syr = 0.5*f1 + 0.5*f2*cos(2*axes.theta); - - // sxr, syr cannot be < 0 (f1 >= f2) - - shape.sx = +1.0 / sqrt(sxr); - shape.sy = +1.0 / sqrt(syr); - shape.sxy = -0.5*f2*sin(2*axes.theta); + + double q = sqrt (PS_SQR(pol.e1) + PS_SQR(pol.e2)); + double p = PS_SQR(pol.e0) + PS_SQR(q) - 2.0*q*pol.e0; + + // double f1 = 4*pol.e0 / p; + // double f2 = 4*q / p; + + double sxr = 2.0*(pol.e0 - pol.e1) / p; + double syr = 2.0*(pol.e0 + pol.e1) / p; + + shape.sx = sqrt(1.0 / sxr); + shape.sy = sqrt(1.0 / syr); + + shape.sxy = -2.0 * pol.e2 / p; assert (isfinite(shape.sx)); @@ -187,33 +242,4 @@ } -// ellipse derotation (sx, sy, sxy) -> (major, minor, theta) -psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape, double maxAR) -{ - psEllipseAxes axes; - - double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx); - double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx); - double f3 = sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy)); - - axes.minor = sqrt (2.0 / (f1 + f3)); - axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2); - - // long, thin objects are likely to have a poorly measured major axis - // the angle and minor axis are likely to be ok. - // restrict the axis ratio - double rAR2 = (f1 - f3) / (f1 + f3); - if (rAR2 < 1.0/PS_SQR(maxAR)) { - axes.major = axes.minor * maxAR; - } else { - axes.major = sqrt (2.0 / (f1 - f3)); - } - - assert (isfinite(axes.theta)); - assert (isfinite(axes.major)); - assert (isfinite(axes.minor)); - - return axes; -} - // XXX keep this construction? // force the axis ratio to be less than 10 Index: trunk/psLib/src/math/psEllipse.h =================================================================== --- trunk/psLib/src/math/psEllipse.h (revision 36364) +++ trunk/psLib/src/math/psEllipse.h (revision 36375) @@ -38,5 +38,5 @@ } psEllipseShape; -/// Ellipse defined in terms of polarisations +/// Ellipse defined in terms of polarizations typedef struct { double e0; ///< Scale (Mxx + Myy) @@ -65,18 +65,20 @@ ); -/// Convert axes to polarisation representation +/// Convert axes to polarization representation psEllipsePol psEllipseAxesToPol(psEllipseAxes axes ///< Axes of ellipse ); -/// Convert shape to polarisation representation +/// Convert polarization to axes representation +psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol, ///< Polarization of ellipse + const float minMinorAxis ///< Minimum allowed minor axis + ); + +/// Convert shape to polarization representation psEllipsePol psEllipseShapeToPol(psEllipseShape shape ///< Shape of ellipse ); -/// Convert polarisation to axes representation -/// -/// XXX This API goes against the PS convention of outputs being first. -psEllipseAxes psEllipsePolToAxes(const psEllipsePol pol, ///< Polarisation of ellipse - const float minMinorAxis ///< Minimum allowed minor axis - ); +/// Convert shape to polarization representation +psEllipseShape psEllipsePolToShape(psEllipsePol pol ///< Shape of ellipse + ); /// @} Index: trunk/psLib/test/imageops/Makefile.am =================================================================== --- trunk/psLib/test/imageops/Makefile.am (revision 36364) +++ trunk/psLib/test/imageops/Makefile.am (revision 36375) @@ -22,4 +22,5 @@ tap_psImageConvolve \ tap_psImageConvolve2 \ + tap_psImageConvolve2dCache \ tap_psImagePixelExtract \ tap_psImageInterpolate2 \ Index: trunk/psLib/test/imageops/tap_psImageConvolve2dCache.c =================================================================== --- trunk/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36375) +++ trunk/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36375) @@ -0,0 +1,140 @@ +/** @file tap_psImageConvolve2dCache.c + * + * @brief Contains the tests for psImageConvolve2dCache.[ch] + * @author Eugene Magnier, IfA + * @date $Date: 2007-05-01 00:08:52 $ + * + * Copyright 2013 IfA, University of Hawaii + */ + +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa); + +int main(int argc, char **argv) { + + plan_tests(141); + + psMemId id = psMemGetId(); + + // generate an image to convolve + tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 3.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 1.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 100, 200, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 200, 100, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 200, 200, 3.0, 2.0, 2.0); + tap_convolve_delta (301, 301, 100, 100, 3.0, 2.0, 2.0); + + tap_convolve_delta (301, 301, 150, 150, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 4.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 4.0, 3.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 5.0, 3.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 5.0, 3.0, 1.0); + + tap_convolve_delta (301, 301, 150, 150, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 200, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 100, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 200, 200, 5.0, 1.0, 1.0); + tap_convolve_delta (301, 301, 100, 100, 5.0, 1.0, 1.0); + + ok(!psMemCheckLeaks (id, NULL, stdout, false), "no memory leaks"); + + return exit_status(); +} + +bool tap_convolve_delta (int Nx, int Ny, int Xo, int Yo, float Nsigma, float sigma, float kappa) { + + bool status; + + psImage *input = psImageAlloc (Nx, Ny, PS_TYPE_F32); + psImageInit (input, 0.0); + input->data.F32[Yo][Xo] = 1.0; + + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (Nsigma); + ok(smdata, "init smdata"); + ok(!smdata->radflux, "init radflux to NULL"); + ok(smdata->Nsigma == Nsigma, "init Nsigma"); + + status = psImageSmooth2dCacheKernel_PS1_V1 (smdata, sigma, kappa); + ok(status, "generate a kernel"); + + status = psImageSmooth2dCache_F32 (input, smdata); + ok(status, "smoothed image"); + + double s = 0.0; + + psImage *model = psImageAlloc (Nx, Ny, PS_TYPE_F32); + psImageInit (model, 0.0); + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { + float z = 0.5*PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + 0.5*PS_SQR((iy + 0.5 - Yo - 0.5) / sigma); + if (z > 0.5*PS_SQR(smdata->Nsigma)) continue; + + float f = 1.0 / (1.0 + z*kappa + pow(z, 1.666)); + + model->data.F32[iy][ix] = f; + s += f; + } + } + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { + model->data.F32[iy][ix] /= s; + } + } + + float min_df = +10.0; + float max_df = -10.0; + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { + float df = input->data.F32[iy][ix] - model->data.F32[iy][ix]; + min_df = PS_MIN (min_df, df); + max_df = PS_MAX (max_df, df); + } + } + ok (fabs(min_df) < 3e-5, "minor deviations"); + ok (fabs(max_df) < 3e-5, "minor deviations"); + + if ((fabs(min_df) >= 3e-5) || (fabs(max_df) >= 3e-5)) { + psFits *fits = NULL; + fits = psFitsOpen ("input.fits", "w"); + psFitsWriteImage (fits, NULL, input, 0, NULL); + psFitsClose (fits); + + fits = psFitsOpen ("model.fits", "w"); + psFitsWriteImage (fits, NULL, model, 0, NULL); + psFitsClose (fits); + } + + psFree (input); + psFree (model); + psFree (smdata); + + return true; +}