Changeset 28017 for branches/czw_branch/20100427/psLib

Timestamp:

May 18, 2010, 5:06:01 PM (16 years ago)

Author:

watersc1

Message:

Almost working version of the mask-stats/software revision code.

Almost.

Location:

branches/czw_branch/20100427

Files:

• . (modified) (1 prop)
• psLib/src/imageops/psImageConvolve.c (modified) (8 diffs)
• psLib/src/imageops/psImageCovariance.c (modified) (3 diffs)
• psLib/src/imageops/psImageCovariance.h (modified) (1 diff)
• psLib/src/mathtypes/psImage.h (modified) (1 diff)
• psLib/src/sys/psConfigure.c (modified) (1 diff)
• psLib/src/sys/psConfigure.h (modified) (1 diff)
• psLib/src/types/psMetadataHeader.c (modified) (2 diffs)

branches/czw_branch/20100427/psLib/src/imageops/psImageConvolve.c

@@ -37 +37 @@
 
 
+
 static bool threaded = false;           // Run image convolution threaded?
-
-
+static pthread_mutex_t threadMutex = PTHREAD_MUTEX_INITIALIZER;
 
 
@@ -871 +871 @@
             psFree(job);
         }
+        if (!psThreadPoolWait(true)) {
+            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+            psFree(gaussNorm);
+            psFree(out);
+            return NULL;
+        }
     } else if (!imageSmoothMaskPixels(out, image, mask, maskVal, x, y,
                                       gaussNorm, minGauss, size, 0, num)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to smooth pixels.");
-        psFree(gaussNorm);
-        psFree(out);
-        return NULL;
-    }
-
-    if (threaded && !psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         psFree(gaussNorm);
         psFree(out);
@@ -1192 +1191 @@
           psImage *calcMask = psImageAlloc(numRows, numCols, PS_TYPE_IMAGE_MASK); /* Mask for calculation image; BW */
 
-          /** 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_SMOOTHMASK_SCANROWS");
-              psArrayAdd(job->args, 1, calculation);
-              psArrayAdd(job->args, 1, calcMask);
-              psArrayAdd(job->args, 1, (psImage *) image); // cast away const
-              psArrayAdd(job->args, 1, (psImage *) mask); // cast away const
-              PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
-              psArrayAdd(job->args, 1, 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);
-              // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, gauss, minGauss, size, rowStart, rowStop);
-
-              // if threading is not active, we simply run the job and return
-              if (!psThreadJobAddPending(job)) {
+          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_SMOOTHMASK_SCANROWS");
+                  psArrayAdd(job->args, 1, calculation);
+                  psArrayAdd(job->args, 1, calcMask);
+                  psArrayAdd(job->args, 1, (psImage *) image); // cast away const
+                  psArrayAdd(job->args, 1, (psImage *) mask); // cast away const
+                  PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+                  psArrayAdd(job->args, 1, 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);
+                  // -> psImageSmoothMask_ScanRows_F32 (calculation, calcMask, image, mask, maskVal, 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(job);
+                      psFree(calculation);
+                      psFree(calcMask);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+                  psFree(job);
+              }
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true)) {
                   psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
-                  psFree(job);
+                  psFree(calculation);
+                  psFree(calcMask);
+                  psFree(gaussNorm);
                   return false;
               }
-              psFree(job);
-
-          }
-
-          // wait here for the threaded jobs to finish (NOP if threading is not active)
-          if (!psThreadPoolWait(true)) {
+          } else if (!psImageSmoothMask_ScanRows_F32(calculation, calcMask, image, mask, maskVal,
+                                                     gaussNorm, minGauss, size, 0, numRows)) {
               psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(calcMask);
+              psFree(gaussNorm);
               return false;
           }
@@ -1229 +1241 @@
 
           /** Smooth in Y direction  **/
-          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_SMOOTHMASK_SCANCOLS");
-              psArrayAdd(job->args, 1, output);
-              psArrayAdd(job->args, 1, calculation);
-              psArrayAdd(job->args, 1, calcMask);
-              PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
-              psArrayAdd(job->args, 1, 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, calcMask, maskVal, gauss, minGauss, size, colStart, colStop);
-
-              // if threading is not active, we simply run the job and return
-              if (!psThreadJobAddPending(job)) {
+          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_SMOOTHMASK_SCANCOLS");
+                  psArrayAdd(job->args, 1, output);
+                  psArrayAdd(job->args, 1, calculation);
+                  psArrayAdd(job->args, 1, calcMask);
+                  PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+                  psArrayAdd(job->args, 1, 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, calcMask, maskVal, 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(job);
+                      psFree(calculation);
+                      psFree(calcMask);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+                  psFree(job);
+              }
+
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true)) {
                   psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
-                  psFree(job);
+                  psFree(calculation);
+                  psFree(calcMask);
+                  psFree(gaussNorm);
                   return false;
               }
