Changeset 36277

Timestamp:

Nov 10, 2013, 6:05:54 AM (13 years ago)

Author:

eugene

Message:

working on 2d convolve

Location:

branches/eam_branches/ipp-20130904

Files:

• psLib/src/imageops/psImageConvolve2dCache.c (modified) (1 diff)
• psModules/src/objects/pmPCMdata.c (modified) (2 diffs)

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

@@ -65 +65 @@
     psMemSetDeallocator(smdata, (psFreeFunc) psImageSmooth2dCacheDataFree);
 
-    smdata->kernel = NULL;
-
-    if (!image) {
-        // relevant terms
-        smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
-        smdata->Nx = 0;
-        smdata->Ny = 0;
-        smdata->resultX = NULL;
-        smdata->resultY = NULL;
-        return smdata;
-    }
+    smdata->radflux = NULL;
+
+    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
+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 Ns2 = Ns * Ns;
+
+    // we are going to use a hard-wired radial profile
+    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;
+    }
+
+    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_NON_NULL(smdata->radflux, false);
+    // assert on data type
 
     // relevant terms
-    smdata->Nrange = sigma*Nsigma + 0.5;    // Number of pixels either side for convolution kernel
-    smdata->Nx = image->numCols;            // Number of columns
-    smdata->Ny = image->numRows;            // Number of rows
-
-    // use a temp running buffer for X and Y directions.
-    smdata->resultX = psAlloc(smdata->Nx * sizeof(psF32));
-    memset (smdata->resultX, 0, smdata->Nx*sizeof(psF32));
-
-    smdata->resultY = psAlloc(smdata->Ny * sizeof(psF32));
-    memset (smdata->resultY, 0, smdata->Ny*sizeof(psF32));
-
-    return smdata;
-}
+    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;
+       
+    /* 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];
+            }
+        }
+    }
+    return true;
+}
+

branches/eam_branches/ipp-20130904/psModules/src/objects/pmPCMdata.c

@@ -59 +59 @@
     psFree (pcm->constraint);
 
-    psFree (pcm->smdata1); // pre-allocated data for psImageSmooth_PreAlloc
-    psFree (pcm->smdata2); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata); // pre-allocated data for psImageSmooth_PreAlloc
+    psFree (pcm->smdata2d); // pre-allocated data for psImageSmooth_PreAlloc
     return;
 }
@@ -90 +90 @@
     }
 
-    pcm->smdata1 = NULL;
-    pcm->smdata2 = NULL;
+    pcm->smdata = NULL;
+    pcm->smdata2d = NULL;
 
     pcm->modelConv = NULL;