IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 18, 2009, 3:19:25 PM (17 years ago)
Author:
Paul Price
Message:

Better method for deciding which image to convolve: generate solutions both ways, and choose the one with the lower residuals. We save the convolutions and least-squares matrices we generate, so that they don't have to be regenerated for the winner. Looks like this works!

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r25101 r25120  
    476476            }
    477477
     478
     479            // Define kernel basis functions
     480            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
     481                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
     482                                                          stamps, footprint, optThreshold, penalty, subMode);
     483                if (!kernels) {
     484                    psErrorClear();
     485                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
     486                }
     487            }
     488            if (kernels == NULL) {
     489                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
     490                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
     491                                                       inner, binning, ringsOrder, penalty, subMode);
     492            }
     493
     494            memCheck("kernels");
     495
    478496            if (subMode == PM_SUBTRACTION_MODE_UNSURE) {
     497#if 0
    479498                // Get backgrounds
    480499                psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background
    … …  
    498517
    499518                pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use
     519#endif
     520                pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej);
    500521                switch (newMode) {
    501522                  case PM_SUBTRACTION_MODE_1:
    … …  
    512533            }
    513534
    514             // Define kernel basis functions
    515             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    516                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    517                                                           stamps, footprint, optThreshold, penalty, subMode);
    518                 if (!kernels) {
    519                     psErrorClear();
    520                     psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
    521                 }
    522             }
    523             if (kernels == NULL) {
    524                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
    525                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    526                                                        inner, binning, ringsOrder, penalty, subMode);
    527             }
    528 
    529             memCheck("kernels");
    530 
    531535            int numRejected = -1;       // Number of rejected stamps in each iteration
    532536            for (int k = 0; k < iter && numRejected != 0; k++) {
    … …  
    565569
    566570                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    567                 numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
     571                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    568572                if (numRejected < 0) {
    569573                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
    … …  
    587591                    goto MATCH_ERROR;
    588592                }
    589                 pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
     593                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    590594                psFree(deviations);
    591595            }
    … …  
    871875
    872876
    873 #if 0
    874 /// A list of stamps
    875 typedef struct {
    876     long num;                           ///< Number of stamps
    877     psArray *stamps;                    ///< The stamps
    878     psArray *regions;                   ///< Regions for each stamp
    879     psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
    880     psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
    881     int footprint;                      ///< Half-size of stamps
    882 } pmSubtractionStampList;
    883 
    884 /// A stamp for image subtraction
    885 typedef struct {
    886     float x, y;                         ///< Position
    887     float flux;                         ///< Flux
    888     float xNorm, yNorm;                 ///< Normalised position
    889     psKernel *image1;                   ///< Reference image postage stamp
    890     psKernel *image2;                   ///< Input image postage stamp
    891     psKernel *variance;                 ///< Variance image postage stamp, or NULL
    892     psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
    893     psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
    894     psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
    895     psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
    896     psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
    897     pmSubtractionStampStatus status;    ///< Status of stamp
    898 } pmSubtractionStamp;
    899 
    900 /// Kernels specification
    901 typedef struct {
    902     pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
    903     psString description;               ///< Description of the kernel parameters
    904     long num;                           ///< Number of kernel components (not including the spatial ones)
    905     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
    906     psVector *widths;                   ///< Gaussian FWHMs (ISIS)
    907     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    908     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
    909     float penalty;                      ///< Penalty for wideness
    910     psVector *penalties;                ///< Penalty for each kernel component
    911     int size;                           ///< The half-size of the kernel
    912     int inner;                          ///< The size of an inner region
    913     int spatialOrder;                   ///< The spatial order of the kernels
    914     int bgOrder;                        ///< The order for the background fitting
    915     pmSubtractionMode mode;             ///< Mode for subtraction
    916     int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    917     psVector *solution1, *solution2;    ///< Solution for the PSF matching
    918     // Quality information
    919     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
    920     int numStamps;                      ///< Number of good stamps
    921 } pmSubtractionKernels;
    922 
    923877// Test a subtraction mode by performing a single iteration
    924878static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode
    925879                                pmSubtractionKernels *kernels, // Kernel description
    926                                 const char *description // Description for trace
     880                                const char *description, // Description for trace
     881                                psImage *subMask,  // Subtraction mask
     882                                float rej               // Rejection threshold
    927883                                )
    928884{
    … …  
    939895    if (!pmSubtractionSolveEquation(kernels, stamps)) {
    940896        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    941         kernels->mode = oldMode;
    942897        return false;
    943898    }
    … …  
    951906
    952907    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
    953     long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
     908    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    954909    if (numRejected < 0) {
    955910        psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
    … …  
    960915
    961916    if (numRejected > 0) {
    962         psTrace("psModules.imcombine", 3, "Solving equation...\n");
     917        // Allow re-fit with reduced stamps set
     918        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    963919        if (!pmSubtractionSolveEquation(kernels, stamps)) {
    964920            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    965921            return false;
    966922        }
     923        psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
    967924        psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    968925        if (!deviations) {
    … …  
    970927            return false;
    971928        }
     929        psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description);
     930        long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
     931        if (numRejected < 0) {
     932            psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
     933            psFree(deviations);
     934            return false;
     935        }
    972936        psFree(deviations);
    973937    }
    … …  
    977941
    978942
    979 pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
    980 {
    981     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, PM_SUBTRACTION_MODE_ERR);
    982     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, PM_SUBTRACTION_MODE_ERR);
    983 
    984     // Copies of the inputs so we can try each way
    985     pmSubtractionStampList *stamps1 = pmSubtractionStampsListCopy(stamps);
    986     pmSubtractionStampList *stamps2 = pmSubtractionStampsListCopy(stamps);
    987 
    988     pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(kernels);
    989     pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(kernels);
    990 
     943pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels,
     944                                        const psImage *subMask, float rej)
     945{
     946    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR);
     947    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR);
     948
     949    // Copies of the inputs
     950    pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps);
     951    pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels);
     952    psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type);
    991953    kernels1->mode = PM_SUBTRACTION_MODE_1;
     954
     955    if (!subtractionModeTest(stamps1, kernels1, "forward", subMask1, rej)) {
     956        psError(PS_ERR_UNKNOWN, false, "Unable to test forward subtraction");
     957        psFree(stamps1);
     958        psFree(kernels1);
     959        psFree(subMask1);
     960        return PM_SUBTRACTION_MODE_ERR;
     961    }
     962    psFree(subMask1);
     963
     964    // Copies of the inputs
     965    pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps);
     966    pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels);
     967    psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type);
    992968    kernels2->mode = PM_SUBTRACTION_MODE_2;
    993969
    994 
    995     subtractionModeTest(stamps1, kernels1, "forward");
    996     subtractionModeTest(stamps2, kernels2, "backward");
    997 
    998     // XXX Compare kernels1->mean, kernels2->mean
    999 }
    1000 
    1001 #endif
     970    if (!subtractionModeTest(stamps2, kernels2, "backward", subMask2, rej)) {
     971        psError(PS_ERR_UNKNOWN, false, "Unable to test backward subtraction");
     972        psFree(stamps2);
     973        psFree(kernels2);
     974        psFree(subMask2);
     975        psFree(stamps1);
     976        psFree(kernels1);
     977        return PM_SUBTRACTION_MODE_ERR;
     978    }
     979    psFree(subMask2);
     980
     981
     982    pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps
     983    pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels
     984    psLogMsg("psModules.imcombine", PS_LOG_INFO,
     985             "Forward: %f +/- %f from %d stamps\nBackward: %f +/- %f from %d stamps\n",
     986             kernels1->mean, kernels1->rms, kernels1->numStamps,
     987             kernels2->mean, kernels2->rms, kernels2->numStamps);
     988
     989    if (kernels1->mean < kernels2->mean) {
     990        bestStamps = stamps1;
     991        bestKernels = kernels1;
     992    } else {
     993        bestStamps = stamps2;
     994        bestKernels = kernels2;
     995    }
     996
     997    psFree(*stamps);
     998    psFree(*kernels);
     999    *stamps = psMemIncrRefCounter(bestStamps);
     1000    *kernels = psMemIncrRefCounter(bestKernels);
     1001
     1002    psFree(stamps1);
     1003    psFree(stamps2);
     1004    psFree(kernels1);
     1005    psFree(kernels2);
     1006
     1007    return bestKernels->mode;
     1008}
     1009
Note: See TracChangeset for help on using the changeset viewer.