Changeset 26893 for trunk/psModules/src/imcombine/pmSubtraction.c
- Timestamp:
- Feb 10, 2010, 7:34:39 PM (16 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/imcombine/pmSubtraction.c (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r26584 r26893 21 21 #include "pmSubtractionStamps.h" 22 22 #include "pmSubtractionEquation.h" 23 #include "pmSubtractionVisual.h" 23 24 #include "pmSubtractionThreads.h" 24 25 … … 63 64 } 64 65 65 // Contribute to an image of the solved kernel component for ISIS66 static void solvedKernel ISIS(psKernel *kernel, // Kernel, updated67 const pmSubtractionKernels *kernels, // Kernel basis functions68 float value, // Normalisation value for basis function69 int index // Index of basis function of interest66 // Contribute to an image of the solved kernel component using the preCalculated image 67 static 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 70 71 ) 71 72 { 72 73 int size = kernels->size; // Kernel half-size 73 p sArray*preCalc = kernels->preCalc->data[index]; // Precalculated values74 pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated values 74 75 #if 0 75 psVector *xKernel = preCalc->data[0]; // Kernel in x76 psVector *yKernel = preCalc->data[1]; // Kernel in y77 76 // Iterating over the kernel 78 77 for (int y = 0, v = -size; v <= size; y++, v++) { 79 float yValue = value * yKernel->data.F32[y];78 float yValue = value * preCalc->yKernel->data.F32[y]; 80 79 for (int x = 0, u = -size; u <= size; x++, u++) { 81 kernel->kernel[v][u] += yValue * xKernel->data.F32[x];80 kernel->kernel[v][u] += yValue * preCalc->xKernel->data.F32[x]; 82 81 } 83 82 } … … 87 86 } 88 87 #else 89 psKernel *k = preCalc->data[2]; // Kernel image90 88 for (int v = -size; v <= size; v++) { 91 89 for (int u = -size; u <= size; u++) { 92 kernel->kernel[v][u] += value * k->kernel[v][u];90 kernel->kernel[v][u] += value * preCalc->kernel->kernel[v][u]; 93 91 } 94 92 } … … 119 117 for (int i = 0; i < numKernels; i++) { 120 118 double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value 119 if (wantDual) { 120 // The model is built with the dual convolution terms added, so to produce zero residual the 121 // equation results in negative coefficients which we must undo. 122 value *= -1.0; 123 } 121 124 122 125 switch (kernels->type) { … … 149 152 case PM_SUBTRACTION_KERNEL_GUNK: { 150 153 if (i < kernels->inner) { 151 solvedKernel ISIS(kernel, kernels, value, i);154 solvedKernelPreCalc(kernel, kernels, value, i); 152 155 } else { 153 156 // Using delta function … … 159 162 break; 160 163 } 161 case PM_SUBTRACTION_KERNEL_ISIS: { 162 solvedKernelISIS(kernel, kernels, value, i); 164 case PM_SUBTRACTION_KERNEL_ISIS: 165 case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: 166 case PM_SUBTRACTION_KERNEL_HERM: 167 case PM_SUBTRACTION_KERNEL_DECONV_HERM: { 168 solvedKernelPreCalc(kernel, kernels, value, i); 163 169 break; 164 170 } 165 171 case PM_SUBTRACTION_KERNEL_RINGS: { 166 psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data 167 psVector *uCoords = preCalc->data[0]; // u coordinates 168 psVector *vCoords = preCalc->data[1]; // v coordinates 169 psVector *poly = preCalc->data[2]; // Polynomial values 170 int num = uCoords->n; // Number of pixels 172 pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; // Precalculated kernels 173 int num = preCalc->uCoords->n; // Number of pixels 171 174 172 175 for (int j = 0; j < num; j++) { 173 int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates 174 kernel->kernel[v][u] += poly->data.F32[j] * value; 176 int u = preCalc->uCoords->data.S32[j]; 177 int v = preCalc->vCoords->data.S32[j]; // Kernel coordinates 178 kernel->kernel[v][u] += preCalc->poly->data.F32[j] * value; 175 179 } 176 180 // Photometric scaling is built into the kernel --- no subtraction! … … 419 423 *target |= maskBad; 420 424 } else if (*source & subConvPoor) { 425 *target &= ~maskBad; 421 426 *target |= maskPoor; 427 } else { 428 *target &= ~maskBad & ~maskPoor; 422 429 } 423 430 } … … 454 461 #endif 455 462 456 // Convolve a stamp using a n ISISkernel basis function457 static psKernel *convolveStamp ISIS(const psKernel *image, // Image to convolve458 const pmSubtractionKernels *kernels, // Kernel basis functions459 int index, // Index of basis function of interest460 int footprint // Half-size of stamp463 // Convolve a stamp using a pre-calculated kernel basis function 464 static psKernel *convolveStampPreCalc(const psKernel *image, // Image to convolve 465 const pmSubtractionKernels *kernels, // Kernel basis functions 466 int index, // Index of basis function of interest 467 int footprint // Half-size of stamp 461 468 ) 462 469 { 463 p sArray*preCalc = kernels->preCalc->data[index]; // Precalculated data470 pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data 464 471 #if 0 465 472 // Convolving using separable convolution 466 psVector *xKernel = preCalc->data[0]; // Kernel in x467 psVector *yKernel = preCalc->data[1]; // Kernel in y468 473 int size = kernels->size; // Size of kernel 469 474 … … 477 482 float value = 0.0; // Value of convolved pixel 478 483 int uMin = x - size, uMax = x + size; // Range for u 479 psF32 *xKernelData = & xKernel->data.F32[xKernel->n - 1]; // Kernel values484 psF32 *xKernelData = &preCalc->xKernel->data.F32[xKernel->n - 1]; // Kernel values 480 485 psF32 *imageData = &image->kernel[y][uMin]; // Image values 481 486 for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) { … … 492 497 float value = 0.0; // Value of convolved pixel 493 498 int vMin = y - size, vMax = y + size; // Range for v 494 psF32 *yKernelData = & yKernel->data.F32[yKernel->n - 1]; // Kernel values499 psF32 *yKernelData = &preCalc->yKernel->data.F32[yKernel->n - 1]; // Kernel values 495 500 psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way! 496 501 for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) { … … 509 514 #else 510 515 // Convolving using precalculated kernel 511 return p_pmSubtractionConvolveStampPrecalc(image, preCalc-> data[2]);516 return p_pmSubtractionConvolveStampPrecalc(image, preCalc->kernel); 512 517 #endif 513 518 } … … 525 530 int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image 526 531 psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version 532 533 // pmSubtractionVisualShowSubtraction(image->image, kernel->image, conv); 534 527 535 psFree(conv); 528 536 return convolved; … … 569 577 } 570 578 579 void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn, 580 int xMin, int xMax, int yMin, int yMax) 581 { 582 float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation 583 *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize; 584 *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize; 585 return; 586 } 587 571 588 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels, 572 int numCols, int numRows, intx, int y)589 int x, int y) 573 590 { 574 591 assert(kernels); 575 assert(numCols > 0 && numRows > 0); 576 577 // Size to use when calculating normalised coordinates (different from actual size when convolving 578 // subimage) 579 int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); 580 int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); 581 582 // Normalised coordinates 583 float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize; 584 float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize; 585 592 593 float xNorm, yNorm; // Normalised coordinates 594 p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y, 595 kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax); 586 596 return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm); 587 597 } … … 664 674 if (index < kernels->inner) { 665 675 // Photometric scaling is already built in to the precalculated kernel 666 return convolveStamp ISIS(image, kernels, index, footprint);676 return convolveStampPreCalc(image, kernels, index, footprint); 667 677 } 668 678 // Using delta function … … 673 683 return convolved; 674 684 } 675 case PM_SUBTRACTION_KERNEL_ISIS: { 676 return convolveStampISIS(image, kernels, index, footprint); 677 } 685 case PM_SUBTRACTION_KERNEL_ISIS: 686 case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: 687 case PM_SUBTRACTION_KERNEL_HERM: 688 case PM_SUBTRACTION_KERNEL_DECONV_HERM: { 689 return convolveStampPreCalc(image, kernels, index, footprint); 690 } 678 691 case PM_SUBTRACTION_KERNEL_RINGS: { 679 psKernel *convolved = psKernelAlloc(-footprint, footprint, 680 -footprint, footprint); // Convolved image 681 psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data 682 psVector *uCoords = preCalc->data[0]; // u coordinates 683 psVector *vCoords = preCalc->data[1]; // v coordinates 684 psVector *poly = preCalc->data[2]; // Polynomial values 685 int num = uCoords->n; // Number of pixels 686 psS32 *uData = uCoords->data.S32, *vData = vCoords->data.S32; // Dereference u,v coordinates 687 psF32 *polyData = poly->data.F32; // Dereference polynomial values 692 psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image 693 pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[index]; // Precalculated data 694 695 int num = preCalc->uCoords->n; // Number of pixels 696 psS32 *uData = preCalc->uCoords->data.S32; // Dereference v coordinate 697 psS32 *vData = preCalc->vCoords->data.S32; // Dereference u coordinate 698 psF32 *polyData = preCalc->poly->data.F32; // Dereference polynomial values 688 699 psF32 **imageData = image->kernel; // Dereference image 689 700 psF32 **convData = convolved->kernel; // Dereference convolved image … … 765 776 766 777 767 768 769 778 int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, 770 779 const psVector *deviations, psImage *subMask, float sigmaRej) … … 860 869 int numGood = 0; // Number of good stamps 861 870 double newMean = 0.0; // New mean 871 psString log = NULL; // Log message 872 psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit); 862 873 for (int i = 0; i < stamps->num; i++) { 863 874 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 872 883 // Mask out the stamp in the image so you it's not found again 873 884 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, 874 (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); 885 (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); 886 psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i, 887 (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), 888 fabsf(deviations->data.F32[i] - mean)); 875 889 numRejected++; 876 890 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { … … 895 909 psFree(stamp->weight); 896 910 stamp->image1 = stamp->image2 = stamp->weight = NULL; 897 psFree(stamp->matrix1); 898 psFree(stamp->matrix2); 899 psFree(stamp->matrixX); 900 stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL; 901 psFree(stamp->vector1); 902 psFree(stamp->vector2); 903 stamp->vector1 = stamp->vector2 = NULL; 911 psFree(stamp->matrix); 912 stamp->matrix = NULL; 913 psFree(stamp->vector); 914 stamp->vector = NULL; 904 915 } else { 905 916 numGood++; … … 910 921 } 911 922 newMean /= numGood; 923 924 if (numRejected == 0) { 925 psStringAppend(&log, "<none>\n"); 926 } 927 psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); 928 psFree(log); 912 929 913 930 if (ds9) { … … 995 1012 psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64); // Backup version 996 1013 997 int num = wantDual ? solution->n - 1 : solution->n;// Number of kernel basis functions1014 int num = kernels->num; // Number of kernel basis functions 998 1015 999 1016 psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // Solved polynomial … … 1060 1077 // Only generate polynomial values every kernel footprint, since we have already assumed 1061 1078 // (with the stamps) that it does not vary rapidly on this scale. 1062 psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows, 1063 xMin + x0 + size + 1, 1064 yMin + y0 + size + 1); 1079 psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1, 1080 yMin + y0 + size + 1); // Polynomial 1065 1081 float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term 1066 1082 … … 1137 1153 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, 1138 1154 psImage *subMask, int stride, psImageMaskType maskBad, psImageMaskType maskPoor, 1139 float poorFrac, float kernelError, const psRegion *region,1155 float poorFrac, float kernelError, float covarFrac, const psRegion *region, 1140 1156 const pmSubtractionKernels *kernels, bool doBG, bool useFFT) 1141 1157 { … … 1146 1162 PM_ASSERT_READOUT_NON_NULL(ro1, false); 1147 1163 PM_ASSERT_READOUT_IMAGE(ro1, false); 1164 PM_ASSERT_READOUT_IMAGE(out1, false); 1148 1165 numCols = ro1->image->numCols; 1149 1166 numRows = ro1->image->numRows; … … 1155 1172 PM_ASSERT_READOUT_NON_NULL(ro2, false); 1156 1173 PM_ASSERT_READOUT_IMAGE(ro2, false); 1174 PM_ASSERT_READOUT_IMAGE(out2, false); 1157 1175 if (numCols == 0 && numRows == 0) { 1158 1176 numCols = ro2->image->numCols; … … 1177 1195 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(kernelError, 0.0, false); 1178 1196 PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(kernelError, 1.0, false); 1197 PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(covarFrac, 0.0, false); 1198 PS_ASSERT_FLOAT_LESS_THAN(covarFrac, 1.0, false); 1179 1199 if (region && psRegionIsNaN(*region)) { 1180 1200 psString string = psRegionToString(*region); … … 1188 1208 bool threaded = pmSubtractionThreaded(); // Running threaded? 1189 1209 1190 // Outputs1191 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {1192 if (!out1->image) {1193 out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);1194 }1195 if (ro1->variance) {1196 if (!out1->variance) {1197 out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);1198 }1199 psImageInit(out1->variance, 0.0);1200 }1201 }1202 if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {1203 if (!out2->image) {1204 out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);1205 }1206 if (ro2->variance) {1207 if (!out2->variance) {1208 out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);1209 }1210 psImageInit(out2->variance, 0.0);1211 }1212 }1213 1210 psImage *convMask = NULL; // Convolved mask image (common to inputs 1 and 2) 1214 1211 if (subMask) { 1215 1212 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1216 if (!out1->mask) {1217 out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);1218 }1219 psImageInit(out1->mask, 0);1220 1213 convMask = out1->mask; 1221 1214 } 1222 1215 if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1223 if (!out2->mask) {1224 out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);1225 }1226 psImageInit(out2->mask, 0);1227 1216 if (!convMask) { 1228 1217 convMask = out2->mask; … … 1244 1233 1245 1234 // Get region for convolution: [xMin:xMax,yMin:yMax] 1246 int xMin = size, xMax = numCols- size;1247 int yMin = size, yMax = numRows- size;1235 int xMin = kernels->xMin + size, xMax = kernels->xMax - size; 1236 int yMin = kernels->yMin + size, yMax = kernels->yMax - size; 1248 1237 if (region) { 1249 1238 xMin = PS_MAX(region->x0, xMin); … … 1336 1325 1337 1326 // Calculate covariances 1338 // This can take a while, so we only do it for a single instance 1339 // XXX psImageCovarianceCalculate could be multithreaded 1327 // This can be fairly involved, so we only do it for a single instance 1328 // Enable threads for covariance calculation, since we're not threading on top of it. 1329 oldThreads = psImageCovarianceSetThreads(true); 1340 1330 if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 1341 1331 psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel 1332 psKernelTruncate(kernel, covarFrac); 1342 1333 out1->covariance = psImageCovarianceCalculate(kernel, ro1->covariance); 1343 1334 psFree(kernel); … … 1346 1337 psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, 1347 1338 kernels->mode == PM_SUBTRACTION_MODE_DUAL); // Conv. kernel 1339 psKernelTruncate(kernel, covarFrac); 1348 1340 out2->covariance = psImageCovarianceCalculate(kernel, ro2->covariance); 1349 1341 psFree(kernel); 1350 1342 } 1343 psImageCovarianceSetThreads(oldThreads); 1351 1344 1352 1345 // Copy anything that wasn't convolved … … 1388 1381 } 1389 1382 1390 1391 1383 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve image: %f sec", 1392 1384 psTimerClear("pmSubtractionConvolve"));
Note:
See TracChangeset
for help on using the changeset viewer.
