Changeset 18287
- Timestamp:
- Jun 23, 2008, 12:41:08 PM (18 years ago)
- Location:
- trunk/psModules/src/imcombine
- Files:
-
- 9 edited
-
pmSubtraction.c (modified) (5 diffs)
-
pmSubtractionEquation.c (modified) (12 diffs)
-
pmSubtractionEquation.h (modified) (1 diff)
-
pmSubtractionKernels.c (modified) (31 diffs)
-
pmSubtractionKernels.h (modified) (10 diffs)
-
pmSubtractionMatch.c (modified) (5 diffs)
-
pmSubtractionMatch.h (modified) (1 diff)
-
pmSubtractionParams.c (modified) (2 diffs)
-
pmSubtractionParams.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtraction.c
r18190 r18287 3 3 * @author Paul Price, IfA 4 4 * @author GLG, MHPCC 5 *6 * @version $Revision: 1.95 $ $Name: not supported by cvs2svn $7 * @date $Date: 2008-06-19 03:03:53 $8 5 * 9 6 * Copyright 2004-2007 Institute for Astronomy, University of Hawaii … … 512 509 PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); 513 510 514 double totalSquareDev = 0.0; // Total square deviation from zero 511 // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but 512 // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like 513 // with large N. Therefore, let's just treat this as a Gaussian distribution. 514 515 double mean = 0.0; // Mean deviation 515 516 int numStamps = 0; // Number of used stamps 516 517 for (int i = 0; i < stamps->num; i++) { … … 519 520 continue; 520 521 } 521 totalSquareDev += PS_SQR(deviations->data.F32[i]);522 mean += deviations->data.F32[i]; 522 523 numStamps++; 523 524 } 524 525 float rms = sqrt(totalSquareDev / (double)numStamps); // Convert to RMS 525 mean /= numStamps; 526 527 double rms = 0.0; // Standard deviation 528 for (int i = 0; i < stamps->num; i++) { 529 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 530 if (stamp->status != PM_SUBTRACTION_STAMP_USED) { 531 continue; 532 } 533 rms += PS_SQR(deviations->data.F32[i] - mean); 534 } 535 rms = sqrt(rms / (numStamps - 1)); 526 536 527 537 if (rmsPtr) { … … 549 559 int numRejected = 0; // Number of stamps rejected 550 560 int numGood = 0; // Number of good stamps 551 double new SquareDev = 0.0; // New square deviation561 double newMean = 0.0; // New mean 552 562 for (int i = 0; i < stamps->num; i++) { 553 563 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 554 564 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 555 if (deviations->data.F32[i] > limit) { 565 // Should we reject stars with low deviation? Well, if this is really a Gaussian-like 566 // distribution and they're low, then we have the right to ask why. Isn't it suspicious that 567 // they're anomalously low, compared to the rest of the population which (we hope) is indicative 568 // of normality? Besides, the standard deviation is going to be blown up by stars that didn't 569 // subtract well, in which case very few (if any) stars will be legitimately rejected for being 570 // low. 571 if (deviations->data.F32[i] - mean > limit) { 556 572 // Mask out the stamp in the image so you it's not found again 557 573 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, … … 587 603 } else { 588 604 numGood++; 589 new SquareDev += PS_SQR(deviations->data.F32[i]);605 newMean += deviations->data.F32[i]; 590 606 } 591 607 } 592 608 } 609 newMean /= numGood; 593 610 594 611 if (numRejected > 0) { 595 612 psLogMsg("psModules.imcombine", PS_LOG_INFO, 596 "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", 597 numGood, numRejected, rms, 598 sqrt(newSquareDev / (double)numGood)); 613 "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n", 614 numGood, numRejected, mean, newMean); 599 615 } else { 600 616 psLogMsg("psModules.imcombine", PS_LOG_INFO, 601 "%d good stamps; 0 rejected.\n RMS deviation: %f\n",602 numGood, rms);617 "%d good stamps; 0 rejected.\nMean deviation: %lf\n", 618 numGood, mean); 603 619 } 604 620 -
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r17825 r18287 30 30 for (int y = - footprint; y <= footprint; y++) { 31 31 for (int x = - footprint; x <= footprint; x++) { 32 sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x];32 sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x]; 33 33 } 34 34 } … … 194 194 for (int y = - footprint; y <= footprint; y++) { 195 195 for (int x = - footprint; x <= footprint; x++) { 196 sumC += conv->kernel[y][x] / weight->kernel[y][x];196 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; 197 197 } 198 198 } … … 217 217 for (int y = - footprint; y <= footprint; y++) { 218 218 for (int x = - footprint; x <= footprint; x++) { 219 double invNoise2 = 1.0 / weight->kernel[y][x];219 double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x]; 220 220 double value = input->kernel[y][x] * invNoise2; 221 221 sumI += value; … … 276 276 for (int y = - footprint; y <= footprint; y++) { 277 277 for (int x = - footprint; x <= footprint; x++) { 278 sumTC += target->kernel[y][x] * conv->kernel[y][x]/ weight->kernel[y][x];278 sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; 279 279 } 280 280 } … … 296 296 for (int y = - footprint; y <= footprint; y++) { 297 297 for (int x = - footprint; x <= footprint; x++) { 298 float value = target->kernel[y][x] / weight->kernel[y][x];298 float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x]; 299 299 sumIT += value * input->kernel[y][x]; 300 300 sumT += value; … … 365 365 for (int y = - footprint; y <= footprint; y++) { 366 366 for (int x = - footprint; x <= footprint; x++) { 367 sumC += conv->kernel[y][x] / weight->kernel[y][x];367 sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x]; 368 368 } 369 369 } … … 382 382 } 383 383 384 385 // Add in penalty term to least-squares vector 386 static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term 387 const pmSubtractionKernels *kernels // Kernel parameters 388 ) 389 { 390 if (kernels->penalty == 0.0) { 391 return true; 392 } 393 394 psVector *penalties = kernels->penalties; // Penalties for each kernel component 395 int spatialOrder = kernels->spatialOrder; // Order of spatial variations 396 int numKernels = kernels->num; // Number of kernel components 397 for (int i = 0; i < numKernels; i++) { 398 for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { 399 for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { 400 vector->data.F64[index] -= penalties->data.F32[i]; 401 } 402 } 403 } 404 405 return true; 406 } 407 384 408 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 385 409 // Semi-public functions 386 // XXX EAM these cannot be inline :: talk to Josh 410 // XXX We might like to define these functions as "extern inline" but gcc currently doesn't handle this in c99 411 // mode. See http://gcc.gnu.org/ml/gcc/2006-11/msg00006.html 387 412 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 388 413 … … 675 700 psVectorInit(sumVector, 0.0); 676 701 psImageInit(sumMatrix, 0.0); 702 int numStamps = 0; // Number of good stamps 677 703 for (int i = 0; i < stamps->num; i++) { 678 704 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 680 706 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); 681 707 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); 682 } 683 } 708 numStamps++; 709 } 710 } 711 calculatePenalty(sumVector, kernels); 684 712 685 713 psVector *permutation = NULL; // Permutation vector, required for LU decomposition … … 716 744 psVectorInit(sumVector2, 0.0); 717 745 746 int numStamps = 0; // Number of good stamps 718 747 for (int i = 0; i < stamps->num; i++) { 719 748 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 724 753 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); 725 754 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); 726 } 727 } 755 numStamps++; 756 } 757 } 758 calculatePenalty(sumVector1, kernels); 759 calculatePenalty(sumVector2, kernels); 728 760 729 761 #if 0 … … 778 810 psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A); 779 811 assert(F->numRows == numParams && F->numCols == numParams); 780 float det = NAN;812 float det = 0.0; 781 813 psImage *Fi = psMatrixInvert(NULL, F, &det); 782 814 assert(Fi->numRows == numParams && Fi->numCols == numParams); -
trunk/psModules/src/imcombine/pmSubtractionEquation.h
r15809 r18287 21 21 22 22 /// Calculate the value of a polynomial, specified by coefficients and polynomial values 23 inlinedouble p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients24 25 26 27 23 double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients 24 const psImage *polyValues, ///< Polynomial values 25 int order, ///< Order of polynomials 26 int index, ///< Index at which to begin 27 int step ///< Step between subsequent indices 28 28 ); 29 29 30 30 /// Return the specified coefficient in the solution 31 inlinedouble p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters32 33 34 31 double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters 32 const psImage *polyValues, ///< Polynomial values 33 int index, ///< Coefficient index to calculate 34 bool wantDual ///< Calculate the coefficient for the dual solution? 35 35 ); 36 36 37 37 /// Return the normalisation in the solution 38 inlinedouble p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters38 double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters 39 39 ); 40 40 41 41 /// Return the background (difference) in the solution 42 inlinedouble p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters43 42 double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters 43 const psImage *polyValues ///< Polynomial values 44 44 ); 45 45 -
trunk/psModules/src/imcombine/pmSubtractionKernels.c
r18146 r18287 24 24 psFree(kernels->vStop); 25 25 psFree(kernels->preCalc); 26 psFree(kernels->penalties); 26 27 psFree(kernels->solution1); 27 28 psFree(kernels->solution2); … … 60 61 kernels->v = psVectorRealloc(kernels->v, start + numNew); 61 62 kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); 63 kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew); 62 64 kernels->inner = start; 63 65 … … 74 76 kernels->v->data.S32[index] = v; 75 77 kernels->preCalc->data[index] = NULL; 78 kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v)); 76 79 77 80 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); … … 84 87 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, 85 88 const psVector *fwhms, const psVector *orders, 86 pmSubtractionMode mode)89 float penalty, pmSubtractionMode mode) 87 90 { 88 91 PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); … … 104 107 } 105 108 106 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, 107 s ize, spatialOrder, mode); // The kernels108 psStringAppend(&kernels->description, "ISIS(%d,%s,%d )", size, params, spatialOrder);109 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, 110 spatialOrder, penalty, mode); // The kernels 111 psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); 109 112 110 113 psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements", … … 122 125 psKernel *preCalc = psKernelAlloc(-size, size, -size, size); 123 126 double sum = 0.0; // Normalisation 127 double moment = 0.0; // Moment, for penalty 124 128 for (int v = -size; v <= size; v++) { 125 129 for (int u = -size; u <= size; u++) { 126 130 sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * 127 131 expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); 132 moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v)); 128 133 } 129 134 } … … 146 151 } 147 152 kernels->preCalc->data[index] = preCalc; 148 149 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, 150 fwhms->data.F32[i], uOrder, vOrder); 153 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 154 155 psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, 156 fwhms->data.F32[i], uOrder, vOrder, fabsf(moment)); 151 157 } 152 158 } … … 161 167 162 168 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, 163 int size, int spatialOrder, pmSubtractionMode mode) 169 int size, int spatialOrder, float penalty, 170 pmSubtractionMode mode) 164 171 { 165 172 pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return … … 173 180 kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); 174 181 kernels->preCalc = psArrayAlloc(numBasisFunctions); 182 kernels->penalty = penalty; 183 kernels->penalties = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); 175 184 kernels->uStop = NULL; 176 185 kernels->vStop = NULL; … … 188 197 } 189 198 190 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) 199 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, 200 pmSubtractionMode mode) 191 201 { 192 202 PS_ASSERT_INT_POSITIVE(size, NULL); … … 195 205 int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions 196 206 197 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, 198 s ize, spatialOrder, mode); // The kernels199 psStringAppend(&kernels->description, "POIS(%d,%d )", size, spatialOrder);207 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, size, 208 spatialOrder, penalty, mode); // The kernels 209 psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); 200 210 psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", 201 211 size, spatialOrder, num); … … 211 221 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, 212 222 const psVector *fwhms, const psVector *orders, 213 pmSubtractionMode mode)214 { 215 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 216 fwhms, orders, mode); // Kernels223 float penalty, pmSubtractionMode mode) 224 { 225 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, 226 penalty, mode); // Kernels 217 227 if (!kernels) { 218 228 return NULL; … … 242 252 243 253 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, 244 pmSubtractionMode mode)254 float penalty, pmSubtractionMode mode) 245 255 { 246 256 PS_ASSERT_INT_POSITIVE(size, NULL); … … 262 272 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 263 273 264 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, 265 size, spatialOrder, mode); // The kernels 266 psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); 274 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, 275 spatialOrder, penalty, mode); // The kernels 276 psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, 277 penalty); 267 278 268 279 psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", … … 324 335 psFree(widths); 325 336 337 psWarning("Kernel penalty for dual-convolution is not configured for SPAM kernels."); 338 psVectorInit(kernels->penalties, 0.0); 339 326 340 return kernels; 327 341 } 328 342 329 343 330 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) 344 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty, 345 pmSubtractionMode mode) 331 346 { 332 347 PS_ASSERT_INT_POSITIVE(size, NULL); … … 354 369 psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); 355 370 356 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, 357 s ize, spatialOrder, mode); // The kernels358 psStringAppend(&kernels->description, "FRIES(%d,%d,%d )", size, inner, spatialOrder);371 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, 372 spatialOrder, penalty, mode); // The kernels 373 psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); 359 374 360 375 psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", … … 414 429 psFree(stop); 415 430 431 psWarning("Kernel penalty for dual-convolution is not configured for FRIES kernels."); 432 psVectorInit(kernels->penalties, 0.0); 433 416 434 return kernels; 417 435 } … … 419 437 // Grid United with Normal Kernel 420 438 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, 421 const psVector *orders, int inner, pmSubtractionMode mode) 439 const psVector *orders, int inner, float penalty, 440 pmSubtractionMode mode) 422 441 { 423 442 PS_ASSERT_INT_POSITIVE(size, NULL); … … 431 450 PS_ASSERT_INT_LESS_THAN(inner, size, NULL); 432 451 433 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 434 fwhms, orders, mode); // Kernels 452 // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed 453 454 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, 455 penalty, mode); // Kernels 435 456 psStringPrepend(&kernels->description, "GUNK="); 436 457 psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); … … 447 468 // RINGS --- just what it says 448 469 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, 449 pmSubtractionMode mode)470 float penalty, pmSubtractionMode mode) 450 471 { 451 472 PS_ASSERT_INT_POSITIVE(size, NULL); … … 477 498 int num = numRings * numPoly; // Total number of basis functions 478 499 479 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, 480 size, spatialOrder, mode); // The kernels 481 psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); 500 pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, 501 spatialOrder, penalty, mode); // The kernels 502 psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, 503 penalty); 482 504 483 505 psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", … … 520 542 psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords 521 543 psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial 544 double moment = 0.0; // Moment, for penalty 522 545 523 546 if (i == 0) { … … 527 550 uCoords->n = vCoords->n = poly->n = 1; 528 551 radiusLast = 0; 552 moment = 0.0; 529 553 } else { 530 554 int j = 0; // Index for data … … 546 570 poly->data.F32[j] = polyVal; 547 571 norm += polyVal; 572 moment += polyVal * (PS_SQR(u) + PS_SQR(v)); 548 573 549 574 psVectorExtend(uCoords, RINGS_BUFFER, 1); … … 571 596 psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); 572 597 } 598 // moment /= norm; 573 599 } 574 600 … … 578 604 kernels->u->data.S32[index] = uOrder; 579 605 kernels->v->data.S32[index] = vOrder; 606 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment); 580 607 581 608 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, … … 590 617 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, 591 618 const psVector *fwhms, const psVector *orders, int inner, 592 int binning, int ringsOrder, pmSubtractionMode mode) 619 int binning, int ringsOrder, float penalty, 620 pmSubtractionMode mode) 593 621 { 594 622 switch (type) { 595 623 case PM_SUBTRACTION_KERNEL_POIS: 596 return pmSubtractionKernelsPOIS(size, spatialOrder, mode);624 return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode); 597 625 case PM_SUBTRACTION_KERNEL_ISIS: 598 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode);626 return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode); 599 627 case PM_SUBTRACTION_KERNEL_SPAM: 600 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode);628 return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode); 601 629 case PM_SUBTRACTION_KERNEL_FRIES: 602 return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode);630 return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode); 603 631 case PM_SUBTRACTION_KERNEL_GUNK: 604 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode);632 return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode); 605 633 case PM_SUBTRACTION_KERNEL_RINGS: 606 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode);634 return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode); 607 635 default: 608 636 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); … … 657 685 int binning = 0; // Binning to use 658 686 int ringsOrder = 0; // Polynomial order for rings 687 float penalty = 0.0; // Penalty for wideness 659 688 660 689 if (strncmp(description, "ISIS", 4) == 0) { … … 684 713 685 714 ptr++; // Eat ',' 686 spatialOrder = parseStringInt(ptr); 715 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 716 penalty = parseStringFloat(ptr); 687 717 } 688 718 } else if (strncmp(description, "RINGS", 5) == 0) { … … 692 722 PARSE_STRING_NUMBER(inner, ptr, ',', parseStringInt); 693 723 PARSE_STRING_NUMBER(ringsOrder, ptr, ',', parseStringInt); 694 PARSE_STRING_NUMBER(spatialOrder, ptr, ')', parseStringInt); 724 PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt); 725 PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt); 695 726 } else { 696 727 psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported."); … … 699 730 700 731 return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, 701 inner, binning, ringsOrder, mode);732 inner, binning, ringsOrder, penalty, mode); 702 733 } 703 734 -
trunk/psModules/src/imcombine/pmSubtractionKernels.h
r18146 r18287 33 33 psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) 34 34 psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) 35 float penalty; ///< Penalty for wideness 36 psVector *penalties; ///< Penalty for each kernel component 35 37 int size; ///< The half-size of the kernel 36 38 int inner; ///< The size of an inner region … … 107 109 int size, ///< Half-size of kernel 108 110 int spatialOrder, ///< Order of spatial variations 111 float penalty, ///< Penalty for wideness 109 112 pmSubtractionMode mode ///< Mode for subtraction 110 113 ); … … 113 116 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) 114 117 int spatialOrder, ///< Order of spatial variations 118 float penalty, ///< Penalty for wideness 115 119 pmSubtractionMode mode ///< Mode for subtraction 116 120 ); … … 121 125 const psVector *fwhms, ///< Gaussian FWHMs 122 126 const psVector *orders, ///< Polynomial order of gaussians 127 float penalty, ///< Penalty for wideness 123 128 pmSubtractionMode mode ///< Mode for subtraction 124 129 ); … … 129 134 const psVector *fwhms, ///< Gaussian FWHMs 130 135 const psVector *orders, ///< Polynomial order of gaussians 136 float penalty, ///< Penalty for wideness 131 137 pmSubtractionMode mode ///< Mode for subtraction 132 138 ); … … 137 143 int inner, ///< Inner radius to preserve unbinned 138 144 int binning, ///< Kernel binning factor 145 float penalty, ///< Penalty for wideness 139 146 pmSubtractionMode mode ///< Mode for subtraction 140 147 ); … … 144 151 int spatialOrder, ///< Order of spatial variations 145 152 int inner, ///< Inner radius to preserve unbinned 153 float penalty, ///< Penalty for wideness 146 154 pmSubtractionMode mode ///< Mode for subtraction 147 155 ); … … 153 161 const psVector *orders, ///< Polynomial order of gaussians 154 162 int inner, ///< Inner radius containing grid of delta functions 163 float penalty, ///< Penalty for wideness 155 164 pmSubtractionMode mode ///< Mode for subtraction 156 165 ); … … 161 170 int inner, ///< Inner radius to preserve unbinned 162 171 int ringsOrder, ///< Polynomial order 172 float penalty, ///< Penalty for wideness 163 173 pmSubtractionMode mode ///< Mode for subtraction 164 174 ); … … 174 184 int binning, ///< Kernel binning factor 175 185 int ringsOrder, ///< Polynomial order for RINGS 186 float penalty, ///< Penalty for wideness 176 187 pmSubtractionMode mode ///< Mode for subtraction 177 188 ); -
trunk/psModules/src/imcombine/pmSubtractionMatch.c
r18202 r18287 117 117 pmSubtractionKernelsType type, int size, int spatialOrder, 118 118 const psVector *isisWidths, const psVector *isisOrders, 119 int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, 120 int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, 121 psMaskType maskBlank, float badFrac, pmSubtractionMode mode) 119 int inner, int ringsOrder, int binning, float penalty, 120 bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, 121 int iter, float rej, psMaskType maskBad, psMaskType maskBlank, 122 float badFrac, pmSubtractionMode mode) 122 123 { 123 124 if (mode != PM_SUBTRACTION_MODE_2) { … … 320 321 if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { 321 322 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, 322 stamps, footprint, optThreshold, mode);323 stamps, footprint, optThreshold, penalty, mode); 323 324 if (!kernels) { 324 325 psErrorClear(); … … 329 330 // Not an ISIS/GUNK kernel, or the optimum kernel search failed 330 331 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, 331 inner, binning, ringsOrder, mode);332 inner, binning, ringsOrder, penalty, mode); 332 333 } 333 334 … … 450 451 // Generate image with convolution kernels 451 452 int fullSize = 2 * size + 1 + 1; // Full size of kernel 452 psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32); 453 int imageSize = (2 * KERNEL_MOSAIC + 1) * fullSize; 454 psImage *convKernels = psImageAlloc((mode == PM_SUBTRACTION_MODE_DUAL ? 2 : 1) * 455 imageSize - 1 + 456 (mode == PM_SUBTRACTION_MODE_DUAL ? 4 : 0), 457 imageSize - 1, PS_TYPE_F32); 453 458 psImageInit(convKernels, NAN); 454 459 for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { … … 471 476 } 472 477 psFree(kernel); 478 479 if (mode == PM_SUBTRACTION_MODE_DUAL) { 480 kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, 481 (float)j / (float)KERNEL_MOSAIC, 482 true); // Image of the kernel 483 if (!kernel) { 484 psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); 485 psFree(convKernels); 486 goto MATCH_ERROR; 487 } 488 489 if (psImageOverlaySection(convKernels, kernel, 490 (2 * KERNEL_MOSAIC + 1 + i + KERNEL_MOSAIC) * fullSize + 491 4, 492 (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { 493 psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); 494 psFree(kernel); 495 psFree(convKernels); 496 goto MATCH_ERROR; 497 } 498 psFree(kernel); 499 } 500 473 501 } 474 502 } -
trunk/psModules/src/imcombine/pmSubtractionMatch.h
r18181 r18287 37 37 int ringsOrder, ///< RINGS polynomial order 38 38 int binning, ///< SPAM kernel binning 39 float penalty, ///< Penalty for wideness 39 40 bool optimum, ///< Search for optimum ISIS kernel? 40 41 const psVector *optFWHMs, ///< FWHMs for optimum search -
trunk/psModules/src/imcombine/pmSubtractionParams.c
r15756 r18287 204 204 int spatialOrder, const psVector *fwhms, int maxOrder, 205 205 const pmSubtractionStampList *stamps, int footprint, 206 float tolerance, pmSubtractionMode mode)206 float tolerance, float penalty, pmSubtractionMode mode) 207 207 { 208 208 if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { … … 231 231 psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders 232 232 psVectorInit(orders, maxOrder); 233 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, 234 fwhms, orders, mode); // Kernels233 pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, 234 penalty, mode); // Kernels 235 235 psFree(orders); 236 236 psFree(kernels->description); -
trunk/psModules/src/imcombine/pmSubtractionParams.h
r15443 r18287 16 16 int footprint, ///< Convolution footprint for stamps 17 17 float tolerance, ///< Maximum difference in chi^2 18 float penalty, ///< Penalty for wideness 18 19 pmSubtractionMode mode // Mode for subtraction 19 20 );
Note:
See TracChangeset
for help on using the changeset viewer.
