Changeset 25101

Timestamp:

Aug 17, 2009, 6:55:35 PM (17 years ago)

Author:

Paul Price

Message:

Working on more robust way of determining the image to convolve. Intend to attempt to measure the convolution kernel each way, and take the one that does the best job. This code is currently #if-ed out.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtractionKernels.c (modified) (1 diff)
• pmSubtractionKernels.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (1 diff)
• pmSubtractionStamps.c (modified) (1 diff)
• pmSubtractionStamps.h (modified) (1 diff)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -765 +765 @@
     return PM_SUBTRACTION_KERNEL_NONE;
 }
+
+pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(in, NULL);
+
+    pmSubtractionKernels *out = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionKernelsFree);
+
+    out->type = in->type;
+    out->description = in->description;
+    out->num = in->num;
+    out->u = psMemIncrRefCounter(in->u);
+    out->v = psMemIncrRefCounter(in->v);
+    out->widths = psMemIncrRefCounter(in->widths);
+    out->preCalc = psMemIncrRefCounter(in->preCalc);
+    out->penalty = in->penalty;
+    out->penalties = psMemIncrRefCounter(in->penalties);
+    out->uStop = psMemIncrRefCounter(in->uStop);
+    out->vStop = psMemIncrRefCounter(in->vStop);
+    out->size = in->size;
+    out->inner = in->inner;
+    out->spatialOrder = in->spatialOrder;
+    out->bgOrder = in->bgOrder;
+    out->mode = in->mode;
+    out->numCols = in->numCols;
+    out->numRows = in->numRows;
+    out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
+    out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
+
+    return out;
+}

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -203 +203 @@
     );
 
+/// Copy kernels
+///
+/// A deep copy is performed on the solution only; the other components are merely pointers.
+pmSubtractionKernels *pmSubtractionKernelsCopy(
+    const pmSubtractionKernels *in      // Kernels to copy
+    );
+
 
 #endif

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -869 +869 @@
     return mode;
 }
+
+
+#if 0
+/// A list of stamps
+typedef struct {
+    long num;                           ///< Number of stamps
+    psArray *stamps;                    ///< The stamps
+    psArray *regions;                   ///< Regions for each stamp
+    psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
+    psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
+    int footprint;                      ///< Half-size of stamps
+} pmSubtractionStampList;
+
+/// A stamp for image subtraction
+typedef struct {
+    float x, y;                         ///< Position
+    float flux;                         ///< Flux
+    float xNorm, yNorm;                 ///< Normalised position
+    psKernel *image1;                   ///< Reference image postage stamp
+    psKernel *image2;                   ///< Input image postage stamp
+    psKernel *variance;                 ///< Variance image postage stamp, or NULL
+    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
+    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
+    psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
+    psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
+    psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
+    pmSubtractionStampStatus status;    ///< Status of stamp
+} pmSubtractionStamp;
+
+/// Kernels specification
+typedef struct {
+    pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
+    psString description;               ///< Description of the kernel parameters
+    long num;                           ///< Number of kernel components (not including the spatial ones)
+    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
+    psVector *widths;                   ///< Gaussian FWHMs (ISIS)
+    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
+    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
+    float penalty;                      ///< Penalty for wideness
+    psVector *penalties;                ///< Penalty for each kernel component
+    int size;                           ///< The half-size of the kernel
+    int inner;                          ///< The size of an inner region
+    int spatialOrder;                   ///< The spatial order of the kernels
+    int bgOrder;                        ///< The order for the background fitting
+    pmSubtractionMode mode;             ///< Mode for subtraction
+    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
+    psVector *solution1, *solution2;    ///< Solution for the PSF matching
+    // Quality information
+    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
+    int numStamps;                      ///< Number of good stamps
+} pmSubtractionKernels;
+
+// Test a subtraction mode by performing a single iteration
+static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
+                                pmSubtractionKernels *kernels, // Kernel description
+                                const char *description // Description for trace
+                                )
+{
+    assert(stamps);
+    assert(kernels);
+
+    psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        kernels->mode = oldMode;
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
+    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
+    if (numRejected < 0) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
+        psFree(deviations);
+        return false;
+    }
+    psFree(deviations);
+
+    if (numRejected > 0) {
+        psTrace("psModules.imcombine", 3, "Solving equation...\n");
+        if (!pmSubtractionSolveEquation(kernels, stamps)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+            return false;
+        }
+        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
+        if (!deviations) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+            return false;
+        }
+        psFree(deviations);
+    }
+
+    return true;
+}
+
+
+pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, PM_SUBTRACTION_MODE_ERR);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, PM_SUBTRACTION_MODE_ERR);
+
+    // Copies of the inputs so we can try each way
+    pmSubtractionStampList *stamps1 = pmSubtractionStampsListCopy(stamps);
+    pmSubtractionStampList *stamps2 = pmSubtractionStampsListCopy(stamps);
+
+    pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(kernels);
+    pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(kernels);
+
+    kernels1->mode = PM_SUBTRACTION_MODE_1;
+    kernels2->mode = PM_SUBTRACTION_MODE_2;
+
+
+    subtractionModeTest(stamps1, kernels1, "forward");
+    subtractionModeTest(stamps2, kernels2, "backward");
+
+    // XXX Compare kernels1->mean, kernels2->mean
+}
+
+#endif

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -225 +225 @@
 }
 
