Changeset 35767 for trunk/psLib/src/imageops/psImageConvolve.c

Timestamp:

Jul 3, 2013, 2:30:22 PM (13 years ago)

Author:

eugene

Message:

add new psMinimizeLMChi2 function with new convergence criterion; add threaded version of psImageSmoothNoMask

File:

• trunk/psLib/src/imageops/psImageConvolve.c (modified) (4 diffs)

trunk/psLib/src/imageops/psImageConvolve.c

@@ -1231 +1231 @@
 }
 
+/*********************** smooth with mask, threaded version ***********************/
+
 bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image,
                                     const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm,
@@ -1484 +1486 @@
 
 
+/*********************** smooth no-mask ***********************/
+
+bool psImageSmoothNoMask_ScanRows_F32(psImage *calculation, const psImage *image,
+                                      psVector *gaussNorm, float minGauss, 
+                                      int size, int rowStart, int rowStop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numCols = image->numCols;
+    int xLast = numCols - 1; // Last index
+
+    for (int j = rowStart; j < rowStop; j++) {
+        for (int i = 0; i < numCols; i++) {
+            int xMin = PS_MAX(i - size, 0);
+            int xMax = PS_MIN(i + size, xLast);
+            const psF32 *imageData = &image->data.F32[j][xMin];
+            int uMin = - PS_MIN(i, size); /* Minimum kernel index */
+            const psF32 *gaussData = &gauss[uMin];
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+            for (int x = xMin; x <= xMax; x++, imageData++, gaussData++) {
+                sumIG += *imageData * *gaussData;
+                sumG += *gaussData;
+            }
+            if (sumG > minGauss) {
+                /* BW */
+                calculation->data.F32[i][j] = sumIG / sumG;
+            } 
+        }
+    }
+    return true;
+}
+
+bool psImageSmoothNoMask_ScanCols_F32(psImage *output, psImage *calculation, 
+                                      psVector *gaussNorm, float minGauss,
+                                      int size, int colStart, int colStop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numRows = output->numRows;
+    int yLast = numRows - 1; // Last index
+
+    for (int i = colStart; i < colStop; i++) {
+        for (int j = 0; j < numRows; j++) {
+            int yMin = PS_MAX(j - size, 0);
+            int yMax = PS_MIN(j + size, yLast);
+            const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
+            int vMin = - PS_MIN(j, size); /* Minimum kernel index */
+            const psF32 *gaussData = &gauss[vMin];
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+            for (int y = yMin; y <= yMax; y++, imageData++, gaussData++) {
+                sumIG += *imageData * *gaussData;
+                sumG += *gaussData;
+            }
+            output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
+        }
+    }
+    return true;
+}
+
+bool psImageSmoothNoMask_ScanRows_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 7, "programming error: job Nargs mismatch");
+
+    psImage *calculation  = job->args->data[0]; // calculation image
+    const psImage *image  = job->args->data[1]; // input image
+
+    psVector *gaussNorm   = job->args->data[2]; // gauss kernel
+    float minGauss        = PS_SCALAR_VALUE(job->args->data[3],F32);
+    int size              = PS_SCALAR_VALUE(job->args->data[4],S32);
+    int rowStart          = PS_SCALAR_VALUE(job->args->data[5],S32);
+    int rowStop           = PS_SCALAR_VALUE(job->args->data[6],S32);
+    return psImageSmoothNoMask_ScanRows_F32(calculation, image,
+                                          gaussNorm, minGauss, size,
+                                          rowStart, rowStop);
+}
+
+bool psImageSmoothNoMask_ScanCols_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 7, "programming error: job Nargs mismatch");
+
+    psImage *output       = job->args->data[0]; // input image
+    psImage *calculation  = job->args->data[1]; // calculation image
+
+    psVector *gaussNorm   = job->args->data[2]; // gauss kernel
+    float minGauss        = PS_SCALAR_VALUE(job->args->data[3],F32);
+    int size              = PS_SCALAR_VALUE(job->args->data[4],S32);
+    int colStart          = PS_SCALAR_VALUE(job->args->data[5],S32);
+    int colStop           = PS_SCALAR_VALUE(job->args->data[6],S32);
+    return psImageSmoothNoMask_ScanCols_F32(output, calculation,
+                                          gaussNorm, minGauss, size,
+                                          colStart, colStop);
+}
+
+psImage *psImageSmoothNoMask_Threaded(psImage *output, const psImage *image, float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    int nThreads = psThreadPoolSize();
+    int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols;
+    int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows;
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    // Generate normalized gaussian
+    psVector *gaussNorm = NULL;
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    switch (image->type.type) {
+        // Smooth an image with masked pixels
+        // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
+        // [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW"
+
+      case PS_TYPE_F32: {
+          psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
+
+          if (threaded) {
+              /** Smooth in X direction **/
+              for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) {
+                  int rowStop = PS_MIN (rowStart + scanRows, numRows);
+
+                  // allocate a job, construct the arguments for this job
+                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS");
+                  psArrayAdd(job->args, 0, calculation);
+                  psArrayAdd(job->args, 0, (psImage *) image); // cast away const
+                  psArrayAdd(job->args, 0, gaussNorm);
+                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, rowStop,  PS_TYPE_S32);
+                  // -> psImageSmoothNoMask_ScanRows_F32 (calculation, image, gauss, minGauss, size, rowStart, rowStop);
+
+                  // if threading is not active, we simply run the job and return
+                  if (!psThreadJobAddPending(job)) {
+                      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                      psFree(calculation);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+              }
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true, true)) {
+                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                  psFree(calculation);
+                  psFree(gaussNorm);
+                  return false;
+              }
+          } else if (!psImageSmoothNoMask_ScanRows_F32(calculation, image, gaussNorm, minGauss, size, 0, numRows)) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(gaussNorm);
+              return false;
+          }
+
+          output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
+
+          /** Smooth in Y direction  **/
+          if (threaded) {
+              for (int colStart = 0; colStart < numCols; colStart+=scanCols) {
+                  int colStop = PS_MIN (colStart + scanCols, numCols);
+
+                  // allocate a job, construct the arguments for this job
+                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS");
+                  psArrayAdd(job->args, 0, output);
+                  psArrayAdd(job->args, 0, calculation);
+                  psArrayAdd(job->args, 0, gaussNorm);
+                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, colStop,  PS_TYPE_S32);
+                  // -> psImageSmoothMask_ScanCols_F32 (output, calculation, gauss, minGauss, size, colStart, colStop);
+
+                  // if threading is not active, we simply run the job and return
+                  if (!psThreadJobAddPending(job)) {
+                      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                      psFree(calculation);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+              }
+
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true, true)) {
+                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                  psFree(calculation);
+                  psFree(gaussNorm);
+                  return false;
+              }
+          } else if (!psImageSmoothNoMask_ScanCols_F32(output, calculation, gaussNorm, minGauss, size, 0, numCols)) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(gaussNorm);
+              return false;
+          }
+
+          psFree(calculation);
+          psFree(gaussNorm);
+
+          return output;
+      }
+      default:
+        psFree(gaussNorm);
+        char *typeStr;
+        PS_TYPE_NAME(typeStr,image->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                _("Specified psImage type, %s, is not supported."),
+                typeStr);
+        return false;
+    }
+
+    psAbort("Should never reach here.");
+    return false;
+}
+
+/****************************************************************************************/
+
 bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal,
                           double  sigma, double  Nsigma)
@@ -1900 +2122 @@
         }
         {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS", 7);
+            task->function = &psImageSmoothNoMask_ScanRows_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS", 7);
+            task->function = &psImageSmoothNoMask_ScanCols_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
             psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
             task->function = &psImageSmoothMaskPixelsThread;
@@ -1909 +2143 @@
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS");
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS");
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
     }