Changeset 25383 for trunk/psLib

Timestamp:

Sep 15, 2009, 12:45:01 PM (17 years ago)

Author:

Paul Price

Message:

Merging branches/pap (detection efficiency development) into trunk.

Location:

trunk

Files:

• . (modified) (1 prop)
• psLib/src/fits/psFits.c (modified) (4 diffs)
• psLib/src/fits/psFits.h (modified) (1 diff)
• psLib/src/imageops/psImageConvolve.c (modified) (4 diffs)
• psLib/src/imageops/psImageConvolve.h (modified) (1 diff)
• psLib/src/types/psMetadata.c (modified) (4 diffs)
• psLib/src/types/psMetadata.h (modified) (1 diff)

trunk/psLib/src/fits/psFits.c

@@ -344 +344 @@
 // Therefore, we implement our own version of moving to an extension specified by name.  The pure cfitsio
 // version is used if "conventions.compression" handling is turned off in the psFits structure.
-bool psFitsMoveExtName(const psFits* fits,
-                       const char* extname)
+static bool fitsMoveExtName(const psFits* fits, // FITS file
+                            const char* extname, // Extension name
+                            bool errors // Generate errors?
+    )
 {
     PS_ASSERT_FITS_NON_NULL(fits, false);
@@ -356 +358 @@
         // User wants to use cfitsio.  Good luck to them!
         if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
-            psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU '%s'"), extname);
+            }
             return false;
         }
@@ -378 +382 @@
         if (fits_movabs_hdu(fits->fd, i, &hdutype, &status)) {
             // We've run off the end
-            psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            if (errors) {
+                psFitsError(status, true, _("Could not find HDU with %s = '%s'"), extword, extname);
+            }
             return false;
         }
@@ -406 +412 @@
     }
     psAbort("Should never reach here.");
+}
+
+
+bool psFitsMoveExtName(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, true);
+}
+
+bool psFitsMoveExtNameClean(const psFits* fits, const char* extname)
+{
+    return fitsMoveExtName(fits, extname, false);
 }
 

trunk/psLib/src/fits/psFits.h

@@ -245 +245 @@
 );
 
