IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27933


Ignore:
Timestamp:
May 12, 2010, 12:56:07 PM (16 years ago)
Author:
Paul Price
Message:

Implementing better Jacobian correction. Previously, was calculating Jacobian from the projections (tangent plane --> sky), which does not account for scales in the transformations from the pixels to the tangent plane. Now using the local linear mapping to provide the Jacobian, so that it's applied on small areas. Also using the (square root of the) Jacobian to fix the plate scale of the covariance (instead of a single value for the entire focal plane) so the noise model should be more accurate as well.

Location:
branches/pap/pswarp/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/pswarp/src/pswarp.h

    r27126 r27933  
    2929#define PSASTRO_RECIPE "PSASTRO" ///< Name of the recipe to use
    3030
    31 #define PSWARP_ANALYSIS_VARFACTOR "PSWARP.VARFACTOR" ///< Name for variance factor in analysis metadata
    3231#define PSWARP_ANALYSIS_GOODPIX   "PSWARP.GOODPIX" ///< Name for number of good pixels in analysis metadata
    3332#define PSWARP_ANALYSIS_COVARIANCES "PSWARP.COVARIANCES" ///< Name for covariance matrices on analysis MD
     33#define PSWARP_ANALYSIS_JACOBIAN "PSWARP.JACOBIAN" ///< Name for Jacobian on analysis MD
    3434
    3535/**
    … …  
    7070    int xMin, xMax, yMin, yMax;         ///< Bounds of tile
    7171    psKernel *covariance;               ///< Covariance matrix
     72    double jacobian;                    ///< (Square root of) local Jacobian
    7273} pswarpTransformTileArgs;
    7374

  • branches/pap/pswarp/src/pswarpLoop.c

    r27733 r27933  
    273273        psAssert(covariances, "Should be there");
    274274        psArray *covars = psListToArray(covariances); // Array of covariance matrices
    275         output->covariance = psImageCovarianceAverage(covars);
     275        psKernel *covar = psImageCovarianceAverage(covars);
    276276        psFree(covars);
    277277        psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
     278
     279        // Correct covariance matrix scale for the mean (square root of the) Jacobian
     280        double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
     281        int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
     282        jacobian /= goodPixels;
     283        output->covariance = psImageCovarianceScale(covar, jacobian);
     284        psFree(covar);
     285
    278286        psImageCovarianceTransfer(output->variance, output->covariance);
    279     }
    280 
    281     // Correct image for change in the plate scale
    282     {
    283         psAssert(input && input->fpa && input->fpa->toSky, "Require astrometry for input");
    284         psAssert(outFPA && outFPA && outFPA->toSky, "Require astrometry for output");
    285 
    286         double inScale = input->fpa->toSky->Xs + input->fpa->toSky->Ys; // Plate scale for input
    287         double outScale = outFPA->toSky->Xs + outFPA->toSky->Ys; // Plate scale for output
    288         float correction = PS_SQR(outScale / inScale); // Correction factor to apply to image
    289         psLogMsg("pswarp", PS_LOG_INFO, "Correcting flux by %f to account for pixel scales", correction);
    290         psBinaryOp(output->image, output->image, "*", psScalarAlloc(correction, PS_TYPE_F32));
    291         if (output->variance) {
    292             psBinaryOp(output->variance, output->variance, "*",
    293                        psScalarAlloc(PS_SQR(correction), PS_TYPE_F32));
    294         }
    295         psKernel *covar = psImageCovarianceScale(output->covariance, outScale / inScale); // Scaled covariance
    296         psFree(output->covariance);
    297         output->covariance = covar;
    298287    }
    299288

  • branches/pap/pswarp/src/pswarpTransformReadout.c

    r27096 r27933  
    148148    psThreadJob *job = NULL;
    149149    int xMin = output->image->numCols, xMax = 0, yMin = output->image->numRows, yMax = 0; // Bounds
    150     int goodPixels = 0;                 // total number of good pixels across all tiles
     150    int goodPixels = psMetadataLookupS32(&mdok, output->analysis, PSWARP_ANALYSIS_GOODPIX); // Number of pixels
    151151    psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
    152152                                              PSWARP_ANALYSIS_COVARIANCES); // Collection of covar. matrices
    … …  
    157157        psFree(covariances);            // Drop reference; still have the copy on the analysis metadata
    158158    }
     159    double jacobian = psMetadataLookupF64(&mdok, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
     160    if (!isfinite(jacobian)) {
     161        jacobian = 0.0;
     162    }
     163
    159164    while ((job = psThreadJobGetDone()) != NULL) {
    160165        if (job->args->n < 1) {
    … …  
    171176                psListAdd(covariances, PS_LIST_TAIL, args->covariance);
    172177            }
     178            if (args->goodPixels > 0 && isfinite(args->jacobian)) {
     179                jacobian += args->jacobian * args->goodPixels;
     180            }
    173181        }
    174182        psFree(job);
    175183    }
    176184    psFree(grid);
     185
     186    psMetadataAddS32(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_GOODPIX, PS_META_REPLACE,
     187                     "Number of good pixels", goodPixels);
     188    psMetadataAddF64(output->analysis, PS_LIST_TAIL, PSWARP_ANALYSIS_JACOBIAN, PS_META_REPLACE,
     189                     "Jacobian of transformation", jacobian);
    177190
    178191    if (xMin < xMax && yMin < yMax) {

  • branches/pap/pswarp/src/pswarpTransformTile.c

    r27096 r27933  
    4444    args->yMax = PS_MIN_S32;
    4545    args->covariance = NULL;
     46    args->jacobian = NAN;
    4647
    4748    return args;
    … …  
    5051bool pswarpTransformTile(pswarpTransformTileArgs *args)
    5152{
    52     psImage *inImage = args->input->image; ///< Input image
    53     psImage *outImage = args->output->image; ///< Output image
     53    pmReadout *input = args->input;        // Input readout
     54    psImage *inImage = input->image, *inMask = input->mask; // Input images
     55    pmReadout *output = args->output;      // Output readout
     56    psImage *outImage = output->image, *outVariance = output->variance, *outMask = output->mask; // Outputs
    5457
    5558    int inNumCols = inImage->numCols, inNumRows = inImage->numRows; ///< Size of input image
    … …  
    6164
    6265    // Dereference images for convenience
    63     psF32 **outImageData     = args->output->image->data.F32;
    64     psF32 **outVarData       = (args->output->variance) ? args->output->variance->data.F32 : NULL;
    65     psImageMaskType **outMaskData = (args->output->mask)   ? args->output->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
    66     psImageMaskType **inMaskData  = (args->input->mask)    ? args->input->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
     66    psF32 **outImageData     = outImage->data.F32;
     67    psF32 **outVarData       = outVariance ? outVariance->data.F32 : NULL;
     68    psImageMaskType **outMaskData = outMask ? outMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
     69    psImageMaskType **inMaskData  = inMask ? inMask->data.PS_TYPE_IMAGE_MASK_DATA : NULL;
    6770
    6871    pswarpMap *map = args->grid->maps[args->gridX][args->gridY]; ///< Map for this tile
    … …  
    7477    int yMin = PS_MAX(minPt.y, 0);
    7578    int yMax = PS_MIN(maxPt.y, outNumRows);
     79
     80    double jacobian = map->Xx * map->Yy - map->Yx * map->Xy; // Jacobian of transformation
     81    double jacobian2 = PS_SQR(jacobian);                     // Square Jacobian
    7682
    7783    // Iterate over the output image pixels (parent frame)
    … …  
    8793            // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
    8894            // both are in the parent frames of the input and output images.
    89             double xInRaw, yInRaw;      // Input raw pixel coordinates
    90             pswarpMapApply(&xInRaw, &yInRaw, map, x + 0.5, y + 0.5);
    91             double xIn = xInRaw - outCol0, yIn = yInRaw - outRow0; // Position on input image
     95            double xIn, yIn;            // Input pixel coordinates
     96            pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
    9297            if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
    9398                continue;
    … …  
    105110
    106111            if (outImageData) {
    107                 outImageData[yOut][xOut] = imageValue;
     112                outImageData[yOut][xOut] = imageValue * jacobian;
    108113            }
    109114            if (outVarData) {
    110                 outVarData[yOut][xOut] = varValue;
     115                outVarData[yOut][xOut] = varValue * jacobian2;
    111116            }
    112117            if (outMaskData) {
    … …  
    122127        double xIn, yIn;                // Position of interest on input
    123128        pswarpMapApply(&xIn, &yIn, map, xOut + 0.5, yOut + 0.5);
    124         // XXX Why are we subtracting the *output* col0,row0 from the *input* coordinates?
    125         // I expect these are zero, so that it makes no difference, but it's distracting.
    126         xIn -= outCol0;
    127         yIn -= outRow0;
    128129        psKernel *kernel = psImageInterpolationKernel(xIn, yIn, args->interp->mode); // Interpolation kernel
    129130        args->covariance = psImageCovarianceCalculate(kernel, args->input->covariance);
     131        args->jacobian = sqrt(jacobian);
    130132        psFree(kernel);
    131133    }
Note: See TracChangeset for help on using the changeset viewer.