IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 28, 2010, 5:41:19 PM (16 years ago)
Author:
eugene
Message:

move the generation of the kernel mosaiced image (used for visualization) to a stand-alone function; do not call the visualization function in pmSubtractionAnalysis

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r29543 r29594  
    1414#include "pmSubtractionVisual.h"
    1515
    16 #define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
    17 
    1816//#define TESTING
    1917
     18// save information about the kernel in the output header.  this function also generates the
     19// image normalization, used by ppSubMatchPSFs.c to rescale the output image.
    2020bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header,
    2121                           pmSubtractionKernels *kernels, psRegion *region,
    … …  
    3636        }
    3737
    38         psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    39                          PS_DATA_REGION | PS_META_DUPLICATE_OK,
    40                          "Region over which subtraction was performed", subRegion);
     38        psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_DATA_REGION | PS_META_DUPLICATE_OK, "Region over which subtraction was performed", subRegion);
    4139
    4240        psString string = psRegionToString(*subRegion);
    4341        psFree(subRegion);
    4442
    45         psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK,
    46                          "Region over which subtraction was performed", string);
     43        psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK, "Region over which subtraction was performed", string);
    4744        psFree(string);
    4845    }
    4946
    5047    // Record kernel
    51     psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    52                      PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    53     psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    54                      PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
    55     psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    56                      PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
     48    psMetadataAddPtr(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
     49    psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
     50    psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode);
    5751
    5852    // Realisations of kernel
    59     {
    60         psTrace("psModules.imcombine", 2, "Generating diagnostics...\n");
    61         // Generate image with convolution kernels
    62         int size = kernels->size;       // Half-size of kernel
    63         int fullSize = 2 * size + 1 + 1; // Full size of kernel
    64         int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize;
    65         psImage *convKernels = psImageAlloc((kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) *
    66                                             imageSize - 1 +
    67                                             (kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0),
    68                                             imageSize - 1, PS_TYPE_F32);
    69         psImageInit(convKernels, NAN);
    70         for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
    71             for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) {
    72                 psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
    73                                                            (float)j / (float)KERNEL_MOSAIC,
    74                                                            false); // Image of the kernel
    75                 if (!kernel) {
    76                     psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
    77                     psFree(convKernels);
    78                     return false;
    79                 }
    80 
    81                 if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
    82                                           (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
    83                     psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
    84                     psFree(kernel);
    85                     psFree(convKernels);
    86                     return false;
    87                 }
    88                 psFree(kernel);
    89 
    90                 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    91                     kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
    92                                                       (float)j / (float)KERNEL_MOSAIC,
    93                                                       true); // Image of the kernel
    94                     if (!kernel) {
    95                         psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
    96                         psFree(convKernels);
    97                         return false;
    98                     }
    99 
    100                     if (psImageOverlaySection(convKernels, kernel,
    101                                               (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 4,
    102                                               (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
    103                         psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
    104                         psFree(kernel);
    105                         psFree(convKernels);
    106                         return false;
    107                     }
    108                     psFree(kernel);
    109                 }
    110             }
    111         }
    112 
    113         pmSubtractionVisualPlotConvKernels(convKernels);
    114         psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
    115                            PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels);
    116         psFree(convKernels);
    117     }
     53    psImage *convKernels = pmSubtractionKernelsImageMosaic(kernels);
     54    psMetadataAddImage(analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", PS_META_DUPLICATE_OK, "Realisations of kernel", convKernels);
     55    psFree(convKernels);
    11856
    11957    // sample difference images
    … …  
    166104        int size = kernels->size;       // Half-size of kernel
    167105        int fullSize = 2 * size + 1;    // Full size of kernel
     106
     107        /* in the CENTRAL_DELTA case, the final convolution kernel replaces the central delta
     108           function which was subtracted to force zero flux for each kernel.  for other cases,
     109           we have to include the zero order component of the normalization here.
     110         */
     111
     112        //# if (CENTRAL_DELTA)
     113# if (1)
    168114        float norm = 0.0;               // Normalisation (kernel sum)
     115# else
     116        float norm = 1.0;               // Normalisation (kernel sum)
     117# endif
    169118        for (int y = 0; y < fullSize; y++) {
    170119            for (int x = 0; x < fullSize; x++) {
    … …  
    172121            }
    173122        }
    174         float max = -INFINITY;          // Maximum fraction
     123        psLogMsg("psModules.imcombine", PS_LOG_INFO, "Kernel Integral: %f", norm);
     124
     125        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_DUPLICATE_OK, "Normalisation", norm);
     126        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, PS_META_DUPLICATE_OK, "Normalisation", norm);
     127
     128        float maxCon =  0.0;          // Maximum +fraction > 0.0
     129        float maxDec =  0.0;          // Maximum -fraction < 0.0
    175130        for (int r = 1; r < size; r++) {
    176131            unsigned long r2 = PS_SQR(r); // r^2
    … …  
    186141            }
    187142            float frac = sum / norm;    // Fraction of flux moving towards centre
    188             psTrace("psModules.imcombine", 5, "Deconvolution fraction at %d: %f\n", r, frac);
    189             max = PS_MAX(max, frac);
     143            psTrace("psModules.imcombine", 5, "(De)Convolution fraction at %d: %f\n", r, frac);
     144            maxCon = PS_MAX(maxCon, +frac);
     145            maxDec = PS_MAX(maxDec, -frac);
    190146        }
    191147        psFree(image);
    192148
    193149        {
    194             psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
    195             if (item) {
    196                 max = item->data.F32 = PS_MAX(item->data.F32, max);
    197             } else {
    198                 psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
    199                                  "Maximum deconvolution fraction", max);
    200             }
    201         }
    202 
     150            psMetadataItem *item = NULL;
     151            item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
     152            if (item) {
     153                maxDec = item->data.F32 = PS_MAX(item->data.F32, maxDec);
     154            } else {
     155                psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, "Maximum deconvolution fraction", maxDec);
     156            }
     157            item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
     158            if (item) {
     159                item->data.F32 = maxDec;
     160            } else {
     161                psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, "Maximum deconvolution fraction", maxDec);
     162            }
     163        }
    203164        {
    204             psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
    205             if (item) {
    206                 item->data.F32 = max;
    207             } else {
    208                 psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
    209                                  "Maximum deconvolution fraction", max);
    210             }
    211         }
    212     }
    213 
    214     // Kernel moments
     165            psMetadataItem *item = NULL;
     166            item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); // Previous
     167            if (item) {
     168                maxCon = item->data.F32 = PS_MAX(item->data.F32, maxCon);
     169            } else {
     170                psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX, 0, "Maximum convolution fraction", maxCon);
     171            }
     172            item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); // Previous
     173            if (item) {
     174                item->data.F32 = maxCon;
     175            } else {
     176                psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX, 0, "Maximum convolution fraction", maxCon);
     177            }
     178        }
     179    }
     180
     181    // Kernel moments : since the kernel can be negative, calculate the absolute flux moments,
     182    // eg \Sum(x|f|) / Sum(|f|), \Sum(x^2|f|) / Sum(|f|)
    215183    {
    216184        psImage *image = pmSubtractionKernelImage(kernels, 0.5, 0.5, false); // Image of the kernel
    … …  
    224192        for (int y = 0, v = -size; y < fullSize; y++, v++) {
    225193            for (int x = 0, u = -size; x < fullSize; x++, u++) {
    226                 float value = image->data.F32[y][x]; // Value of kernel
     194                float value = fabs(image->data.F32[y][x]); // Value of kernel
    227195                m00 += value;
    228196                m10 += u * value;
    … …  
    244212        m11 = m11 / m00 - m10 * m01;
    245213
    246         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM,
    247                          PS_META_DUPLICATE_OK, "Normalisation", m00);
    248214        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX,
    249215                         PS_META_DUPLICATE_OK, "Moment in x", m10);
    … …  
    257223                         PS_META_DUPLICATE_OK, "Moment in yy", m02);
    258224
    259         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM,
    260                          PS_META_DUPLICATE_OK, "Normalisation", m00);
    261225        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX,
    262226                         PS_META_DUPLICATE_OK, "Moment in x", m10);
    … …  
    285249    // Quality of fit
    286250    {
    287         psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps",
    288                          kernels->numStamps);
    289         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation",
    290                          kernels->mean);
    291         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
    292                          kernels->rms);
    293 
    294         psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps",
    295                          kernels->numStamps);
    296         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation",
    297                          kernels->mean);
    298         psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
    299                          kernels->rms);
    300 
    301         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
    302         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
    303         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
    304         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
    305         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
    306         psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
    307 
    308         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  0, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
    309         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, 0, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
    310         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  0, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
    311         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, 0, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
    312         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  0, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
    313         psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, 0, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
     251        psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS,   PS_META_REPLACE, "Number of stamps", kernels->numStamps);
     252        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, PS_META_REPLACE, "Mean stamp deviation", kernels->mean);
     253        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS,  PS_META_REPLACE, "RMS stamp deviation", kernels->rms);
     254        psMetadataAddS32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS,   PS_META_REPLACE, "Number of stamps", kernels->numStamps);
     255        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, PS_META_REPLACE, "Mean stamp deviation", kernels->mean);
     256        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS,  PS_META_REPLACE, "RMS stamp deviation", kernels->rms);
     257
     258        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  PS_META_REPLACE, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
     259        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, PS_META_REPLACE, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
     260        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
     261        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
     262        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
     263        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
     264
     265        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_MEAN,  PS_META_REPLACE, "Fractional Sigma of Residuals (Mean)", kernels->fResSigmaMean);
     266        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_SIGMA_STDEV, PS_META_REPLACE, "Fractional Sigma of Residuals (Stdev)", kernels->fResSigmaStdev);
     267        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_MEAN,  PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 2 pix)", kernels->fResOuterMean);
     268        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_OUTER_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 2 pix)", kernels->fResOuterStdev);
     269        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_MEAN,  PS_META_REPLACE, "Fractional Residual Flux (Mean, R > 0 pix)", kernels->fResTotalMean);
     270        psMetadataAddF32(header,   PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FRES_TOTAL_STDEV, PS_META_REPLACE, "Fractional Residual Flux (Stdev, R > 0 pix)", kernels->fResTotalStdev);
    314271    }
    315272
Note: See TracChangeset for help on using the changeset viewer.