+/** Moves the FITS HDU to the specified extension name without generating errors.
+ *
+ *  @return bool        TRUE if the extension name was found and move was
+ *                      successful, otherwise FALSE
+ */
+bool psFitsMoveExtNameClean(
+    const psFits* fits,                ///< the psFits object to move
+    const char* extname                ///< the extension name
+);
+
 /** Moves the FITS HDU to the specified extension number
  *

trunk/psLib/src/imageops/psImageConvolve.c

@@ -678 +678 @@
 }
 
+static bool imageSmoothMaskPixels(psVector *out, const psImage *image, const psImage *mask,
+                                  psImageMaskType maskVal, const psVector *x, const psVector *y,
+                                  const psVector *gaussNorm, float minGauss, int size, int start, int stop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int xLast = numCols - 1, yLast = numRows - 1; // Last index
+    for (int i = start; i < stop; i++) {
+        int xPix = x->data.S32[i], yPix = y->data.S32[i]; // Pixel coordinates for smoothing
+
+        int yMin = PS_MAX(yPix - size, 0);
+        int yMax = PS_MIN(yPix + size, yLast);
+        int xMin = PS_MAX(xPix - size, 0);
+        int xMax = PS_MIN(xPix + size, xLast);
+
+        const float *yGauss = &gauss[yMin - yPix];
+
+        double ySumIG = 0.0, ySumG = 0.0;
+        for (int v = yMin; v <= yMax; v++, yGauss++) {
+            const float *xGauss = &gauss[xMin - xPix];
+            double xSumIG = 0.0, xSumG = 0.0;
+            const psImageMaskType *maskData = &mask->data.PS_TYPE_IMAGE_MASK_DATA[v][xMin];
+            const psF32 *imageData = &image->data.F32[v][xMin];
+            for (int u = xMin; u <= xMax; u++, xGauss++, imageData++, maskData++) {
+                if (*maskData & maskVal) {
+                    continue;
+                }
+                xSumIG += *imageData * *xGauss;
+                xSumG += *xGauss;
+            }
+            if (xSumG > minGauss) {
+                ySumIG += xSumIG * *yGauss;
+                ySumG += xSumG * *yGauss;
+            }
+        }
+
+        out->data.F32[i] = ySumG > minGauss ? ySumIG / ySumG : NAN;
+    }
+
+    return true;
+}
+
+static bool psImageSmoothMaskPixelsThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "programming error: no job arguments");
+    psAssert(job->args->n == 11, "programming error: wrong number of job arguments");
+
+    psVector *out = job->args->data[0]; // Output vector
+    const psImage *image  = job->args->data[1]; // Input image
+    const psImage *mask   = job->args->data[2]; // Input mask
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3], PS_TYPE_IMAGE_MASK_DATA);
+    const psVector *x = job->args->data[4];
+    const psVector *y = job->args->data[5];
+    const psVector *gaussNorm = job->args->data[6];
+    float minGauss = PS_SCALAR_VALUE(job->args->data[7], F32);
+    int size = PS_SCALAR_VALUE(job->args->data[8], S32);
+    int start = PS_SCALAR_VALUE(job->args->data[9], S32);
+    int stop = PS_SCALAR_VALUE(job->args->data[10], S32);
+    return imageSmoothMaskPixels(out, image, mask, maskVal, x, y, gaussNorm,
+                                 minGauss, size, start, stop);
+}
+
+
+psVector *psImageSmoothMaskPixels(const psImage *image, const psImage *mask, psImageMaskType maskVal,
+                                  const psVector *x, const psVector *y,
+                                  float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(image, mask, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    int num = x->n;                     // Number of pixels to smooth
+    psVector *out = psVectorAlloc(num, PS_TYPE_F32); // Output results
+
+    // Generate normalized gaussian
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    // Columns
+    if (threaded) {
+        int numThreads = psThreadPoolSize(); // Number of threads
+        int delta = (numThreads) ? num / numThreads + 1 : num; // Block of cols to do at once
+        for (int start = 0; start < num; start += delta) {
+            int stop = PS_MIN(start + delta, num);  // End of block
+
+            psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
+            psArrayAdd(job->args, 1, out);
+            psArrayAdd(job->args, 1, (psImage*)image);
+            psArrayAdd(job->args, 1, (psImage*)mask);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
+            psArrayAdd(job->args, 1, (psVector*)x);
+            psArrayAdd(job->args, 1, (psVector*)y);
+            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, start, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, stop, PS_TYPE_S32);
+            if (!psThreadJobAddPending(job)) {
+                psFree(job);
+                psFree(gaussNorm);
+                psFree(out);
+                return NULL;
+            }
+            psFree(job);
+        }
+    } 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);
+        return NULL;
+    }
+
+    psFree(gaussNorm);
+    return out;
+}
+
+
 // Smooth an image with masked pixels
 // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
@@ -1325 +1458 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaRows = (numThreads) ? numRows / numThreads : numRows; // Block of rows to do at once
+        int deltaRows = (numThreads) ? numRows / numThreads + 1 : numRows; // Block of rows to do at once
         for (int start = 0; start < numRows; start += deltaRows) {
-            int stop = PS_MIN (start + deltaRows, numRows);  // end of row block
+            int stop = PS_MIN(start + deltaRows, numRows);  // end of row block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1363 +1496 @@
     if (threaded) {
         int numThreads = psThreadPoolSize(); // Number of threads
-        int deltaCols = (numThreads) ? numCols / numThreads : numCols; // Block of cols to do at once
+        int deltaCols = (numThreads) ? numCols / numThreads + 1 : numCols; // Block of cols to do at once
         for (int start = 0; start < numCols; start += deltaCols) {
-            int stop = PS_MIN (start + deltaCols, numCols);  // end of col block
+            int stop = PS_MIN(start + deltaCols, numCols);  // end of col block
 
             psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_CONVOLVE_MASK");
@@ -1486 +1619 @@
             psFree(task);
         }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
+            task->function = &psImageSmoothMaskPixelsThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
     } else if (!set && threaded) {
         psThreadTaskRemove("PSLIB_IMAGE_CONVOLVE_MASK");

trunk/psLib/src/imageops/psImageConvolve.h

@@ -222 +222 @@
     );
 
+/// Smooth particular pixels on an image, allowing for masked pixels
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
+psVector *psImageSmoothMaskPixels(
+    const psImage *image,               ///< Input image (F32)
+    const psImage *mask,                ///< Mask image
+    psImageMaskType maskVal,            ///< Value to mask
+    const psVector *x,                  ///< x coordinates
+    const psVector *y,                  ///< y coordinates
+    float sigma,                        ///< Width of the smoothing kernel (pixels)
+    float numSigma,                     ///< Size of the smoothing box (sigma)
+    float minGauss                      ///< Minimum fraction of Gaussian to accept
+    );
 
 psImage *psImageSmoothMask_Threaded(psImage *output,

trunk/psLib/src/types/psMetadata.c

@@ -996 +996 @@
     if (metadataItem->type == PS_DATA_METADATA_MULTI) {
         // if multiple keys found, use the first.
-        //        metadataItem = (psMetadataItem*)((metadataItem->data.list)->head);
         metadataItem = (psMetadataItem*)(metadataItem->data.list->head->data);
         if (status) {
@@ -1082 +1081 @@
     } \
     \
-    /* psFree(metadataItem); currently, the lookup doesn't increment the ref count */ \
     return value; \
 }
