IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26667


Ignore:
Timestamp:
Jan 22, 2010, 10:58:17 AM (16 years ago)
Author:
Paul Price
Message:

Move scaling into psModules so ppStack can use it too.

Location:
branches/eam_branches/20091201
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/ppSub/src/ppSubMatchPSFs.c

    r26649 r26667  
    131131    }
    132132
    133 
    134133    float inFWHM = subImagePSF(data, inRO, inSources); // FWHM for input
    135134    float refFWHM = subImagePSF(data, refRO, refSources); // FWHM for reference
    … …  
    139138        return true;
    140139    }
    141 
    142 //    float diff = sqrtf(PS_SQR(PS_MAX(inFWHM, refFWHM)) - PS_SQR(PS_MIN(inFWHM, refFWHM))); // Difference
    143     float diff = PS_MAX(inFWHM, refFWHM);
    144140
    145141    float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling
    … …  
    153149    }
    154150
    155     float scale = diff / scaleRef;      // Scaling factor
    156     if (scale < scaleMin) {
    157         scale = scaleMin;
    158     }
    159     if (scale > scaleMax) {
    160         scale = scaleMax;
    161     }
    162 
    163     for (int i = 0; i < kernelWidths->n; i++) {
    164         kernelWidths->data.F32[i] *= scale;
    165     }
    166     *kernelSize = *kernelSize * scale + 0.5;
    167     *stampSize = *stampSize * scale + 0.5;
    168 
    169     psLogMsg("ppSub", PS_LOG_INFO, "Scaled kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     151    if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths, inFWHM, refFWHM,
     152                                  scaleRef, scaleMin, scaleMax)) {
     153        psError(PS_ERR_UNKNOWN, false, "Unable to scale parameters.");
     154        return false;
     155    }
    170156
    171157    return true;

  • branches/eam_branches/20091201/ppSub/src/ppSubReadoutJpeg.c

    r26597 r26667  
    5050}
    5151
    52 bool ppSubResidualSampleJpeg(pmConfig *config) {
     52bool ppSubResidualSampleJpeg(pmConfig *config)
     53{
     54    return true;
    5355
    5456    pmFPAview *view = ppSubViewReadout(); // View to readout
    … …  
    6769    psMetadataItem *item = NULL;// Item from iteration
    6870    while ((item = psMetadataGetAndIncrement(iter))) {
    69         assert(item->type == PS_DATA_ARRAY);
    70         psArray *sampleStamps = item->data.V;
    71         samples = psArrayAdd(samples, 16, sampleStamps);
     71        assert(item->type == PS_DATA_ARRAY);
     72        psArray *sampleStamps = item->data.V;
     73        samples = psArrayAdd(samples, 16, sampleStamps);
    7274    }
    7375    psFree(iter);
    … …  
    108110    for (int i = 0; i < samples->n; i++) {
    109111
    110         int xBlock = i % NXblock;
    111         int yBlock = i / NYblock;
    112              
    113         psArray *kernels = samples->data[i];
     112        int xBlock = i % NXblock;
     113        int yBlock = i / NYblock;
    114114
    115         for (int j = 0; j < kernels->n; j++) {
     115        psArray *kernels = samples->data[i];
    116116
    117             psImage *kernel = kernels->data[j];
     117        for (int j = 0; j < kernels->n; j++) {
    118118
    119             int xSub = j % 3;
    120             int ySub = j / 3;
     119            psImage *kernel = kernels->data[j];
    121120
    122             int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder);
    123             int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder);
     121            int xSub = j % 3;
     122            int ySub = j / 3;
    124123
    125             for (int y = 0; y < kernel->numRows; y++) {
    126                 for (int x = 0; x < kernel->numCols; x++) {
    127                     readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x];
    128                 }
    129             }
    130         }
     124            int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder);
     125            int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder);
     126
     127            for (int y = 0; y < kernel->numRows; y++) {
     128                for (int x = 0; x < kernel->numCols; x++) {
     129                    readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x];
     130                }
     131            }
     132        }
    131133    }
    132134
    133135    {
    134         psFits *fits = psFitsOpen ("resid.stamps.fits", "w");
    135         psFitsWriteImage (fits, NULL, readout->image, 0, NULL);
    136         psFitsClose (fits);
     136        psFits *fits = psFitsOpen ("resid.stamps.fits", "w");
     137        psFitsWriteImage (fits, NULL, readout->image, 0, NULL);
     138        psFitsClose (fits);
    137139    }
    138140

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c

    r26579 r26667  
    6666// Contribute to an image of the solved kernel component using the preCalculated image
    6767static void solvedKernelPreCalc(psKernel *kernel, // Kernel, updated
    68                                 const pmSubtractionKernels *kernels, // Kernel basis functions
    69                                 float value,                         // Normalisation value for basis function
    70                                 int index                  // Index of basis function of interest
     68                                const pmSubtractionKernels *kernels, // Kernel basis functions
     69                                float value,                         // Normalisation value for basis function
     70                                int index                  // Index of basis function of interest
    7171    )
    7272{
    … …  
    162162              break;
    163163          }
    164           case PM_SUBTRACTION_KERNEL_ISIS: 
    165           case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: 
    166           case PM_SUBTRACTION_KERNEL_HERM:               
     164          case PM_SUBTRACTION_KERNEL_ISIS:
     165          case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
     166          case PM_SUBTRACTION_KERNEL_HERM:
    167167          case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
    168               solvedKernelPreCalc(kernel, kernels, value, i);             
    169               break;
     168              solvedKernelPreCalc(kernel, kernels, value, i);
     169              break;
    170170          }
    171171          case PM_SUBTRACTION_KERNEL_RINGS: {
    … …  
    460460// Convolve a stamp using a pre-calculated kernel basis function
    461461static psKernel *convolveStampPreCalc(const psKernel *image, // Image to convolve
    462                                       const pmSubtractionKernels *kernels, // Kernel basis functions
    463                                       int index,                            // Index of basis function of interest
    464                                       int footprint                         // Half-size of stamp
     462                                      const pmSubtractionKernels *kernels, // Kernel basis functions
     463                                      int index,                            // Index of basis function of interest
     464                                      int footprint                         // Half-size of stamp
    465465    )
    466466{
    … …  
    678678          return convolved;
    679679      }
    680       case PM_SUBTRACTION_KERNEL_ISIS: 
    681       case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: 
     680      case PM_SUBTRACTION_KERNEL_ISIS:
     681      case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
    682682      case PM_SUBTRACTION_KERNEL_HERM:
    683683      case PM_SUBTRACTION_KERNEL_DECONV_HERM: {
    684             return convolveStampPreCalc(image, kernels, index, footprint);
    685         }
     684            return convolveStampPreCalc(image, kernels, index, footprint);
     685        }
    686686      case PM_SUBTRACTION_KERNEL_RINGS: {
    687687          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
    … …  
    864864    int numGood = 0;                    // Number of good stamps
    865865    double newMean = 0.0;               // New mean
     866    psString log = NULL;                // Log message
     867    psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit);
    866868    for (int i = 0; i < stamps->num; i++) {
    867869        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    … …  
    877879                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
    878880                        (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
     881                psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i,
     882                               (int)(stamp->x + 0.5), (int)(stamp->y + 0.5),
     883                               fabsf(deviations->data.F32[i] - mean));
    879884                numRejected++;
    880885                for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
    … …  
    911916    }
    912917    newMean /= numGood;
     918
     919    if (numRejected == 0) {
     920        psStringAppend(&log, "<none>\n");
     921    }
     922    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
     923    psFree(log);
    913924
    914925    if (ds9) {

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

    r26609 r26667  
    14251425    }
    14261426
     1427    psString log = psStringCopy("Deviations:\n");               // Log message with deviations
    14271428    for (int i = 0; i < stamps->num; i++) {
    14281429        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
    … …  
    15561557        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
    15571558                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
     1559        psStringAppend(&log, "Stamp %d (%d,%d): %f\n",
     1560                       i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
    15581561        if (!isfinite(deviations->data.F32[i])) {
    15591562            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    … …  
    16001603
    16011604    }
     1605
     1606    psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log);
     1607    psFree(log);
    16021608
    16031609    // calculate and report the normalization and background for the image center

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c

    r26580 r26667  
    496496
    497497
    498             // generate the window function from the set of stamps
    499             if (!pmSubtractionStampsGetWindow(stamps, size)) {
    500                 goto MATCH_ERROR;
    501             }
     498            // generate the window function from the set of stamps
     499            if (!pmSubtractionStampsGetWindow(stamps, size)) {
     500                goto MATCH_ERROR;
     501            }
    502502
    503503            // Define kernel basis functions
    … …  
    514514                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
    515515                                                       inner, binning, ringsOrder, penalty, subMode);
    516                 // pmSubtractionVisualShowKernels(kernels);
     516                // pmSubtractionVisualShowKernels(kernels);
    517517            }
    518518
    … …  
    568568                }
    569569
    570                 // generate the window function from the set of stamps
    571                 if (!pmSubtractionStampsGetWindow(stamps, size)) {
    572                     goto MATCH_ERROR;
    573                 }
    574 
    575                 // XXX step 1: calculate normalization
     570                // generate the window function from the set of stamps
     571                if (!pmSubtractionStampsGetWindow(stamps, size)) {
     572                    goto MATCH_ERROR;
     573                }
     574
     575                // XXX step 1: calculate normalization
    576576                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    577577                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    … …  
    588588
    589589# if (SEPARATE)
    590                 // set USED -> CALCULATE
    591                 pmSubtractionStampsResetStatus (stamps);
    592 
    593                 // XXX step 2: calculate kernel parameters
     590                // set USED -> CALCULATE
     591                pmSubtractionStampsResetStatus (stamps);
     592
     593                // XXX step 2: calculate kernel parameters
    594594                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
    595595                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
    … …  
    626626            }
    627627
    628             // if we hit the max number of iterations and we have rejected stamps, re-solve
     628            // if we hit the max number of iterations and we have rejected stamps, re-solve
    629629            if (numRejected > 0) {
    630                 // XXX step 1: calculate normalization
     630                // XXX step 1: calculate normalization
    631631                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    632632                if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    … …  
    635635                }
    636636
    637                 // solve normalization
     637                // solve normalization
    638638                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
    639639                if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    … …  
    643643
    644644# if (SEPARATE)
    645                 // set USED -> CALCULATE
    646                 pmSubtractionStampsResetStatus (stamps);
    647 
    648                 // XXX step 2: calculate kernel parameters
     645                // set USED -> CALCULATE
     646                pmSubtractionStampsResetStatus (stamps);
     647
     648                // XXX step 2: calculate kernel parameters
    649649                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    650650                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
    … …  
    653653                }
    654654
    655                 // solve kernel parameters
     655                // solve kernel parameters
    656656                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
    657657                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
    … …  
    992992        return false;
    993993    }
    994 # endif 
     994# endif
    995995
    996996    psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
    … …  
    10121012    if (numRejected > 0) {
    10131013        // Allow re-fit with reduced stamps set
    1014         psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1015         if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
    1016             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    1017             return false;
    1018         }
     1014        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1015        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
     1016            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1017            return false;
     1018        }
    10191019
    10201020        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    … …  
    10261026
    10271027# if (SEPARATE)
    1028         // set USED -> CALCULATE
    1029         pmSubtractionStampsResetStatus (stamps);
    1030        
    1031         psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1032         if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
    1033             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    1034             return false;
    1035         }
     1028        // set USED -> CALCULATE
     1029        pmSubtractionStampsResetStatus (stamps);
     1030
     1031        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
     1032        if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
     1033            psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     1034            return false;
     1035        }
    10361036
    10371037        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    … …  
    11291129}
    11301130
     1131
     1132bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
     1133                              float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
     1134{
     1135    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
     1136    PS_ASSERT_PTR_NON_NULL(stampSize, false);
     1137    PS_ASSERT_VECTOR_NON_NULL(widths, false);
     1138    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
     1139    PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
     1140    PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
     1141    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
     1142    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
     1143    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false);
     1144    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
     1145
     1146//    float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
     1147    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
     1148
     1149    if (isfinite(scaleMin) && scale < scaleMin) {
     1150        scale = scaleMin;
     1151    }
     1152    if (isfinite(scaleMax) && scale > scaleMax) {
     1153        scale = scaleMax;
     1154    }
     1155
     1156    for (int i = 0; i < widths->n; i++) {
     1157        widths->data.F32[i] *= scale;
     1158    }
     1159    *kernelSize = *kernelSize * scale + 0.5;
     1160    *stampSize = *stampSize * scale + 0.5;
     1161
     1162    psLogMsg("psModules.imcombine", PS_LOG_INFO,
     1163             "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize);
     1164
     1165    return true;
     1166}

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.h

    r26505 r26667  
    9595    );
    9696
     97
     98/// Scale subtraction parameters according to the FWHMs of the inputs
     99bool pmSubtractionParamsScale(
     100    int *kernelSize,                    ///< Half-size of the kernel
     101    int *stampSize,                     ///< Half-size of the stamp (footprint)
     102    psVector *widths,                   ///< ISIS widths
     103    float fwhm1, float fwhm2,           ///< FWHMs for inputs
     104    float scaleRef,                     ///< Reference width for scaling
     105    float scaleMin,                     ///< Minimum scaling ratio, or NAN
     106    float scaleMax                      ///< Maximum scaling ratio, or NAN
     107    );
     108
     109
    97110#endif
Note: See TracChangeset for help on using the changeset viewer.