IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:45:22 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simmosaic_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/simmosaic_branches

  • branches/simmosaic_branches/psModules

  • branches/simmosaic_branches/psModules/src/imcombine/pmSubtractionAnalysis.c

    r23780 r27839  
    1616#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
    1717
    18 
    19 bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region,
     18//#define TESTING
     19
     20bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header,
     21                           pmSubtractionKernels *kernels, psRegion *region,
    2022                           int numCols, int numRows)
    2123{
    22     if (analysis) {
    23         PS_ASSERT_METADATA_NON_NULL(analysis, false);
    24     }
     24    PS_ASSERT_METADATA_NON_NULL(analysis, false);
     25    PS_ASSERT_METADATA_NON_NULL(header, false);
    2526    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    2627    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false);
    … …  
    3839                         PS_DATA_REGION | PS_META_DUPLICATE_OK,
    3940                         "Region over which subtraction was performed", subRegion);
     41
     42        psString string = psRegionToString(*subRegion);
    4043        psFree(subRegion);
     44
     45        psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK,
     46                         "Region over which subtraction was performed", string);
     47        psFree(string);
    4148    }
    4249
    … …  
    4552                     PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
    4653    psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE,
    47                      PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->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);
    4857
    4958    // Realisations of kernel
    … …  
    6574                                                           false); // Image of the kernel
    6675                if (!kernel) {
    67                     psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
     76                    psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
    6877                    psFree(convKernels);
    6978                    return false;
    … …  
    7281                if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
    7382                                          (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
    74                     psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
     83                    psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
    7584                    psFree(kernel);
    7685                    psFree(convKernels);
    … …  
    8493                                                      true); // Image of the kernel
    8594                    if (!kernel) {
    86                         psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
     95                        psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
    8796                        psFree(convKernels);
    8897                        return false;
    … …  
    92101                                              (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 4,
    93102                                              (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
    94                         psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
     103                        psError(psErrorCodeLast(), false, "Unable to overlay kernel image.");
    95104                        psFree(kernel);
    96105                        psFree(convKernels);
    … …  
    108117    }
    109118
    110 
    111 #if 0
     119    // sample difference images
     120    {
     121        psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
     122    }
     123
     124#ifdef TESTING
    112125    // Generate images of the kernel components
    113126    {
    114127        psMetadata *header = psMetadataAlloc(); // Header
    115         for (int i = 0; i < solution->n; i++) {
     128        for (int i = 0; i < kernels->solution1->n; i++) {
    116129            psString name = NULL;       // Header keyword
    117130            psStringAppend(&name, "SOLN%04d", i);
    118             psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]);
     131            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution1->data.F64[i]);
    119132            psFree(name);
    120133        }
    121         psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0);
    122         psFits *kernelFile = psFitsOpen("kernels.fits", "w");
     134        psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false);
     135        psFits *kernelFile = psFitsOpen("kernels1.fits", "w");
    123136        (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
    124137        psFitsClose(kernelFile);
    … …  
    126139        psFree(header);
    127140    }
     141    if (kernels->solution2) {
     142        psMetadata *header = psMetadataAlloc(); // Header
     143        for (int i = 0; i < kernels->solution2->n; i++) {
     144            psString name = NULL;       // Header keyword
     145            psStringAppend(&name, "SOLN%04d", i);
     146            psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution2->data.F64[i]);
     147            psFree(name);
     148        }
     149        psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, true);
     150        psFits *kernelFile = psFitsOpen("kernels2.fits", "w");
     151        (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL);
     152        psFitsClose(kernelFile);
     153        psFree(kernelImages);
     154        psFree(header);
     155    }
    128156#endif
    129157
    130 
    131     // Set the variance factors
    132     float vf1 = 1.0, vf2 = 1.0;         // Variance factors for each image
    133     switch (kernels->mode) {
    134       case PM_SUBTRACTION_MODE_1:
    135         vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
    136         break;
    137       case PM_SUBTRACTION_MODE_2:
    138         vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
    139         break;
    140       case PM_SUBTRACTION_MODE_DUAL:
    141         vf1 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, false);
    142         vf2 = pmSubtractionVarianceFactor(kernels, 0.5, 0.5, true);
    143         break;
    144       default:
    145         psAbort("Invalid subtraction mode: %x", kernels->mode);
    146     }
    147 
    148     // Weight by the area
    149     if (region) {
    150         float norm = (region->x1 - region->x0 + 1) * (region->y1 - region->y0 + 1) / (numCols * numRows);
    151         vf1 *= norm;
    152         vf2 *= norm;
    153     }
    154 
    155     // Update the variance factor
    156 #define UPDATE_VARFACTOR(VF, ANALYSIS) { \
    157     psMetadataItem *vfItem = psMetadataLookup(analysis, ANALYSIS); \
    158     if (vfItem) { \
    159         psAssert(vfItem->type == PS_TYPE_F32, "Should be the type we said."); \
    160         vfItem->data.F32 += VF; \
    161     } else { \
    162         psMetadataAddF32(analysis, PS_LIST_TAIL, ANALYSIS, 0, "Variance factor weighted by the area", VF); \
    163     } \
    164 }
    165 
    166     UPDATE_VARFACTOR(vf1, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
    167     UPDATE_VARFACTOR(vf2, PM_SUBTRACTION_ANALYSIS_VARFACTOR_2);
    168158
    169159    // Kernel shape
    … …  
    171161        psImage *image = pmSubtractionKernelImage(kernels, 0.5, 0.5, false); // Image of the kernel
    172162        if (!image) {
    173             psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     163            psError(psErrorCodeLast(), false, "Unable to generate image of kernel.");
    174164            return false;
    175165        }
    … …  
    201191        psFree(image);
    202192
    203         psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous
    204         if (item) {
    205             item->data.F32 = PS_MAX(item->data.F32, max);
    206         } else {
    207             psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0,
    208                              "Maximum deconvolution fraction", max);
     193        {
     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
     203        {
     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            }
    209211        }
    210212    }
    … …  
    214216        psImage *image = pmSubtractionKernelImage(kernels, 0.5, 0.5, false); // Image of the kernel
    215217        if (!image) {
    216             psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
     218            psError(psErrorCodeLast(), false, "Unable to generate image of kernel.");
    217219            return false;
    218220        }
    … …  
    254256        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
    255257                         PS_META_DUPLICATE_OK, "Moment in yy", m02);
     258
     259        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM,
     260                         PS_META_DUPLICATE_OK, "Normalisation", m00);
     261        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX,
     262                         PS_META_DUPLICATE_OK, "Moment in x", m10);
     263        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY,
     264                         PS_META_DUPLICATE_OK, "Moment in y", m01);
     265        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX,
     266                         PS_META_DUPLICATE_OK, "Moment in xx", m20);
     267        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY,
     268                         PS_META_DUPLICATE_OK, "Moment in xy", m11);
     269        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY,
     270                         PS_META_DUPLICATE_OK, "Moment in yy", m02);
    256271    }
    257272
    258273    // Difference in background
    259274    {
    260         psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
    261                                                                   numCols / 2.0, numRows / 2.0); // Polynomial
     275        psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0); // Polynomial
    262276        float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference
    263277
    264278        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
     279                         PS_META_DUPLICATE_OK, "Background difference", bg);
     280        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF,
    265281                         PS_META_DUPLICATE_OK, "Background difference", bg);
    266282        psFree(polyValues);
    … …  
    275291        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation",
    276292                         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_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     302        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     303        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     304        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     305        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     306        psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     307
     308        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     309        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FSIGMA_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     310        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     311        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMIN_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
     312        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_MEAN,  0, "Mean Fractional Sigma of Residuals", kernels->fSigResMean);
     313        psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_FMAX_RES_STDEV, 0, "Mean Fractional Sigma of Residuals", kernels->fSigResStdev);
    277314    }
    278315
Note: See TracChangeset for help on using the changeset viewer.