Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/imageops/Makefile.am (revision 36277) +++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/Makefile.am (revision 36278) @@ -7,4 +7,5 @@ psImageBackground.c \ psImageConvolve.c \ + psImageConvolve2dCache.c \ psImageCovariance.c \ psImageGeomManip.c \ Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h (revision 36277) +++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h (revision 36278) @@ -326,4 +326,8 @@ bool psImageConvolveGetThreads(void); +psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (void); +bool psImageSmooth2dCacheKernel_PS1_V1 (psImageSmooth2dCacheData *smdata, float sigma, float kappa); +bool psImageSmooth2dCache_F32(psImage *image, psImageSmooth2dCacheData *smdata); + /// @} #endif // #ifndef PS_IMAGE_CONVOLVE_H Index: /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c (revision 36277) +++ /branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c (revision 36278) @@ -50,15 +50,38 @@ }; +# 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]); +# 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 (psImage *image, double sigma, double Nsigma) { +psImageSmooth2dCacheData *psImageSmooth2dCacheAlloc (void) { psImageSmooth2dCacheData *smdata = psAlloc(sizeof(psImageSmooth2dCacheData)); @@ -66,34 +89,35 @@ smdata->radflux = NULL; + smdata->Nsigma = 3; + smdata->Ns = -1; return smdata; } -// generate a 2D smoothing kernel for supplied sigma & kappa (PS1_V1 profile). sigma here -// does not need to match that used to allocate the structure, but it is recommended +// 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 size = smdata->Nrange; - - smdata->sigma = sigma; - - smdata->Ns = (int)(smdata->Nsigma * sigma); - smdata->Ns = MAX (3, MIN (smdata->Ns, 10)); + int Ns = (int)(smdata->Nsigma * sigma); + Ns = PS_MAX (3, PS_MIN (Ns, 10)); + smdata->Ns = Ns; int Ns2 = Ns * Ns; - // we are going to use a hard-wired radial profile + // we are going to use a hard-wired set of radial points smdata->radflux = psAlloc(sizeof(float)*NRAD_MAX); float sum = 0.0; - for (i = 0; i < NRAD_MAX; i++) { - float z = radii2[i] / SQ(sigma); - smdata->radflux[i] = 1.0 / (1.0 + kappa*z + pow(z,1.666)); - if (radii2[i] > Ns2) continue; - sum += radiiN[i] * smdata->radflux[i]; - } - for (i = 0; i < NRAD_MAX; i++) { - smdata->radflux[i] = smdata->radflux[i] / sum; + for (int i = 0; i < NRAD_MAX; i++) { + if (radii2[i] > Ns2) { + smdata->radflux[i] = 0.0; + continue; + } + float z = 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; } @@ -105,116 +129,140 @@ { PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_NON_NULL(smdata->radflux, false); + PS_ASSERT_PTR_NON_NULL(smdata->radflux, false); // assert on data type // relevant terms - int Nrange = smdata->Nrange; // Number of pixels either side for convolution kernel - int Nx = smdata->Nx; // Number of columns - int Ny = smdata->Ny; // Number of rows - - psF32 *gauss = &smdata->kernel->data.F32[Nrange]; - psF32 *resultX = smdata->resultX; - psF32 *resultY = smdata->resultY; + int Ns = smdata->Ns; // Number of pixels either side for convolution kernel - /* Smooth in X direction */ - { - for (int j = 0; j < Ny; j++) { - psF32 *vi = image->data.F32[j]; - int xMax = PS_MIN(Nrange, Nx); - /* Smooth first Nrange pixels, with renorm */ - for (int i = 0; i < xMax; i++, vi++) { - int convRange = PS_MIN(Nrange + 1, Nx - i); - psF32 *vr = vi - i; - psF32 *vg = gauss - i; - double g = 0.0; - double s = 0.0; - for (int n = -i; n < convRange; n++, vr++, vg++) { - s += *vg * *vr; - g += *vg; - } - resultX[i] = s / g; - } - /* If that's all the pixels we have, then we're done already */ - if (Nx > Nrange) { - /* Smooth middle pixels; if Nx < 2*Nrange, this pass is skipped */ - for (int i = Nrange; i < Nx - Nrange; i++, vi++) { - psF32 *vr = vi - Nrange; - psF32 *vg = gauss - Nrange; - double s = 0; - for (int n = -Nrange; n < Nrange + 1; n++, vr++, vg++) { - s += *vg * *vr; - } - resultX[i] = s; - } - /* Smooth last Nrange pixels, with renorm */ - // if Nx < 2*Nrange, this pass starts at i == Nrange - int xMin = PS_MAX(Nx - Nrange, Nrange); - for (int i = xMin; i < Nx; i++, vi++) { - psF32 *vr = vi - Nrange; - psF32 *vg = gauss - Nrange; - double g = 0.0; - double s = 0.0; - for (int n = -Nrange; n < Nx - i; n++, vr++, vg++) { - s += *vg * *vr; - g += *vg; - } - resultX[i] = s / g; - } - } - memcpy(image->data.F32[j], resultX, Nx*sizeof(psF32)); - } - } - - // this section probably hits the cache poorly for large images, but is probably OK for small ones - /* Smooth in Y direction */ - { - for (int i = 0; i < Nx; i++) { - int yMax = PS_MIN(Nrange, Ny); - /* Smooth first Nrange pixels, with renorm */ - for (int j = 0; j < yMax; j++) { - int convRange = PS_MIN(Nrange + 1, Ny - j); - psF32 *vg = gauss - j; - double g = 0.0; - double s = 0.0; - for (int n = -j; n < convRange; n++, vg++) { - psF32 vr = image->data.F32[j+n][i]; - s += *vg * vr; - g += *vg; - } - resultY[j] = s / g; - } - /* If that's all the pixels we have, then we're done already */ - if (Ny > Nrange) { - /* Smooth middle pixels */ - for (int j = Nrange; j < Ny - Nrange; j++) { - psF32 *vg = gauss - Nrange; - double s = 0; - for (int n = -Nrange; n < Nrange + 1; n++, vg++) { - psF32 vr = image->data.F32[j+n][i]; - s += *vg * vr; - } - resultY[j] = s; - } - /* Smooth last Nrange pixels, with renorm */ - // if Ny < 2*Nrange, this pass starts at j == Nrange - int yMin = PS_MAX(Ny - Nrange, Nrange); - for (int j = yMin; j < Ny; j++) { - psF32 *vg = gauss - Nrange; - double g = 0.0; - double s = 0.0; - for (int n = -Nrange; n < Ny - j; n++, vg++) { - psF32 vr = image->data.F32[j+n][i]; - s += *vg * vr; - g += *vg; - } - resultY[j] = s / g; - } - } - // loop here - for (int j = 0; j < Ny; j++) { - image->data.F32[j][i] = resultY[j]; - } - } - } + 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; + + 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; } Index: /branches/eam_branches/ipp-20130904/psLib/test/imageops/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/test/imageops/Makefile.am (revision 36277) +++ /branches/eam_branches/ipp-20130904/psLib/test/imageops/Makefile.am (revision 36278) @@ -22,4 +22,5 @@ tap_psImageConvolve \ tap_psImageConvolve2 \ + tap_psImageConvolve2dCache \ tap_psImagePixelExtract \ tap_psImageInterpolate2 \ Index: /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36278) +++ /branches/eam_branches/ipp-20130904/psLib/test/imageops/tap_psImageConvolve2dCache.c (revision 36278) @@ -0,0 +1,99 @@ +/** @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" + +int main(int argc, char **argv) { + + bool status; + + plan_tests(62); + + psMemId id = psMemGetId(); + + // generate an image to convolve + + int Nx = 301; + int Ny = 301; + + int Xo = (Nx - 1) / 2; + int Yo = (Ny - 1) / 2; + + psImage *input = psImageAlloc (Nx, Ny, PS_TYPE_F32); + psImageInit (input, 0.0); + input->data.F32[Xo][Yo] = 1.0; + + psImageSmooth2dCacheData *smdata = psImageSmooth2dCacheAlloc (); + ok(smdata, "init smdata"); + ok(!smdata->radflux, "init radflux to NULL"); + ok(smdata->Nsigma == 3, "init Nsigma to 3"); + + float sigma = 2.0; + float kappa = 1.0; + status = psImageSmooth2dCacheKernel_PS1_V1 (smdata, sigma, kappa); + ok(status, "generate a kernel"); + + status = psImageSmooth2dCache_F32 (input, smdata); + ok(status, "smoothed image"); + + float 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 = PS_SQR((ix + 0.5 - Xo - 0.5) / sigma) + PS_SQR((iy + 0.5 - Yo - 0.5) / sigma); + if (z > 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) < 1e-5, "minor deviations"); + ok (fabs(max_df) < 1e-5, "minor deviations"); + + /* + 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); + ok(!psMemCheckLeaks (id, NULL, stdout, false), "no memory leaks"); + + return exit_status(); +}