@@ -1101 +1099 @@
 psMetadataLookupNumTYPE(VectorMaskType,VectorMask)
 psMetadataLookupNumTYPE(ImageMaskType,ImageMask)
+
+#define psMetadataLookupPtrTYPE(TYPENAME,NAME,TYPE,VAL) \
+TYPENAME psMetadataLookup##NAME(bool *status, const psMetadata *md, const char *key) \
+{ \
+    PS_ASSERT_METADATA_NON_NULL(md, NULL); \
+    PS_ASSERT_STRING_NON_EMPTY(key, NULL); \
+    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); /* The item of interest */ \
+    if (!item) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (item->type == PS_DATA_METADATA_MULTI) { \
+        /* if multiple keys found, use the first. */ \
+        item = item->data.list->head->data; \
+    } \
+    if (item->type != TYPE) { \
+        if (status) { \
+            *status = false; \
+        } else { \
+            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key); \
+        } \
+        return NULL; \
+    } \
+    \
+    if (status) { \
+        *status = true; \
+    } \
+    return item->data.VAL; \
+}
+
+psMetadataLookupPtrTYPE(psMetadata*, Metadata, PS_DATA_METADATA, md)
+psMetadataLookupPtrTYPE(psString, Str, PS_DATA_STRING, str)
+psMetadataLookupPtrTYPE(psTime*, Time, PS_DATA_TIME, V)
+psMetadataLookupPtrTYPE(psVector*, Vector, PS_DATA_VECTOR, V)
+
 
 psMetadataItem* psMetadataGet(const psMetadata *md,
@@ -1257 +1295 @@
 }
 
-psMetadata *psMetadataLookupMetadata(bool *status,
-                                     const psMetadata *md,
-                                     const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    psMetadata *value = NULL;  // The value to return
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_METADATA) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.md; // The requested metadata
-    }
-    return value;
-}
-
-psTime *psMetadataLookupTime(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key);
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-        return NULL;
-    }
-
-    if (item->type != PS_DATA_TIME) {
-        // The value at the key isn't metadata
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_TIME, as expected.\n", key);
-        }
-        return NULL;
-    }
-
-    // We have the requested metadata
-    if (status) {
-        *status = true;
-    }
-
-    return item->data.V;
-}
-
-
-psString psMetadataLookupStr(bool *status,
-                             const psMetadata *md,
-                             const char *key)
-{
-    PS_ASSERT_METADATA_NON_NULL(md,NULL);
-
-    psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments
-    //    char *value = NULL;   // The value to return
-    psString value = NULL;
-    if (!item) {
-        // The given key isn't in the metadata
-        if (status) {
-            *status = false;
-        } else {
-            psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key);
-        }
-    } else if (item->type != PS_DATA_STRING) {
-        // The value at the key isn't of the desired type
-        if (status) {
-            *status = false;
-        } else {
-            psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_STRING, as expected.\n", key);
-        }
-        //        value = NULL;
-    } else {
-        // We have the requested metadata
-        if (status) {
-            *status = true;
-        }
-        value = item->data.V;
-    }
-    return value;
-}
+
+
+
 
 psList *psMetadataKeys(psMetadata *md)

trunk/psLib/src/types/psMetadata.h

@@ -717 +717 @@
 PS_METADATA_LOOKUP_TYPE_DECL(Ptr, psPtr);
 PS_METADATA_LOOKUP_TYPE_DECL(Str, psString);
+PS_METADATA_LOOKUP_TYPE_DECL(Vector, psVector*);
 PS_METADATA_LOOKUP_TYPE_DECL(Metadata, psMetadata*);
 PS_METADATA_LOOKUP_TYPE_DECL(Time, psTime*);