IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.