IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19532 for trunk/psModules


Ignore:
Timestamp:
Sep 11, 2008, 6:12:42 PM (18 years ago)
Author:
Paul Price
Message:

Correcting stack variances for the imperfect convolution. Changed (again) the way the deviations are calculated; they're now straight chi2, which allows us to estimate how the variance factor needs to be changed in order to correct for it. Moving additional values into pmSubtractionKernels so that this can be calculated from a kernel that has been read in (very useful for debugging, also nice to have a record of the quality of the fit).

Location:
trunk/psModules/src/imcombine
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r19487 r19532  
    88 *  @author GLG, MHPCC
    99 *
    10  *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
    11  *  @date $Date: 2008-09-11 20:52:16 $
     10 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
     11 *  @date $Date: 2008-09-12 04:12:42 $
    1212 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    1313 *
    … …  
    3030#define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
    3131#define PIXEL_MAP_BUFFER 2              // Number of entries to add to pixel map at a time
    32 //#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
     32#define VARIANCE_FACTORS                // Use variance factors when calculating the variances?
    3333#define NUM_DIRECT_STDEV 5              // For less than this number of values, measure stdev directly
    3434
    … …  
    312312              // Use variance to check that the two are consistent
    313313              float diff = pixelData->data.F32[0] - pixelData->data.F32[1];
    314 #if VARIANCE_FACTORS
     314#ifdef VARIANCE_FACTORS
    315315              float sigma2 = pixelVariances->data.F32[0] * varFactors->data.F32[pixelSources->data.U16[0]] +
    316316                  pixelVariances->data.F32[1] * varFactors->data.F32[pixelSources->data.U16[1]];
    … …  
    346346                  float rej2 = PS_SQR(rej); // Rejection level squared
    347347                  for (int i = 0; i < num; i++) {
    348 #if VARIANCE_FACTORS
     348#ifdef VARIANCE_FACTORS
    349349                      pixelVariances->data.F32[i] *= rej2 * varFactors->data.F32[pixelSources->data.U16[i]];
    350350#else
    … …  
    538538
    539539/// Constructor
    540 pmStackData *pmStackDataAlloc(pmReadout *readout, float weight)
     540pmStackData *pmStackDataAlloc(pmReadout *readout, float weight, float addVariance)
    541541{
    542542    pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return
    … …  
    547547    data->inspect = NULL;
    548548    data->weight = weight;
     549    data->addVariance = addVariance;
    549550
    550551    return data;
    … …  
    602603        }
    603604        varFactors->data.F32[i] = vf;
     605        if (isfinite(data->addVariance)) {
     606            varFactors->data.F32[i] *= data->addVariance;
     607        }
    604608        if (!haveRejects && !data->inspect) {
    605609            data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);

  • trunk/psModules/src/imcombine/pmStack.h

    r17005 r19532  
    88 * @author GLG, MHPCC
    99 *
    10  * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
    11  * @date $Date: 2008-03-17 21:38:43 $
     10 * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
     11 * @date $Date: 2008-09-12 04:12:42 $
    1212 *
    1313 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    3131    psPixels *inspect;                  ///< Pixels to inspect
    3232    float weight;                       ///< Relative weighting for image
     33    float addVariance;                  ///< Additional variance when rejecting
    3334} pmStackData;
    3435
    3536/// Constructor
    3637pmStackData *pmStackDataAlloc(pmReadout *readout, ///< Warped readout (sky cell)
    37                               float weight ///< Weight to apply
     38                              float weight, ///< Weight to apply
     39                              float addVariance ///< Additional variance when rejecting
    3840    );
    3941

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

    r19482 r19532  
    665665
    666666
    667 int pmSubtractionRejectStamps(float *rmsPtr, int *numPtr, pmSubtractionStampList *stamps,
    668                               const psVector *deviations, psImage *subMask, float sigmaRej,
    669                               int footprint)
    670 {
     667int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
     668                              const psVector *deviations, psImage *subMask, float sigmaRej, int footprint)
     669{
     670    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    671671    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
    672672    PS_ASSERT_VECTOR_NON_NULL(deviations, -1);
    … …  
    678678    // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
    679679    // with large N.  Therefore, let's just treat this as a Gaussian distribution.
     680
     681    kernels->mean = NAN;
     682    kernels->rms = NAN;
     683    kernels->numStamps = -1;
    680684
    681685    int numStamps = 0;                  // Number of used stamps
    … …  
    721725    psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
    722726
    723     if (rmsPtr) {
    724         *rmsPtr = rms;
    725     }
    726     if (numPtr) {
    727         *numPtr = numStamps;
    728     }
    729 
    730     psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation from %d stamps: %lf", numStamps, rms);
     727    kernels->mean = mean;
     728    kernels->rms = rms;
     729    kernels->numStamps = numStamps;
     730
     731    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",
     732             numStamps, mean, rms);
    731733
    732734    if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
    733         // User just wanted to calculate and record the RMS for posterity
     735        // User just wanted to calculate and record the deviation for posterity
    734736        return 0;
    735737    }

  • trunk/psModules/src/imcombine/pmSubtraction.h

    r19299 r19532  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2008-08-30 02:28:07 $
     8 * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2008-09-12 04:12:42 $
    1010 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    1111 */
    … …  
    6464
    6565/// Reject stamps
    66 int pmSubtractionRejectStamps(float *rms, ///< RMS deviation, to return
    67                               int *num, ///< Number of good stamps, to return
     66int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, ///< Kernel parameters to update
    6867                              pmSubtractionStampList *stamps, ///< Stamps
    6968                              const psVector *deviations, ///< Deviations for each stamp

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

    r19299 r19532  
    10851085            }
    10861086        }
    1087         deviations->data.F32[i] = sqrtf(devNorm * deviation);
     1087        deviations->data.F32[i] = devNorm * deviation;
    10881088        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    10891089                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);

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

    r19230 r19532  
    3232#define NAME_SOL1 "SOLUTION_1"          // Solution for convolving image 1
    3333#define NAME_SOL2 "SOLUTION_2"          // Solution for convolving image 2
     34#define NAME_MEAN "MEAN"                // Mean of chi^2 from stamps
     35#define NAME_RMS  "RMS"                 // RMS of chi^2 from stamps
     36#define NAME_NUMSTAMPS "NUMSTAMPS"      // Number of good stamps
    3437
    3538#define FALLBACK_EXTNAME "SUBTRACTION_KERNEL" // Fallback extension name
    … …  
    159162            psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL2, 0, "Solution vector 2", kernel->solution2);
    160163        }
     164        psMetadataAddF32(row, PS_LIST_TAIL, NAME_MEAN,  0, "Mean of chi^2 from stamps", kernel->mean);
     165        psMetadataAddF32(row, PS_LIST_TAIL, NAME_RMS,  0, "RMS of chi^2 from stamps", kernel->rms);
     166        psMetadataAddF32(row, PS_LIST_TAIL, NAME_NUMSTAMPS,  0, "Number of good stamps", kernel->numStamps);
    161167    }
    162168    psFree(regions);
    … …  
    167173    // CVS tags, used to identify the version of this file (in case incompatibilities are introduced)
    168174    psString cvsFile = psStringCopy("$RCSfile: pmSubtractionIO.c,v $");
    169     psString cvsRev  = psStringCopy("$Revision: 1.3 $");
    170     psString cvsDate = psStringCopy("$Date: 2008-08-27 02:54:44 $");
     175    psString cvsRev  = psStringCopy("$Revision: 1.4 $");
     176    psString cvsDate = psStringCopy("$Date: 2008-09-12 04:12:42 $");
    171177    psStringSubstitute(&cvsFile, NULL, "RCSfile: ");
    172178    psStringSubstitute(&cvsRev,  NULL, "Revision: ");
    … …  
    338344        TABLE_LOOKUP(int, S32, numRows, NAME_ROWS);
    339345
     346        TABLE_LOOKUP(float, F32, mean,      NAME_MEAN);
     347        TABLE_LOOKUP(float, F32, rms,       NAME_RMS);
     348        TABLE_LOOKUP(int,   S32, numStamps, NAME_NUMSTAMPS);
     349
    340350        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode);
    341351        kernels->numCols = numCols;
    342352        kernels->numRows = numRows;
     353        kernels->mean = mean;
     354        kernels->rms = rms;
     355        kernels->numStamps = numStamps;
    343356
    344357        bool mdok;                      // Status of MD lookup

  • trunk/psModules/src/imcombine/pmSubtractionKernels.h

    r18287 r19532  
    4242    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
    4343    psVector *solution1, *solution2;    ///< Solution for the PSF matching
     44    // Quality information
     45    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
     46    int numStamps;                      ///< Number of good stamps
    4447} pmSubtractionKernels;
    4548

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

    r19482 r19532  
    370370            memCheck("kernels");
    371371
    372             float rmsStamps = NAN;      // RMS for stamps
    373             int numStamps = 0;          // Number of good stamps
    374372            int numRejected = -1;       // Number of rejected stamps in each iteration
    375373            for (int k = 0; k < iter && numRejected != 0; k++) {
    … …  
    407405
    408406                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    409                 numRejected = pmSubtractionRejectStamps(&rmsStamps, &numStamps, stamps, deviations,
    410                                                         subMask, rej, footprint);
     407                numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint);
    411408                if (numRejected < 0) {
    412409                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
    … …  
    430427                    goto MATCH_ERROR;
    431428                }
    432                 pmSubtractionRejectStamps(&rmsStamps, &numStamps, stamps, deviations,
    433                                           subMask, NAN, footprint);
     429                pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint);
    434430                psFree(deviations);
    435431            }
    436432            psFree(stamps);
    437433            stamps = NULL;
    438 
    439             if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    440                 psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM,
    441                                  PS_META_DUPLICATE_OK, "Number of good stamps", numStamps);
    442                 psMetadataAddF32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_RMS,
    443                                  PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps);
    444             }
    445             if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) {
    446                 psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_NUM,
    447                                  PS_META_DUPLICATE_OK, "Number of good stamps", numStamps);
    448                 psMetadataAddF32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS_RMS,
    449                                  PS_META_DUPLICATE_OK, "RMS deviation of stamps", rmsStamps);
    450             }
    451434
    452435            memCheck("solution");

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

    r19299 r19532  
    99#include <pmSubtractionStamps.h>
    1010
    11 #define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata
    12 #define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata
    13 #define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD
    14 #define PM_SUBTRACTION_ANALYSIS_STAMPS_RMS "SUBTRACTION.RMS" // Name of stamp rms in the analysis metadata
    15 #define PM_SUBTRACTION_ANALYSIS_STAMPS_NUM "SUBTRACTION.NUM"// Name of the number of stamps in the analysis MD
    16 #define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Name of variance factor in the analysis MD
     11// Names for things put on the readout analysis metadata
     12#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Kernel used for convolving
     13#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Subtraction mode
     14#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Subtraction region
     15#define PM_SUBTRACTION_ANALYSIS_STAMPS_MEAN "SUBTRACTION.MEAN" // Stamp mean deviation
     16#define PM_SUBTRACTION_ANALYSIS_STAMPS_RMS "SUBTRACTION.RMS" // Stamp RMS deviation
     17#define PM_SUBTRACTION_ANALYSIS_STAMPS_NUM "SUBTRACTION.NUM" // Number of stamps
     18#define PM_SUBTRACTION_ANALYSIS_VARFACTOR "SUBTRACTION.VARFACTOR"// Variance factor
    1719
    1820
Note: See TracChangeset for help on using the changeset viewer.