IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14709


Ignore:
Timestamp:
Aug 30, 2007, 11:45:26 AM (19 years ago)
Author:
Paul Price
Message:

Cleaning up, consolidating some duplicated code. Fixing photometric preservation for optimum kernels.

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

Legend:

Unmodified
Added
Removed

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

    r14701 r14709  
    44 *  @author GLG, MHPCC
    55 *
    6  *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
    7  *  @date $Date: 2007-08-30 03:50:28 $
     6 *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
     7 *  @date $Date: 2007-08-30 21:45:26 $
    88 *
    99 *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
    … …  
    249249}
    250250
    251 // Generate the convolution given a precalculated kernel
    252 static psKernel *convolvePrecalc(const psKernel *image, // Image to convolve (a kernel for convenience)
    253                                  const psKernel *kernel, // Kernel by which to convolve
    254                                  int footprint // Size of region of interest
    255                                  )
    256 {
    257     psImage *conv = psImageConvolveFFT(image->image, kernel, 0.0); // Convolved image
    258     int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
    259     psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
    260     psFree(conv);
    261     return convolved;
    262 }
    263 
    264251// Generate the convolved reference image
    265252static psKernel *convolveRef(const pmSubtractionKernels *kernels, // Kernel basis functions
    … …  
    308295          if (index < kernels->inner) {
    309296              // Photometric scaling is already built in to the precalculated kernel
    310               return convolvePrecalc(image, kernels->preCalc->data[index], footprint);
     297              return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
    311298          }
    312299          // Using delta function
    … …  
    321308      case PM_SUBTRACTION_KERNEL_ISIS: {
    322309          // Photometric scaling is already built in to the precalculated kernel
    323           return convolvePrecalc(image, kernels->preCalc->data[index], footprint);
     310          return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
    324311      }
    325312      case PM_SUBTRACTION_KERNEL_RINGS: {
    … …  
    355342}
    356343
     344// Convolve an image using FFT
     345static void convolveFFT(psImage *target,// Place the result in here
     346                        const psImage *image, // Image to convolve
     347                        const psKernel *kernel, // Kernel by which to convolve
     348                        psRegion region,// Region of interest
     349                        float background, // Background to add
     350                        int size        // Size of (square) kernel
     351                        )
     352{
     353    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
     354                                  region.y0 - size, region.y1 + size); // Add a border
     355
     356    // Casting away const so psImageSubset can add the child
     357    psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
     358    psImage *convolved = psImageConvolveFFT(subImage, kernel, 0.0); // Convolution
     359    psFree(subImage);
     360
     361    // Now, we have to chop off the borders, and stick it in where it belongs
     362    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size,
     363                                                            size, -size)); // Cut off the edges
     364    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
     365    if (background != 0.0) {
     366        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
     367    } else {
     368        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
     369        for (int y = 0; y < subTarget->numRows; y++) {
     370            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
     371        }
     372    }
     373    psFree(subConv);
     374    psFree(convolved);
     375    psFree(subTarget);
     376    return;
     377}
     378
     379// Convolve an image directly
     380static void convolveDirect(psImage *target, // Put the result here
     381                           const psImage *image, // Image to convolve
     382                           const psKernel *kernel, // Kernel by which to convolve
     383                           psRegion region,// Region of interest
     384                           float background, // Background to add
     385                           int size        // Size of (square) kernel
     386                           )
     387{
     388    for (int y = region.y0; y < region.y1; y++) {
     389        for (int x = region.x0; x < region.x1; x++) {
     390            target->data.F32[y][x] = background;
     391            for (int v = -size; v <= size; v++) {
     392                for (int u = -size; u <= size; u++) {
     393                    target->data.F32[y][x] += kernel->kernel[v][u] *
     394                        image->data.F32[y - v][x - u];
     395                }
     396            }
     397        }
     398    }
     399    return;
     400}
     401
    357402// Mark a pixel as blank in the image, mask and weight
    358403static inline void markBlank(psImage *image, // Image to mark as blank
    … …  
    372417    return;
    373418}
     419
     420//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     421// Semi-public functions
     422//////////////////////////////////////////////////////////////////////////////////////////////////////////////
     423
     424// Generate the convolution given a precalculated kernel
     425psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel)
     426{
     427    PS_ASSERT_KERNEL_NON_NULL(image, NULL);
     428    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
     429
     430    psImage *conv = psImageConvolveFFT(image->image, kernel, 0.0); // Convolved image
     431    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
     432    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
     433    psFree(conv);
     434    return convolved;
     435}
     436
    374437
    375438//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    … …  
    916979}
    917980
    918 // Convolve an image using FFT
    919 static void convolveFFT(psImage *target,// Place the result in here
    920                         const psImage *image, // Image to convolve
    921                         const psKernel *kernel, // Kernel by which to convolve
    922                         psRegion region,// Region of interest
    923                         float background, // Background to add
    924                         int size        // Size of (square) kernel
    925                         )
    926 {
    927     psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
    928                                   region.y0 - size, region.y1 + size); // Add a border
    929 
    930     // Casting away const so psImageSubset can add the child
    931     psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
    932     psImage *convolved = psImageConvolveFFT(subImage, kernel, 0.0); // Convolution
    933     psFree(subImage);
    934 
    935     // Now, we have to chop off the borders, and stick it in where it belongs
    936     psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size,
    937                                                             size, -size)); // Cut off the edges
    938     psImage *subTarget = psImageSubset(target, region); // Target for this subregion
    939     if (background != 0.0) {
    940         psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
    941     } else {
    942         int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
    943         for (int y = 0; y < subTarget->numRows; y++) {
    944             memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
    945         }
    946     }
    947     psFree(subConv);
    948     psFree(convolved);
    949     psFree(subTarget);
    950     return;
    951 }
    952 
    953 // Convolve an image directly
    954 static void convolveDirect(psImage *target, // Put the result here
    955                            const psImage *image, // Image to convolve
    956                            const psKernel *kernel, // Kernel by which to convolve
    957                            psRegion region,// Region of interest
    958                            float background, // Background to add
    959                            int size        // Size of (square) kernel
    960                            )
    961 {
    962     for (int y = region.y0; y < region.y1; y++) {
    963         for (int x = region.x0; x < region.x1; x++) {
    964             target->data.F32[y][x] = background;
    965             for (int v = -size; v <= size; v++) {
    966                 for (int u = -size; u <= size; u++) {
    967                     target->data.F32[y][x] += kernel->kernel[v][u] *
    968                         image->data.F32[y - v][x - u];
    969                 }
    970             }
    971         }
    972     }
    973     return;
    974 }
    975 
    976981bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
    977982                           const psImage *inImage, const psImage *inWeight, const psImage *subMask,

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

    r14701 r14709  
    66 * @author GLG, MHPCC
    77 *
    8  * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
    9  * @date $Date: 2007-08-30 03:50:28 $
     8 * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
     9 * @date $Date: 2007-08-30 21:45:26 $
    1010 *
    1111 * Copyright 2004-207 Institute for Astronomy, University of Hawaii
    … …  
    9393    );
    9494
     95/// Generate the convolution of an image, given a precalculated kernel
     96///
     97/// The 'image' is a kernel for convenience --- intended to be a stamp
     98psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, ///< Image to convolve
     99                                              const psKernel *kernel ///< Kernel by which to convolve
     100    );
     101
    95102/// @}
    96103#endif

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

    r14701 r14709  
    216216
    217217    // Subtract a reference kernel from all the others to maintain flux scaling
    218 //    if (spatialOrder > 0) {
     218    if (spatialOrder > 0) {
    219219        int subIndex = 0;                   // Index of kernel to subtract from others
    220220        psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others
    … …  
    230230            }
    231231        }
    232 //    }
     232    }
    233233
    234234    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
    … …  
    454454    // Subtract unity from the kernels to maintain photometric flux scaling
    455455    int numGaussians = fwhms->n;       // Number of Gaussians
    456 //    if (spatialOrder > 0) {
     456    if (spatialOrder > 0) {
    457457        for (int i = 0, index = 0; i < numGaussians; i++) {
    458458            for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
    … …  
    465465            }
    466466        }
    467 //    }
     467    }
    468468
    469469    int numISIS = kernels->num;         // Number of ISIS kernels

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

    r14671 r14709  
    99
    1010#include "pmSubtractionStamps.h"
     11#include "pmSubtraction.h"
     12#include "pmSubtractionParams.h"
    1113
    1214//////////////////////////////////////////////////////////////////////////////////////////////////////////////
    1315
     16#if 0
    1417// Convolve the reference stamp by the kernel
    1518static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved
    … …  
    4548    return convolution;
    4649}
    47 
     50#endif
    4851
    4952// Accumulate cross-term sums for a stamp
    … …  
    5356                            const pmSubtractionStamp *stamp, // Stamp with weight
    5457                            const psKernel *input, // Input image, I(x)
    55                             const psKernel *convolution // Convolution, conv(x)
    56     )
    57 {
    58     // Range of convolution
    59     int xMin = convolution->xMin;
    60     int xMax = convolution->xMax;
    61     int yMin = convolution->yMin;
    62     int yMax = convolution->yMax;
    63 
     58                            int kernelIndex, // Index for kernel component
     59                            int footprint // Size of region of interest
     60    )
     61{
    6462    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
    65 
    66     for (int y = yMin; y <= yMax; y++) {
    67         psF32 *in = &input->kernel[y][xMin]; // Dereference input
    68         psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
    69         psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
    70         for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) {
     63    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     64
     65    for (int y = -footprint; y <= footprint; y++) {
     66        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     67        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     68        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
     69        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    7170            double temp = *in / *wt; // Temporary product
    7271            *sumI += temp;
    … …  
    8281                                   double *sumCC, // Sum of conv(x)^2/sigma(x)^2
    8382                                   const pmSubtractionStamp *stamp, // Stamp with input and weight
    84                                    const psKernel *convolution // Convolution, conv(x)
    85     )
    86 {
    87     // Range of convolution
    88     int xMin = convolution->xMin;
    89     int xMax = convolution->xMax;
    90     int yMin = convolution->yMin;
    91     int yMax = convolution->yMax;
    92 
     83                                   int kernelIndex, // Index for kernel component
     84                                   int footprint // Size of region of interest
     85    )
     86{
    9387    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
    94 
    95     for (int y = yMin; y <= yMax; y++) {
    96         psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
    97         psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
    98         for (int x = xMin; x <= xMax; x++, wt++, conv++) {
     88    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     89
     90    for (int y = -footprint; y <= footprint; y++) {
     91        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     92        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
     93        for (int x = -footprint; x <= footprint; x++, wt++, conv++) {
    9994            double convNoise = *conv / *wt; // Temporary product
    10095            *sumC += convNoise;
    … …  
    106101
    107102static double accumulateChi2(psKernel *input, // Input stamp
    108                              psKernel *convolution, // Convolution, R(x)*k(u)
    109103                             pmSubtractionStamp *stamp, // Stamp with weight
     104                             int kernelIndex, // Index for kernel component
    110105                             double coeff, // Coefficient of convolution
    111                              double bg  // Background term
    112     )
    113 {
    114     // Range of convolution
    115     int xMin = convolution->xMin;
    116     int xMax = convolution->xMax;
    117     int yMin = convolution->yMin;
    118     int yMax = convolution->yMax;
    119 
     106                             double bg,  // Background term
     107                             int footprint // Size of region of interest
     108    )
     109{
    120110    double chi2 = 0.0;
    121     psKernel *weight = stamp->weight;
    122 
    123     for (int y = yMin; y <= yMax; y++) {
    124         psF32 *in = &input->kernel[y][xMin]; // Dereference input
    125         psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
    126         psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
    127         for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) {
     111    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
     112    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     113
     114    for (int y = -footprint; y <= footprint; y++) {
     115        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     116        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
     117        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
     118        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
    128119            chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt;
    129120        }
    … …  
    159150static void subtractConvolution(psKernel *input, // Input stamp
    160151                                const pmSubtractionStamp *stamp, // Stamp with weight
    161                                 const psKernel *convolution, // Convolution, R(x)*k(u)
     152                                int kernelIndex, // Index for kernel component
    162153                                float coeff, // Coefficient of subtraction
    163                                 float bg // Background term
    164     )
    165 {
    166     // Range of convolution
    167     int xMin = convolution->xMin;
    168     int xMax = convolution->xMax;
    169     int yMin = convolution->yMin;
    170     int yMax = convolution->yMax;
    171 
    172     psKernel *weight = stamp->weight;   // Weight map
    173 
    174     for (int y = yMin; y <= yMax; y++) {
    175         psF32 *in = &input->kernel[y][xMin]; // Dereference input
    176         psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
    177         psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
    178         for (int x = xMin; x <= xMax; x++, in++, conv++, wt++) {
     154                                float bg, // Background term
     155                                int footprint // Size of region of interest
     156    )
     157{
     158    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
     159
     160    for (int y = -footprint; y <= footprint; y++) {
     161        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
     162        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
     163        for (int x = -footprint; x <= footprint; x++, in++, conv++) {
    179164            *in -= *conv * coeff + bg;
    180165        }
    … …  
    229214        psKernel *input = stamp->input; // Input image of interest
    230215        psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image
    231         inputs->data[i] = psKernelAllocFromImage(copy, size + footprint + 1, size + footprint + 1);
     216        inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
    232217        psFree(copy);                   // Drop reference
    233218    }
    234219
    235220    // Generate the convolutions
    236     psArray *convolutions = psArrayAlloc(numKernels); // Convolutions per kernel component and stamp
    237     for (int i = 0; i < numKernels; i++) {
    238         psKernel *kernel = kernels->preCalc->data[i]; // Precalculated kernel
    239 
    240         psArray *stampConv = psArrayAlloc(numStamps); // Convolutions for each stamp
    241         convolutions->data[i] = stampConv;
    242 
    243         for (int j = 0; j < numStamps; j++) {
    244             stampConv->data[j] = convolveStamp(stamps->data[j], kernel, footprint);
     221    for (int i = 0; i < numStamps; i++) {
     222        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
     223        if (!stamp->convolutions) {
     224            stamp->convolutions = psArrayAlloc(numKernels);
     225        }
     226
     227        for (int j = 0; j < numKernels; j++) {
     228            psKernel *kernel = kernels->preCalc->data[j]; // Precalculated kernel
     229            stamp->convolutions->data[j] = p_pmSubtractionConvolveStampPrecalc(stamp->reference, kernel);
    245230        }
    246231    }
    … …  
    264249    psVectorInit(sumCC, 0.0);
    265250    for (int i = 0; i < numKernels; i++) {
    266         psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp
    267251        for (int j = 0; j < numStamps; j++) {
    268             accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j],
    269                                    stampsConv->data[j]);
     252            accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], i, footprint);
    270253        }
    271254    }
    … …  
    298281            double sumIC = 0.0;         // sum of I(x)C(x)/sigma(x)^2
    299282
    300             psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp
    301 
    302283            for (int j = 0; j < numStamps; j++) {
    303                 accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], stampsConv->data[j]);
     284                accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], i, footprint);
    304285            }
    305286
    … …  
    310291            double chi2 = 0.0;          // Chi^2
    311292            for (int j = 0; j < numStamps; j++) {
    312                 chi2 += accumulateChi2(inputs->data[j], stampsConv->data[j], stamps->data[j], coeff, bg);
     293                chi2 += accumulateChi2(inputs->data[j], stamps->data[j], i, coeff, bg, footprint);
    313294            }
    314295
    … …  
    329310        ranking->data.S32[bestIndex] = iter;
    330311        // Remove its contribution, and don't include it in the future.
    331         psArray *stampsConv = convolutions->data[bestIndex]; // Convolutions per stamp
    332312        for (int j = 0; j < numStamps; j++) {
    333             subtractConvolution(inputs->data[j], stamps->data[j], stampsConv->data[j],
    334                                            bestCoeff, bestBG);
     313            subtractConvolution(inputs->data[j], stamps->data[j], bestIndex, bestCoeff, bestBG, footprint);
    335314        }
    336315
    … …  
    349328    }
    350329    psFree(inputs);
    351     psFree(convolutions);
    352330    psFree(sumC);
    353331    psFree(sumCC);
    … …  
    392370        psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
    393371        for (int i = 1; i < newSize; i++) {
    394             psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    395             psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
     372            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
     373                psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
     374                psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
     375            }
    396376        }
    397377    } else if (type == PM_SUBTRACTION_KERNEL_GUNK) {
    … …  
    400380
    401381        for (int i = 0; i < newSize; i++) {
    402             psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
    403             kernel->kernel[0][0] -= 1.0;
     382            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
     383                psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
     384                kernel->kernel[0][0] -= 1.0;
     385            }
    404386        }
    405387
Note: See TracChangeset for help on using the changeset viewer.