Changeset 36278

Timestamp:

Nov 11, 2013, 10:13:29 AM (13 years ago)

Author:

eugene

Message:

finish psImageConvolve2dCache, and add test

Location:

branches/eam_branches/ipp-20130904/psLib

Files:

• src/imageops/Makefile.am (modified) (1 diff)
• src/imageops/psImageConvolve.h (modified) (1 diff)
• src/imageops/psImageConvolve2dCache.c (modified) (3 diffs)
• test/imageops/Makefile.am (modified) (1 diff)
• test/imageops/tap_psImageConvolve2dCache.c (added)

branches/eam_branches/ipp-20130904/psLib/src/imageops/Makefile.am

@@ -7 +7 @@
         psImageBackground.c \
         psImageConvolve.c \
+        psImageConvolve2dCache.c \
         psImageCovariance.c \
         psImageGeomManip.c \

branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve.h

@@ -326 +326 @@
 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

branches/eam_branches/ipp-20130904/psLib/src/imageops/psImageConvolve2dCache.c

@@ -50 +50 @@
 };
 
+# 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 +89 @@
 
     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 +129 @@
 {
     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;
 }

branches/eam_branches/ipp-20130904/psLib/test/imageops/Makefile.am

@@ -22 +22 @@
         tap_psImageConvolve \
         tap_psImageConvolve2 \
+        tap_psImageConvolve2dCache \
         tap_psImagePixelExtract \
         tap_psImageInterpolate2 \