IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26717


Ignore:
Timestamp:
Jan 28, 2010, 5:32:15 PM (16 years ago)
Author:
Paul Price
Message:

Attempting to use a wide function for the normalisation component, rather than delta-function.

Location:
branches/pap/psModules/src/imcombine
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/psModules/src/imcombine/pmSubtraction.c

    r26667 r26717  
    188188    if (normalise) {
    189189        // Put in the normalisation component
    190         kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
     190//        kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
     191        double value = wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels);
     192        for (int v = -size; v <= size; v++) {
     193            for (int u = -size; u <= size; u++) {
     194                kernel->kernel[v][u] += value * kernels->norm->kernel[v][u];
     195            }
     196        }
    191197    }
    192198
    … …  
    754760      case PM_SUBTRACTION_MODE_1:
    755761        stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint);
     762        stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image1, kernels->norm);
    756763        break;
    757764      case PM_SUBTRACTION_MODE_2:
    758765        stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint);
     766        stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image2, kernels->norm);
    759767        break;
    760768      case PM_SUBTRACTION_MODE_UNSURE:
    … …  
    762770        stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint);
    763771        stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint);
     772        stamp->normConv = p_pmSubtractionConvolveStampPrecalc(stamp->image1, kernels->norm);
    764773        break;
    765774      default:
    … …  
    908917                psFree(stamp->vector);
    909918                stamp->vector = NULL;
     919
     920                psFree(stamp->normConv);
     921                stamp->normConv = NULL;
    910922            } else {
    911923                numGood++;

  • branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

    r26703 r26717  
    229229                                      const psArray *convolutions1, // Convolutions of image 1 for each kernel
    230230                                      const psArray *convolutions2, // Convolutions of image 2 for each kernel
     231                                      const psKernel *normConv,     // Normalisation, convolved
    231232                                      const pmSubtractionKernels *kernels, // Kernels
    232233                                      const psImage *polyValues, // Spatial polynomial values
    … …  
    370371        double sumB = 0.0;              // Sum of B products (for matrix, background)
    371372        double sumI2 = 0.0;             // Sum of I_2 (for vector, background)
     373
     374        double sumAN = 0.0;             // Matrix: A*norm
     375        double sumBN = 0.0;             // Matrix: B*norm
     376
    372377        for (int y = - footprint; y <= footprint; y++) {
    373378            for (int x = - footprint; x <= footprint; x++) {
    374379                double a = iConv1->kernel[y][x];
    375380                double b = iConv2->kernel[y][x];
     381                double n = normConv->kernel[y][x];
    376382                float i1 = image1->kernel[y][x];
    377383                float i2 = image2->kernel[y][x];
    … …  
    381387                double ai1 = a * i1;
    382388                double bi1 = b * i1;
     389
     390                double an = a * n;
     391                double bn = b * n;
    383392
    384393                if (weight) {
    … …  
    391400                    b *= wtVal;
    392401                    i2 *= wtVal;
     402
     403                    an *= wtVal;
     404                    bn *= wtVal;
    393405                }
    394406                if (window) {
    … …  
    401413                    b *= wtVal;
    402414                    i2 *= wtVal;
     415
     416                    an *= wtVal;
     417                    bn *= wtVal;
    403418                }
    404419                sumAI2 += ai2;
    … …  
    409424                sumB += b;
    410425                sumI2 += i2;
     426
     427                sumAN += an;
     428                sumBN += bn;
    411429            }
    412430        }
    … …  
    420438            double b   = sumB * poly[iTerm];
    421439
     440            double an  = sumAN * poly[iTerm];
     441            double bn  = sumBN * poly[iTerm];
     442
    422443            if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    423                 matrix->data.F64[iIndex][normIndex] = ai1;
    424                 matrix->data.F64[normIndex][iIndex] = ai1;
    425                 matrix->data.F64[iIndex + numParams][normIndex] = bi1;
    426                 matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     444//                matrix->data.F64[iIndex][normIndex] = ai1;
     445//                matrix->data.F64[normIndex][iIndex] = ai1;
     446//                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
     447//                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
     448                matrix->data.F64[iIndex][normIndex] = an;
     449                matrix->data.F64[normIndex][iIndex] = an;
     450                matrix->data.F64[iIndex + numParams][normIndex] = bn;
     451                matrix->data.F64[normIndex][iIndex + numParams] = bn;
    427452            }
    428453            if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
    … …  
    453478    double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
    454479    double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
     480
     481    double sumN = 0.0;              // Matrix: bg*norm
     482    double sumN2 = 0.0;             // Matrix: norm*norm
     483    double sumNI2 = 0.0;            // Vector: I2*norm
     484
    455485    for (int y = - footprint; y <= footprint; y++) {
    456486        for (int x = - footprint; x <= footprint; x++) {
    … …  
    461491            double one = 1.0;
    462492            double i1i2 = i1 * i2;
     493
     494            double n = normConv->kernel[y][x];
     495            double n2 = n * n;
     496            double ni2 = n * i2;
    463497
    464498            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
    … …  
    474508                i2 *= wtVal;
    475509                i1i2 *= wtVal;
     510
     511                n *= wtVal;
     512                n2 *= wtVal;
     513                ni2 *= wtVal;
    476514            }
    477515            if (window) {
    … …  
    482520                i2 *= wtVal;
    483521                i1i2 *= wtVal;
     522
     523                n *= wtVal;
     524                n2 *= wtVal;
     525                ni2 *= wtVal;
    484526            }
    485527            sumI1 += i1;
    … …  
    488530            sumI2 += i2;
    489531            sumI1I2 += i1i2;
     532
     533            sumN += n;
     534            sumN2 += n2;
     535            sumNI2 += ni2;
    490536        }
    491537    }
    … …  
    496542
    497543    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
    498         matrix->data.F64[normIndex][normIndex] = sumI1I1;
    499         vector->data.F64[normIndex] = sumI1I2;
     544//        matrix->data.F64[normIndex][normIndex] = sumI1I1;
     545//        vector->data.F64[normIndex] = sumI1I2;
     546        matrix->data.F64[normIndex][normIndex] = sumN2;
     547        vector->data.F64[normIndex] = sumNI2;
     548
    500549    }
    501550    if (mode & PM_SUBTRACTION_EQUATION_BG) {
    … …  
    504553    }
    505554    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
    506         matrix->data.F64[bgIndex][normIndex] = sumI1;
    507         matrix->data.F64[normIndex][bgIndex] = sumI1;
     555//        matrix->data.F64[bgIndex][normIndex] = sumI1;
     556//        matrix->data.F64[normIndex][bgIndex] = sumI1;
     557        matrix->data.F64[bgIndex][normIndex] = sumN;
     558        matrix->data.F64[normIndex][bgIndex] = sumN;
    508559    }
    509560    return true;
    … …  
    729780                                           stamp->image1, stamp->image2,
    730781                                           weight, window, stamp->convolutions1, stamp->convolutions2,
    731                                            kernels, polyValues, footprint, stamps->normWindow, mode);
     782                                           stamp->normConv, kernels, polyValues, footprint,
     783                                           stamps->normWindow, mode);
    732784        break;
    733785      default:

  • branches/pap/psModules/src/imcombine/pmSubtractionKernels.c

    r26686 r26717  
    1515
    1616#define RINGS_BUFFER 10                 // Buffer size for RINGS data
     17
     18#define NORM_WIDTH 7.5
    1719
    1820// Free function for pmSubtractionKernels
    … …  
    3032    psFree(kernels->solution2);
    3133    psFree(kernels->sampleStamps);
     34    psFree(kernels->norm);
    3235}
    3336
    … …  
    242245
    243246    if (zeroNull) {
    244         preCalc->kernel->kernel[0][0] -= 1.0;
     247        //        preCalc->kernel->kernel[0][0] -= 1.0;
     248        for (int v = -size; v <= size; v++) {
     249            for (int u = -size; u <= size; u++) {
     250                preCalc->kernel->kernel[v][u] -= kernels->norm->kernel[v][u];
     251            }
     252        }
    245253    }
    246254
    … …  
    582590    kernels->fMinResMean  = NAN;
    583591    kernels->fMinResStdev = NAN;
     592
     593    kernels->norm = psKernelAlloc(-size, size, -size, size);
     594    double sumNorm = 0.0;
     595    float sigma = NORM_WIDTH / 2.35;
     596    for (int y = -size; y <= size; y++) {
     597        for (int x = -size; x <= size; x++) {
     598            float z = 0.5 * (PS_SQR(x) + PS_SQR(y)) / PS_SQR(sigma);
     599            sumNorm += kernels->norm->kernel[y][x] = expf(-z);
     600        }
     601    }
     602    for (int y = -size; y <= size; y++) {
     603        for (int x = -size; x <= size; x++) {
     604            kernels->norm->kernel[y][x] /= sumNorm;
     605        }
     606    }
    584607
    585608    return kernels;

  • branches/pap/psModules/src/imcombine/pmSubtractionKernels.h

    r26593 r26717  
    1212    PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
    1313    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
    14     PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
     14    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
    1515    PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
    1616    PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
    … …  
    3737    psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
    3838    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM)
     39    psKernel *norm;                     ///< Normalisation component
    3940    float penalty;                      ///< Penalty for wideness
    4041    psVector *penalties;                ///< Penalty for each kernel component
    … …  
    4950    float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
    5051    int numStamps;                      ///< Number of good stamps
    51     float fSigResMean;                  ///< mean fractional stdev of residuals
    52     float fSigResStdev;                 ///< stdev of fractional stdev of residuals
    53     float fMaxResMean;                  ///< mean fractional positive swing in residuals
    54     float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
    55     float fMinResMean;                  ///< mean fractional negative swing in residuals
    56     float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
    57     psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
     52    float fSigResMean;                  ///< mean fractional stdev of residuals
     53    float fSigResStdev;                 ///< stdev of fractional stdev of residuals
     54    float fMaxResMean;                  ///< mean fractional positive swing in residuals
     55    float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
     56    float fMinResMean;                  ///< mean fractional negative swing in residuals
     57    float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
     58    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
    5859} pmSubtractionKernels;
    5960
    6061// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
    6162typedef struct {
    62     psVector *uCoords;                  // used by RINGS
    63     psVector *vCoords;                  // used by RINGS
    64     psVector *poly;                     // used by RINGS
    65 
    66     psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
    67     psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
    68     psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
     63    psVector *uCoords;                  // used by RINGS
     64    psVector *vCoords;                  // used by RINGS
     65    psVector *poly;                     // used by RINGS
     66
     67    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
     68    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
     69    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
    6970} pmSubtractionKernelPreCalc;
    7071
    … …  
    168169pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(
    169170    pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated)
    170     int uOrder,                    ///< order in x-direction
    171     int vOrder,                    ///< order in x-direction
    172     int size,                      ///< Half-size of the kernel
    173     float sigma                    ///< sigma of gaussian kernel
     171    int uOrder,                    ///< order in x-direction
     172    int vOrder,                    ///< order in x-direction
     173    int size,                      ///< Half-size of the kernel
     174    float sigma                    ///< sigma of gaussian kernel
    174175    );
    175176
    … …  
    202203/// Generate ISIS + RADIAL_HERM kernels
    203204pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, ///< Half-size of the kernel
    204                                                       int spatialOrder, ///< Order of spatial variations
    205                                                       const psVector *fwhms, ///< Gaussian FWHMs
    206                                                       const psVector *orders, ///< Polynomial order of gaussians
    207                                                       float penalty, ///< Penalty for wideness
    208                                                       pmSubtractionMode mode ///< Mode for subtraction
     205                                                      int spatialOrder, ///< Order of spatial variations
     206                                                      const psVector *fwhms, ///< Gaussian FWHMs
     207                                                      const psVector *orders, ///< Polynomial order of gaussians
     208                                                      float penalty, ///< Penalty for wideness
     209                                                      pmSubtractionMode mode ///< Mode for subtraction
    209210                                               );
    210211
    … …  
    220221/// Generate DECONV_HERM kernels
    221222pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel
    222                                                      int spatialOrder, ///< Order of spatial variations
    223                                                      const psVector *fwhms, ///< Gaussian FWHMs
    224                                                      const psVector *orders, ///< order of hermitian polynomials
    225                                                      float penalty, ///< Penalty for wideness
    226                                                      pmSubtractionMode mode ///< Mode for subtraction
     223                                                     int spatialOrder, ///< Order of spatial variations
     224                                                     const psVector *fwhms, ///< Gaussian FWHMs
     225                                                     const psVector *orders, ///< order of hermitian polynomials
     226                                                     float penalty, ///< Penalty for wideness
     227                                                     pmSubtractionMode mode ///< Mode for subtraction
    227228    );
    228229

  • branches/pap/psModules/src/imcombine/pmSubtractionStamps.c

    r26703 r26717  
    312312    stamp->vector = NULL;
    313313    stamp->norm = NAN;
     314    stamp->normConv = NULL;
    314315
    315316    return stamp;

  • branches/pap/psModules/src/imcombine/pmSubtractionStamps.h

    r26703 r26717  
    8383    psVector *vector;                   ///< Least-squares vector, or NULL
    8484    double norm;                        ///< Normalisation difference
     85    psKernel *normConv;                 ///< Convolution of normalisation
    8586    pmSubtractionStampStatus status;    ///< Status of stamp
    8687} pmSubtractionStamp;
Note: See TracChangeset for help on using the changeset viewer.