Changeset 29543 for trunk/psModules/src/imcombine/pmSubtractionMatch.c
- Timestamp:
- Oct 25, 2010, 2:45:41 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionMatch.c
r29004 r29543 28 28 static bool useFFT = true; // Do convolutions using FFT 29 29 30 # define SEPARATE 031 # if (SEPARATE)32 # define SUBMODE PM_SUBTRACTION_EQUATION_NORM33 # else34 30 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL 35 # endif36 31 37 32 //#define TESTING … … 280 275 } 281 276 282 bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) {283 284 if (!readout) return true;285 286 psImage *image = readout->image;287 psImage *mask = readout->mask;288 psImage *variance = readout->variance;289 for (int y = 0; y < image->numRows; y++) {290 for (int x = 0; x < image->numCols; x++) {291 if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue;292 bool valid = false;293 valid = isfinite(image->data.F32[y][x]);294 if (variance) {295 valid &= isfinite(variance->data.F32[y][x]);296 }297 if (valid) continue;298 mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal;299 }300 }301 302 return true;303 }277 // bool pmSubtractionMaskInvalid (const pmReadout *readout, psImageMaskType maskVal) { 278 // 279 // if (!readout) return true; 280 // 281 // psImage *image = readout->image; 282 // psImage *mask = readout->mask; 283 // psImage *variance = readout->variance; 284 // for (int y = 0; y < image->numRows; y++) { 285 // for (int x = 0; x < image->numCols; x++) { 286 // if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) continue; 287 // bool valid = false; 288 // valid = isfinite(image->data.F32[y][x]); 289 // if (variance) { 290 // valid &= isfinite(variance->data.F32[y][x]); 291 // } 292 // if (valid) continue; 293 // mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = maskVal; 294 // } 295 // } 296 // 297 // return true; 298 // } 304 299 305 300 bool pmSubtractionMatchPrecalc(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, … … 392 387 } 393 388 394 pmSubtractionMaskInvalid(ro1, maskVal); 395 pmSubtractionMaskInvalid(ro2, maskVal); 389 // XXX this is done before calling this function 390 // pmSubtractionMaskInvalid(ro1, maskVal); 391 // pmSubtractionMaskInvalid(ro2, maskVal); 396 392 397 393 // General background subtraction, since this is done in pmSubtractionMatch … … 565 561 memCheck("start"); 566 562 567 pmSubtractionMaskInvalid(ro1, maskVal);568 pmSubtractionMaskInvalid(ro2, maskVal);563 // pmSubtractionMaskInvalid(ro1, maskVal); 564 // pmSubtractionMaskInvalid(ro2, maskVal); 569 565 570 566 psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels … … 672 668 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, 673 669 inner, binning, ringsOrder, penalty, bounds, subMode); 674 // pmSubtractionVisualShowKernels(kernels);675 670 } 676 671 … … 678 673 679 674 if (subMode == PM_SUBTRACTION_MODE_UNSURE) { 680 #if 0681 // Get backgrounds682 psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background683 psVector *buffer = NULL;// Buffer for stats684 if (!psImageBackground(bgStats, &buffer, ro1->image, ro1->mask, maskVal, rng)) {685 psError(PM_ERR_DATA, false, "Unable to measure background of image 1.");686 psFree(bgStats);687 psFree(buffer);688 goto MATCH_ERROR;689 }690 float bg1 = psStatsGetValue(bgStats, BG_STAT); // Background for image 1691 if (!psImageBackground(bgStats, &buffer, ro2->image, ro2->mask, maskVal, rng)) {692 psError(PM_ERR_DATA, false, "Unable to measure background of image 2.");693 psFree(bgStats);694 psFree(buffer);695 goto MATCH_ERROR;696 }697 float bg2 = psStatsGetValue(bgStats, BG_STAT); // Background for image 2698 psFree(bgStats);699 psFree(buffer);700 701 pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use702 #endif703 675 pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); 704 676 switch (newMode) { … … 732 704 } 733 705 734 // XXX step 1: calculate normalization 735 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); 706 // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue 707 psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); 708 if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { 709 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); 710 goto MATCH_ERROR; 711 } 712 713 // step 0b : calculate the moments, pass along to the next steps via stamps->normValue 714 // XXX this step is now skipped -- delete 715 psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); 716 if (!pmSubtractionCalculateMoments(kernels, stamps)) { 717 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); 718 goto MATCH_ERROR; 719 } 720 721 // step 1: generate the elements of the matrix equation Ax = B 722 psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); 736 723 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { 737 724 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); … … 739 726 } 740 727 741 psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); 728 // step 2: solve the matrix equation Ax = B 729 psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); 742 730 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { 743 731 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); … … 746 734 memCheck(" solve equation"); 747 735 748 # if (SEPARATE)749 // set USED -> CALCULATE750 pmSubtractionStampsResetStatus (stamps);751 752 // XXX step 2: calculate kernel parameters753 psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");754 if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {755 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");756 goto MATCH_ERROR;757 }758 memCheck(" calculate equation");759 760 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");761 if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {762 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");763 goto MATCH_ERROR;764 }765 memCheck(" solve equation");766 # endif767 736 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations 768 737 if (!deviations) { … … 770 739 goto MATCH_ERROR; 771 740 } 772 773 741 memCheck(" calculate deviations"); 774 742 … … 787 755 // if we hit the max number of iterations and we have rejected stamps, re-solve 788 756 if (numRejected > 0) { 789 // XXX step 1: calculate normalization 757 758 // step 1: generate the elements of the matrix equation Ax = B 790 759 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); 791 760 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { … … 794 763 } 795 764 796 // s olve normalization765 // step 2: solve the matrix equation Ax = B 797 766 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); 798 767 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { … … 801 770 } 802 771 803 # if (SEPARATE)804 // set USED -> CALCULATE805 pmSubtractionStampsResetStatus (stamps);806 807 // XXX step 2: calculate kernel parameters808 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");809 if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {810 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");811 goto MATCH_ERROR;812 }813 814 // solve kernel parameters815 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");816 if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {817 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");818 goto MATCH_ERROR;819 }820 memCheck(" solve equation");821 # endif822 772 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations 823 773 if (!deviations) { … … 837 787 goto MATCH_ERROR; 838 788 } 839 840 789 memCheck("diag outputs"); 841 790 … … 1134 1083 } 1135 1084 1136 # if (SEPARATE)1137 // set USED -> CALCULATE1138 pmSubtractionStampsResetStatus (stamps);1139 1140 psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);1141 if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {1142 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");1143 return false;1144 }1145 1146 psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);1147 if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {1148 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");1149 return false;1150 }1151 # endif1152 1153 1085 psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); 1154 1086 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations … … 1181 1113 } 1182 1114 psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); 1183 1184 # if (SEPARATE)1185 // set USED -> CALCULATE1186 pmSubtractionStampsResetStatus (stamps);1187 1188 psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);1189 if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {1190 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");1191 return false;1192 }1193 1194 psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);1195 if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {1196 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");1197 return false;1198 }1199 psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);1200 # endif1201 1115 1202 1116 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations … … 1288 1202 1289 1203 bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, 1290 float fwhm1, float fwhm2, floatscaleRef, float scaleMin, float scaleMax)1204 float scaleRef, float scaleMin, float scaleMax) 1291 1205 { 1292 1206 PS_ASSERT_PTR_NON_NULL(kernelSize, false); … … 1294 1208 PS_ASSERT_VECTOR_NON_NULL(widths, false); 1295 1209 PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false); 1296 PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);1297 PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);1298 1210 PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); 1299 1211 PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); … … 1301 1213 PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); 1302 1214 1303 // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference 1215 float fwhm1; 1216 float fwhm2; 1217 1218 pmSubtractionGetFWHMs(&fwhm1, &fwhm2); 1219 psAssert(isfinite(fwhm1), "fwhm 1 not set"); 1220 psAssert(isfinite(fwhm2), "fwhm 2 not set"); 1221 1222 // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference 1304 1223 float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor 1305 1306 // XXX save these values in a static for later use1307 pmSubtractionSetFWHMs(fwhm1, fwhm2);1308 1224 1309 1225 if (isfinite(scaleMin) && scale < scaleMin) {
Note:
See TracChangeset
for help on using the changeset viewer.