+pmSubtractionStampList *pmSubtractionStampListCopy(const pmSubtractionStampList *in)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(in, NULL);
+
+    pmSubtractionStampList *out = psAlloc(sizeof(pmSubtractionStampList)); // Copied stamp list to return
+    psMemSetDeallocator(out, (psFreeFunc)subtractionStampListFree);
+
+    int num = out->num = in->num;       // Number of stamps
+    out->stamps = psArrayAlloc(num);
+    out->regions = psArrayAlloc(num);
+    out->x = NULL;
+    out->y = NULL;
+    out->flux = NULL;
+    out->footprint = in->footprint;
+
+    for (int i = 0; i < num; i++) {
+        psRegion *inRegion = in->regions->data[i]; // Input region
+        out->regions->data[i] = psRegionAlloc(inRegion->x0, inRegion->x1, inRegion->y0, inRegion->y1);
+
+        pmSubtractionStamp *inStamp = in->stamps->data[i]; // Input stamp
+        pmSubtractionStamp *outStamp = psAlloc(sizeof(pmSubtractionStamp));
+        psMemSetDeallocator(outStamp, (psFreeFunc)subtractionStampFree);
+        outStamp->x = inStamp->x;
+        outStamp->y = inStamp->y;
+        outStamp->flux = inStamp->flux;
+        outStamp->xNorm = inStamp->xNorm;
+        outStamp->yNorm = inStamp->yNorm;
+        outStamp->status = inStamp->status;
+
+        outStamp->image1 = psKernelCopy(inStamp->image1);
+        outStamp->image2 = psKernelCopy(inStamp->image2);
+        outStamp->variance = psKernelCopy(inStamp->variance);
+
+        if (inStamp->convolutions1) {
+            int size = inStamp->convolutions1->n; // Size of array
+            outStamp->convolutions1 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                outStamp->convolutions1->data[j] = psKernelCopy(inStamp->convolutions1->data[j]);
+            }
+        }
+        if (inStamp->convolutions2) {
+            int size = inStamp->convolutions2->n; // Size of array
+            outStamp->convolutions2 = psArrayAlloc(size);
+            for (int j = 0; j < size; j++) {
+                outStamp->convolutions2->data[j] = psKernelCopy(inStamp->convolutions2->data[j]);
+            }
+        }
+
+        outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
+        outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
+        outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
+        outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
+        outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
+
+        out->stamps->data[i] = outStamp;
+    }
+
+    return out;
+}
+
 pmSubtractionStamp *pmSubtractionStampAlloc(void)
 {

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -52 +52 @@
     } \
 }
+
+/// Copy a list of stamps
+///
+/// A deep copy is performed of the stamp list and the component stamps
+pmSubtractionStampList *pmSubtractionStampListCopy(
+    const pmSubtractionStampList *in    // Stamp list to copy
+    );
+
 
 /// A stamp for image subtraction