-              psFree(job);
-          }
-
-          // wait here for the threaded jobs to finish (NOP if threading is not active)
-          if (!psThreadPoolWait(true)) {
+          } else if (!psImageSmoothMask_ScanCols_F32(output, calculation, calcMask, maskVal,
+                                                     gaussNorm, minGauss, size, 0, numCols)) {
               psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(calcMask);
+              psFree(gaussNorm);
               return false;
           }
+
           psFree(calculation);
           psFree(calcMask);
@@ -1559 +1587 @@
             psFree(job);
         }
+        if (!psThreadPoolWait(true)) {
+            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+            psFree(conv);
+            psFree(out);
+            return NULL;
+        }
     } else if (!imageConvolveMaskColumns(conv, mask, 0, numRows, maskVal, xMin, xMax)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask columns.");
-        psFree(conv);
-        psFree(out);
-        return NULL;
-    }
-
-    if (threaded && !psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         psFree(conv);
         psFree(out);
@@ -1597 +1624 @@
             psFree(job);
         }
+        if (!psThreadPoolWait(true)) {
+            psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+            psFree(conv);
+            psFree(out);
+            return NULL;
+        }
     } else if (!imageConvolveMaskRows(out, conv, 0, numCols, setVal, yMin, yMax)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask columns.");
-        psFree(conv);
-        psFree(out);
-        return NULL;
-    }
-
-    if (threaded && !psThreadPoolWait(true)) {
-        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
         psFree(conv);
         psFree(out);
@@ -1679 +1705 @@
 bool psImageConvolveSetThreads(bool set)
 {
+    pthread_mutex_lock(&threadMutex);
     bool old = threaded;                // Old value
     if (set && !threaded) {
@@ -1711 +1738 @@
     }
     threaded = set;
+    pthread_mutex_unlock(&threadMutex);
     return old;
 }

branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.c

@@ -11 +11 @@
 #include "psMemory.h"
 #include "psConstants.h"
+#include "psImageStructManip.h"
+#include "psImagePixelManip.h"
 #include "psImageConvolve.h"
 #include "psTrace.h"
@@ -16 +18 @@
 #include "psScalar.h"
 #include "psThread.h"
+#include "psImageInterpolate.h"
 
 #include "psImageCovariance.h"
 
 static bool threaded = false;           // Run threaded?
-
 
 psKernel *psImageCovarianceNone(void)
@@ -530 +532 @@
 
 
+psKernel *psImageCovarianceScale(const psKernel *in, float scale)
+{
+    // Trivial cases
+    if (!in) {
+        psKernel *out = psKernelAlloc(0, 0, 0, 0); // Output covariance
+        out->kernel[0][0] = 1.0;
+        return out;
+    }
+    PS_ASSERT_KERNEL_NON_NULL(in, NULL);
+    if (scale == 1.0) {
+        psImage *copy = psImageCopy(NULL, in->image, PS_TYPE_F32); // Copy of input covariance
+        psKernel *out = psKernelAllocFromImage(copy, -in->xMin, -in->yMin); // Output covariance
+        psFree(copy);
+        return out;
+    }
+
+    int xMinIn = in->xMin, xMaxIn = in->xMax, yMinIn = in->yMin, yMaxIn = in->yMax; // Input size
+    int xMinOut = (float)xMinIn / scale - 0.5, xMaxOut = (float)xMaxIn / scale + 0.5;     // Output size in x
+    int yMinOut = (float)yMinIn / scale - 0.5, yMaxOut = (float)yMaxIn / scale + 0.5;     // Output size in y
+
+    // Over-fill the covariance matrix so we're not troubled by edge effects
+    psKernel *overfill = psKernelAlloc(xMinIn - 1, xMaxIn + 1, yMinIn - 1, yMaxIn + 1); // Overfilled covar
+    psImageInit(overfill->image, 0.0);
+    int numOverlay = (xMaxIn - xMinIn + 1) * (yMaxIn - yMinIn + 1); // Number of pixels to overlay
+    if (psImageOverlaySection(overfill->image, in->image, 1, 1, "=") != numOverlay) {
+        psError(psErrorCodeLast(), false, "Unable to overfill covariance matrix.");
+        psFree(overfill);
+        return NULL;
+    }
+
+    psImageInterpolation *interp = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, overfill->image,
+                                                             NULL, NULL, 0, NAN, NAN, 0xFF, 0xFF,
+                                                             0.0, 0); // Interpolation
+    psFree(overfill);
+
+    // In transforming the positions, we get +0.5 to account for the centre of the pixels being at 0.5
+    // and +1 to account for the overfill.
+
+    psKernel *out = psKernelAlloc(xMinOut, xMaxOut, yMinOut, yMaxOut); // Output covariance
+    for (int y = yMinOut; y <= yMaxOut; y++) {
+        float yIn = y * scale + 0.5 - yMinIn + 1; // Position on input image (not the kernel)
+        for (int x = xMinOut; x <= xMaxOut; x++) {
+            float xIn = x * scale + 0.5 - xMinIn + 1; // Position on input (not the kernel)
+            double value;                                     // Value on output
+            if (!psImageInterpolate(&value, NULL, NULL, xIn, yIn, interp)) {
+                psError(psErrorCodeLast(), false, "Unable to interpolate kernel.");
+                return false;
+            }
+            out->kernel[y][x] = value;
+        }
+    }
+
+    psFree(interp);
+
+    return out;
+}
+
+
 bool psImageCovarianceSetThreads(bool set)
 {

branches/czw_branch/20100427/psLib/src/imageops/psImageCovariance.h

@@ -90 +90 @@
     );
 
+
+/// Rescale a covariance matrix following a change in plate scale
+///
+/// The covariance matrix is stretched or shrunk to match the new plate scale.
+psKernel *psImageCovarianceScale(
+    const psKernel *in,                 ///< Input covariance pseudo-matrix
+    float scale                         ///< Scale factor (output plate scale relative to input plate scale)
+    );
+
 /// Control threading for image covariance functions
 ///

branches/czw_branch/20100427/psLib/src/mathtypes/psImage.h

@@ -66 +66 @@
 #define P_PSIMAGE_SET_ROW0(img,r0) {*(int*)&img->row0 = r0;}
 #define P_PSIMAGE_SET_TYPE(img,t) {*(psMathType*)&img->type = t;}
+#define P_PSIMAGE_GET_TYPE(img) ((img)->type->type)
 
 /** Create an image of the specified size and type.

branches/czw_branch/20100427/psLib/src/sys/psConfigure.c

@@ -66 +66 @@
 }
 
+psString psLibRevision(void)
+{
+    char *value = NULL;
+    psStringAppend(&value, "%s", PSLIB_VERSION);
+    return value;
+}
+
 psString psLibSource(void)
 {

branches/czw_branch/20100427/psLib/src/sys/psConfigure.h

@@ -32 +32 @@
  */
 psString psLibVersion(void);
+
+/** Get current psLib revision number
+ *
+ *  Returns the current psLib revision number as a string.
+ *
+ *  @return psString: String with revision number.
+ */
+psString psLibRevision(void);
 
 /** Get current psLib source

branches/czw_branch/20100427/psLib/src/types/psMetadataHeader.c

@@ -17 +17 @@
     psString version = psLibVersion();  // Software version
     psString source = psLibSource();    // Software source
-
-    psMetadataAddStr(header, PS_LIST_TAIL, "PSLIB_V", 0, NULL, source);
+    psString revision = psLibRevision();
+    psMetadataAddStr(header, PS_LIST_TAIL, "PSLIB_V", PS_META_REPLACE, NULL, revision);
     
     psStringPrepend(&version, "psLib version: ");
@@ -28 +28 @@
     psFree(version);
     psFree(source);
-
+    psFree(revision);
     return true;
 }