Changeset 29165
- Timestamp:
- Sep 15, 2010, 5:32:21 PM (16 years ago)
- Location:
- branches/eam_branches/ipp-20100823/psModules/src/imcombine
- Files:
-
- 5 edited
-
pmSubtractionEquation.c (modified) (38 diffs)
-
pmSubtractionEquation.h (modified) (1 diff)
-
pmSubtractionMatch.c (modified) (14 diffs)
-
pmSubtractionStamps.c (modified) (2 diffs)
-
pmSubtractionStamps.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c
r29148 r29165 27 27 28 28 // Calculate the least-squares matrix and vector 29 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated30 psVector *vector, // Least-squares vector, updated31 double *norm, // Normalisation, updated29 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated 30 psVector *vector, // Least-squares vector, updated 31 double normValue, // Normalisation, supplied 32 32 const psKernel *input, // Input image (target) 33 33 const psKernel *reference, // Reference image (convolution source) … … 37 37 const pmSubtractionKernels *kernels, // Kernels 38 38 const psImage *polyValues, // Spatial polynomial values 39 int footprint, // (Half-)Size of stamp 40 int normWindow1, // Window (half-)size for normalisation measurement 41 int normWindow2, // Window (half-)size for normalisation measurement 42 const pmSubtractionEquationCalculationMode mode 39 int footprint // (Half-)Size of stamp 43 40 ) 44 41 { … … 50 47 51 48 int numKernels = kernels->num; // Number of kernels 52 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation53 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background54 49 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 55 50 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 56 51 double poly[numPoly]; // Polynomial terms 57 52 double poly2[numPoly][numPoly]; // Polynomial-polynomial values 53 int numParams = numKernels * numPoly; 54 55 psAssert(matrix && 56 matrix->type.type == PS_TYPE_F64 && 57 matrix->numCols == numParams && 58 matrix->numRows == numParams, 59 "Least-squares matrix is bad."); 60 psAssert(vector && 61 vector->type.type == PS_TYPE_F64 && 62 vector->n == numParams, 63 "Least-squares vector is bad."); 58 64 59 65 // Evaluate polynomial-polynomial terms 60 // XXX we can skip this if we are not calculating kernel coeffs61 66 for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { 62 67 for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { … … 84 89 // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] 85 90 // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] 86 // normalization87 // bg 0, bg 1, bg 2 (only 0 is currently used?)88 91 89 92 for (int i = 0; i < numKernels; i++) { … … 107 110 108 111 // Spatial variation of kernel coeffs 109 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 110 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 111 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 112 double value = sumCC * poly2[iTerm][jTerm]; 113 matrix->data.F64[iIndex][jIndex] = value; 114 matrix->data.F64[jIndex][iIndex] = value; 115 } 116 } 117 } 112 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 113 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 114 double value = sumCC * poly2[iTerm][jTerm]; 115 matrix->data.F64[iIndex][jIndex] = value; 116 matrix->data.F64[jIndex][iIndex] = value; 117 } 118 } 118 119 } 119 120 120 121 double sumRC = 0.0; // Sum of the reference-convolution products 121 122 double sumIC = 0.0; // Sum of the input-convolution products 122 double sumC = 0.0; // Sum of the convolution123 123 for (int y = - footprint; y <= footprint; y++) { 124 124 for (int x = - footprint; x <= footprint; x++) { … … 128 128 double ic = in * conv; 129 129 double rc = ref * conv; 130 double c = conv;131 130 if (weight) { 132 131 float wtVal = weight->kernel[y][x]; 133 132 ic *= wtVal; 134 133 rc *= wtVal; 135 c *= wtVal;136 134 } 137 135 if (window) { … … 139 137 ic *= winVal; 140 138 rc *= winVal; 141 c *= winVal;142 139 } 143 140 sumIC += ic; 144 141 sumRC += rc; 145 sumC += c;146 147 } 142 } 143 } 144 148 145 // Spatial variation 149 146 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 150 double normTerm = sumRC * poly[iTerm]; 151 double bgTerm = sumC * poly[iTerm]; 152 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 153 matrix->data.F64[iIndex][normIndex] = normTerm; 154 matrix->data.F64[normIndex][iIndex] = normTerm; 155 } 156 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 157 matrix->data.F64[iIndex][bgIndex] = bgTerm; 158 matrix->data.F64[bgIndex][iIndex] = bgTerm; 159 } 160 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 161 vector->data.F64[iIndex] = sumIC * poly[iTerm]; 162 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 163 // subtract norm * sumRC * poly[iTerm] 164 psAssert (kernels->solution1, "programming error: define solution first!"); 165 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 166 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 167 vector->data.F64[iIndex] -= norm * normTerm; 168 } 169 } 170 } 171 } 172 173 double sumRR = 0.0; // Sum of the reference product 174 double sumIR = 0.0; // Sum of the input-reference product 175 double sum1 = 0.0; // Sum of the background 176 double sumR = 0.0; // Sum of the reference 177 double sumI = 0.0; // Sum of the input 178 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 179 for (int y = - footprint; y <= footprint; y++) { 180 for (int x = - footprint; x <= footprint; x++) { 181 double in = input->kernel[y][x]; 182 double ref = reference->kernel[y][x]; 183 double ir = in * ref; 184 double rr = PS_SQR(ref); 185 double one = 1.0; 186 187 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 188 normI1 += ref; 189 } 190 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 191 normI2 += in; 192 } 193 194 if (weight) { 195 float wtVal = weight->kernel[y][x]; 196 rr *= wtVal; 197 ir *= wtVal; 198 in *= wtVal; 199 ref *= wtVal; 200 one *= wtVal; 201 } 202 if (window) { 203 float winVal = window->kernel[y][x]; 204 rr *= winVal; 205 ir *= winVal; 206 in *= winVal; 207 ref *= winVal; 208 one *= winVal; 209 } 210 sumRR += rr; 211 sumIR += ir; 212 sumR += ref; 213 sumI += in; 214 sum1 += one; 215 } 216 } 217 218 *norm = normI2 / normI1; 219 220 fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); 221 222 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 223 matrix->data.F64[normIndex][normIndex] = sumRR; 224 vector->data.F64[normIndex] = sumIR; 225 // subtract sum over kernels * kernel solution 226 } 227 if (mode & PM_SUBTRACTION_EQUATION_BG) { 228 matrix->data.F64[bgIndex][bgIndex] = sum1; 229 vector->data.F64[bgIndex] = sumI; 230 } 231 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 232 matrix->data.F64[normIndex][bgIndex] = sumR; 233 matrix->data.F64[bgIndex][normIndex] = sumR; 147 vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm]; 148 } 234 149 } 235 150 … … 255 170 256 171 // Calculate the least-squares matrix and vector for dual convolution 257 static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated258 psVector *vector, // Least-squares vector, updated259 double *norm,// Normalisation, updated172 static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated 173 psVector *vector, // Least-squares vector, updated 174 double normValue, // Normalisation, updated 260 175 const psKernel *image1, // Image 1 261 176 const psKernel *image2, // Image 2 … … 266 181 const pmSubtractionKernels *kernels, // Kernels 267 182 const psImage *polyValues, // Spatial polynomial values 268 int footprint, // (Half-)Size of stamp 269 int normWindow1, // Window (half-)size for normalisation measurement 270 int normWindow2, // Window (half-)size for normalisation measurement 271 const pmSubtractionEquationCalculationMode mode 183 int footprint // (Half-)Size of stamp 272 184 ) 273 185 { 274 186 int numKernels = kernels->num; // Number of kernels 275 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation276 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background277 187 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 278 188 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms … … 280 190 double poly2[numPoly][numPoly]; // Polynomial-polynomial values 281 191 282 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 283 int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters 284 int numParams2 = numKernels * numPoly; // Number of additional parameters for dual 285 int numDual = numParams + numParams2; // Total number of parameters for dual 192 int numParams = numKernels * numPoly; // Number of regular parameters 193 int numParams2 = numKernels * numPoly; // Number of additional parameters for dual 194 int numDual = numParams + numParams2; // Total number of parameters for dual 286 195 287 196 psAssert(matrix && … … 317 226 matrix->data.F64[i][i] = 1.0; 318 227 } 228 229 // the order of the elements in the matrix and vector is: 230 // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] 231 // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] 232 // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] 319 233 320 234 for (int i = 0; i < numKernels; i++) { … … 352 266 353 267 // Spatial variation of kernel coeffs 354 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {355 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {356 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 357 double aa = sumAA* poly2[iTerm][jTerm];358 double bb = sumBB * poly2[iTerm][jTerm];359 double ab = sumAB * poly2[iTerm][jTerm]; 360 361 matrix->data.F64[iIndex][jIndex] = aa;362 matrix->data.F64[jIndex][iIndex] = aa; 363 364 matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb;365 matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; 366 367 matrix->data.F64[iIndex][jIndex + numParams] = ab;368 matrix->data.F64[jIndex + numParams][iIndex] = ab; 369 370 371 } 372 } 268 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 269 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 270 double aa = sumAA * poly2[iTerm][jTerm]; 271 double bb = sumBB * poly2[iTerm][jTerm]; 272 double ab = sumAB * poly2[iTerm][jTerm]; 273 274 matrix->data.F64[iIndex][jIndex] = aa; 275 matrix->data.F64[jIndex][iIndex] = aa; 276 277 matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; 278 matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; 279 280 matrix->data.F64[iIndex][jIndex + numParams] = ab; 281 matrix->data.F64[jIndex + numParams][iIndex] = ab; 282 } 283 } 284 } 285 286 // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B 373 287 for (int j = 0; j < i; j++) { 374 288 psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j … … 388 302 389 303 // Spatial variation of kernel coeffs 390 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 391 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 392 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 393 double ab = sumAB * poly2[iTerm][jTerm]; 394 matrix->data.F64[iIndex][jIndex + numParams] = ab; 395 matrix->data.F64[jIndex + numParams][iIndex] = ab; 396 } 397 } 304 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 305 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 306 double ab = sumAB * poly2[iTerm][jTerm]; 307 matrix->data.F64[iIndex][jIndex + numParams] = ab; 308 matrix->data.F64[jIndex + numParams][iIndex] = ab; 309 } 398 310 } 399 311 } … … 402 314 double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) 403 315 double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) 404 double sumA = 0.0; // Sum of A (for matrix, background)405 316 double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) 406 double sumB = 0.0; // Sum of B products (for matrix, background)407 double sumI2 = 0.0; // Sum of I_2 (for vector, background)408 317 for (int y = - footprint; y <= footprint; y++) { 409 318 for (int x = - footprint; x <= footprint; x++) { … … 424 333 ai1 *= wtVal; 425 334 bi1 *= wtVal; 426 a *= wtVal;427 b *= wtVal;428 i2 *= wtVal;429 335 } 430 336 if (window) { … … 434 340 ai1 *= wtVal; 435 341 bi1 *= wtVal; 436 a *= wtVal;437 b *= wtVal;438 i2 *= wtVal;439 342 } 440 343 sumAI2 += ai2; 441 344 sumBI2 += bi2; 442 345 sumAI1 += ai1; 443 sumA += a;444 346 sumBI1 += bi1; 445 sumB += b;446 sumI2 += i2;447 347 } 448 348 } … … 452 352 double bi2 = sumBI2 * poly[iTerm]; 453 353 double ai1 = sumAI1 * poly[iTerm]; 454 double a = sumA * poly[iTerm];455 354 double bi1 = sumBI1 * poly[iTerm]; 456 double b = sumB * poly[iTerm]; 457 458 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 459 matrix->data.F64[iIndex][normIndex] = ai1; 460 matrix->data.F64[normIndex][iIndex] = ai1; 461 matrix->data.F64[iIndex + numParams][normIndex] = bi1; 462 matrix->data.F64[normIndex][iIndex + numParams] = bi1; 463 } 464 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 465 matrix->data.F64[iIndex][bgIndex] = a; 466 matrix->data.F64[bgIndex][iIndex] = a; 467 matrix->data.F64[iIndex + numParams][bgIndex] = b; 468 matrix->data.F64[bgIndex][iIndex + numParams] = b; 469 } 470 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 471 vector->data.F64[iIndex] = ai2; 472 vector->data.F64[iIndex + numParams] = bi2; 473 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 474 // subtract norm * sumRC * poly[iTerm] 475 psAssert (kernels->solution1, "programming error: define solution first!"); 476 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 477 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 478 vector->data.F64[iIndex] -= norm * ai1; 479 vector->data.F64[iIndex + numParams] -= norm * bi1; 480 } 481 } 482 } 483 } 484 485 double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) 486 double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) 487 double sum1 = 0.0; // Sum of 1 (for matrix, background-background) 488 double sumI2 = 0.0; // Sum of I_2 (for vector, background) 489 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) 490 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 491 for (int y = - footprint; y <= footprint; y++) { 492 for (int x = - footprint; x <= footprint; x++) { 493 double i1 = image1->kernel[y][x]; 494 double i2 = image2->kernel[y][x]; 495 496 double i1i1 = i1 * i1; 497 double one = 1.0; 498 double i1i2 = i1 * i2; 499 500 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 501 normI1 += i1; 502 } 503 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 504 normI2 += i2; 505 } 506 507 if (weight) { 508 float wtVal = weight->kernel[y][x]; 509 i1 *= wtVal; 510 i1i1 *= wtVal; 511 one *= wtVal; 512 i2 *= wtVal; 513 i1i2 *= wtVal; 514 } 515 if (window) { 516 float wtVal = window->kernel[y][x]; 517 i1 *= wtVal; 518 i1i1 *= wtVal; 519 one *= wtVal; 520 i2 *= wtVal; 521 i1i2 *= wtVal; 522 } 523 sumI1 += i1; 524 sumI1I1 += i1i1; 525 sum1 += one; 526 sumI2 += i2; 527 sumI1I2 += i1i2; 528 } 529 } 530 531 *norm = normI2 / normI1; 532 fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); 533 534 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 535 matrix->data.F64[normIndex][normIndex] = sumI1I1; 536 vector->data.F64[normIndex] = sumI1I2; 537 } 538 if (mode & PM_SUBTRACTION_EQUATION_BG) { 539 matrix->data.F64[bgIndex][bgIndex] = sum1; 540 vector->data.F64[bgIndex] = sumI2; 541 } 542 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 543 matrix->data.F64[bgIndex][normIndex] = sumI1; 544 matrix->data.F64[normIndex][bgIndex] = sumI1; 355 vector->data.F64[iIndex] = ai2 - normValue * ai1; 356 vector->data.F64[iIndex + numParams] = bi2 - normValue * bi1; 357 } 545 358 } 546 359 … … 565 378 } 566 379 567 #if 1568 380 // Add in penalty term to least-squares vector 569 381 bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term 570 382 psVector *vector, // Vector to which to add in penalty term 571 383 const pmSubtractionKernels *kernels, // Kernel parameters 572 float norm // Normalisation 384 float normSquare1, // Normalisation for image 1 385 float normSquare2 // Normalisation for image 2 573 386 ) 574 387 { 388 psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution"); 389 575 390 if (kernels->penalty == 0.0) { 576 391 return true; … … 589 404 // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] 590 405 // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] 591 // [norm] 592 // [bg] 406 593 407 // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] 594 408 // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] … … 600 414 // Contribution to chi^2: a_i^2 P_i 601 415 psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); 602 fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); 603 matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; 604 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 605 fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 606 matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); 607 matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; 608 // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) 609 // penalties scale with second moments 610 // 611 } 612 } 613 } 614 } 615 616 return true; 617 } 618 # endif 416 fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]); 417 matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i]; 418 419 fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 420 matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]); 421 matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i]; 422 // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) 423 // penalties scale with second moments 424 // 425 } 426 } 427 } 428 429 return true; 430 } 619 431 620 432 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 647 459 int index, bool wantDual) 648 460 { 649 #if 0650 461 // This is probably in a tight loop, so don't check inputs 651 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);652 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);653 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);654 PS_ASSERT_INT_POSITIVE(index, NAN);655 #endif656 657 462 psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector 658 463 return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, … … 662 467 double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) 663 468 { 664 #if 0665 469 // This is probably in a tight loop, so don't check inputs 666 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);667 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);668 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);669 #endif670 671 470 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 672 471 return kernels->solution1->data.F64[normIndex]; … … 676 475 const psImage *polyValues) 677 476 { 678 #if 0679 477 // This is probably in a tight loop, so don't check inputs 680 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);681 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);682 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);683 #endif684 685 478 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 686 479 return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); … … 690 483 // Public functions 691 484 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 485 486 bool pmSubtractionCalculateNormalization( 487 pmSubtractionStampList *stamps, 488 const pmSubtractionMode mode) 489 { 490 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 491 492 psTimerStart("pmSubtractionCalculateNormalization"); 493 494 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 495 496 int footprint = stamps->footprint; // Half-size of stamps 497 498 // Loop over each stamp and calculate its normalization factor 499 for (int i = 0; i < stamps->num; i++) { 500 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 501 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; 502 if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; 503 504 // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic) 505 // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this 506 if (mode == PM_SUBTRACTION_MODE_2) { 507 pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1); 508 } else { 509 pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2); 510 } 511 psVectorAppend(norms, stamp->norm); 512 } 513 514 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 515 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 516 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 517 psFree(stats); 518 psFree(norms); 519 return false; 520 } 521 stamps->normValue = stats->robustMedian; 522 523 fprintf (stderr, "norm (1): %f\n", stamps->normValue); 524 525 psFree(stats); 526 psFree(norms); 527 528 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization")); 529 530 return true; 531 } 532 533 bool pmSubtractionCalculateNormalizationStamp( 534 pmSubtractionStamp *stamp, // stamp on which to save normalization) 535 const psKernel *image1, // Input image (target) 536 const psKernel *image2, // Reference image (convolution source) 537 int footprint, // (Half-)Size of stamp 538 int normWindow1, // Window (half-)size for normalisation measurement 539 int normWindow2 // Window (half-)size for normalisation measurement 540 ) 541 { 542 double normI1 = 0.0; // Sum of I_1 within the normalisation window (aperture) 543 double normI2 = 0.0; // Sum of I_2 within the normalisation window (aperture) 544 double normSquare1 = 0.0; // Sum of (I_1)^2 within the normalisation window (aperture) 545 double normSquare2 = 0.0; // Sum of (I_2)^2 within the normalisation window (aperture) 546 for (int y = - footprint; y <= footprint; y++) { 547 for (int x = - footprint; x <= footprint; x++) { 548 double im1 = image1->kernel[y][x]; 549 double im2 = image2->kernel[y][x]; 550 551 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 552 normI1 += im1; 553 normSquare1 += PS_SQR(im1); 554 } 555 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 556 normI2 += im2; 557 normSquare2 += PS_SQR(im2); 558 } 559 } 560 } 561 stamp->norm = normI2 / normI1; 562 stamp->normSquare1 = normSquare1; 563 stamp->normSquare2 = normSquare2; 564 565 fprintf (stderr, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); 566 567 return true; 568 } 692 569 693 570 bool pmSubtractionCalculateEquationThread(psThreadJob *job) … … 715 592 int numKernels = kernels->num; // Number of kernel basis functions 716 593 int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations 717 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms718 594 719 595 // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). … … 722 598 // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the 723 599 // number of coefficients for the spatial polynomial, normalisation and a constant background offset. 724 int numParams = numKernels * numSpatial + 1 + numBackground; 600 // XXX we no longer solve for the normalization and background in the matrix inversion 601 int numParams = numKernels * numSpatial; 725 602 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 726 603 // An additional image is convolved … … 790 667 switch (kernels->mode) { 791 668 case PM_SUBTRACTION_MODE_1: 792 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, 793 weight, window, stamp->convolutions1, kernels, 794 polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); 669 status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1, 670 weight, window, stamp->convolutions1, kernels, polyValues, footprint); 795 671 break; 796 672 case PM_SUBTRACTION_MODE_2: 797 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, 798 weight, window, stamp->convolutions2, kernels, 799 polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); 673 status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, 674 weight, window, stamp->convolutions2, kernels, polyValues, footprint); 800 675 break; 801 676 case PM_SUBTRACTION_MODE_DUAL: 802 status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, 803 stamp->image1, stamp->image2, 677 status = calculateDualMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, 804 678 weight, window, stamp->convolutions1, stamp->convolutions2, 805 kernels, polyValues, footprint , stamps->normWindow1, stamps->normWindow2, mode);679 kernels, polyValues, footprint); 806 680 break; 807 681 default: … … 928 802 int numKernels = kernels->num; // Number of kernel basis functions 929 803 int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations 930 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 931 int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for 932 int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution 804 // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 805 // XXX int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for 806 int numParams = numKernels * numSpatial; // Number of parameters being solved for 807 int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution 933 808 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 934 809 // An additional image is convolved … … 960 835 961 836 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { 962 // Accumulate the least-squares matricies and vectors 837 // Accumulate the least-squares matrices and vectors. These are generated for the 838 // kernel elements, excluding the background and normalization. 963 839 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix 964 840 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector 965 841 psVectorInit(sumVector, 0.0); 966 842 psImageInit(sumMatrix, 0.0); 967 968 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations969 843 970 844 int numStamps = 0; // Number of good stamps … … 976 850 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 977 851 978 psVectorAppend(norms, stamp->norm);979 980 852 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 981 853 numStamps++; … … 985 857 } 986 858 987 #if 0988 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background989 calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);990 #endif991 992 859 psVector *solution = NULL; // Solution to equation! 993 860 solution = psVectorAlloc(numParams, PS_TYPE_F64); 994 861 psVectorInit(solution, 0); 995 862 996 #if 0 997 // Regular, straight-forward solution 998 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 999 #else 1000 { 1001 // Solve normalisation and background separately 1002 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1003 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1004 1005 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 1006 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 1007 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 1008 psFree(stats); 1009 psFree(sumMatrix); 1010 psFree(sumVector); 1011 psFree(norms); 1012 return false; 1013 } 1014 1015 // double normValue = 1.0; 1016 double normValue = stats->robustMedian; 1017 // double bgValue = 0.0; 1018 1019 psFree(stats); 1020 1021 #ifdef TESTING 1022 fprintf(stderr, "Norm: %lf\n", normValue); 1023 #endif 1024 // Solve kernel components 1025 for (int i = 0; i < numSolution1; i++) { 1026 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1027 1028 sumMatrix->data.F64[i][normIndex] = 0.0; 1029 sumMatrix->data.F64[normIndex][i] = 0.0; 1030 } 1031 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1032 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1033 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1034 1035 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1036 sumVector->data.F64[normIndex] = 0.0; 1037 1038 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1039 1040 solution->data.F64[normIndex] = normValue; 1041 } 1042 # endif 1043 1044 #if (1) 863 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 864 1045 865 for (int i = 0; i < solution->n; i++) { 1046 fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); 1047 } 1048 #endif 866 psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]); 867 } 1049 868 1050 869 if (!kernels->solution1) { 1051 kernels->solution1 = psVectorAlloc(sumVector->n , PS_TYPE_F64);870 kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg 1052 871 psVectorInit(kernels->solution1, 0.0); 1053 872 } 1054 873 1055 // only update the solutions that we chose to calculate: 1056 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1057 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1058 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1059 } 1060 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1061 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1062 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1063 } 1064 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1065 int numKernels = kernels->num; 1066 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 1067 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 1068 for (int i = 0; i < numKernels * numPoly; i++) { 1069 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1070 } 1071 } 1072 1073 psFree(norms); 874 int numKernels = kernels->num; 875 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 876 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 877 for (int i = 0; i < numKernels * numPoly; i++) { 878 kernels->solution1->data.F64[i] = solution->data.F64[i]; 879 } 880 881 // Apply the normalisation and background separately 882 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 883 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 884 kernels->solution1->data.F64[normIndex] = stamps->normValue; 885 kernels->solution1->data.F64[bgIndex] = 0.0; 886 1074 887 psFree(solution); 1075 888 psFree(sumVector); 1076 889 psFree(sumMatrix); 1077 890 1078 #ifdef TESTING1079 // XXX double-check for NAN in data:1080 for (int ix = 0; ix < kernels->solution1->n; ix++) {1081 if (!isfinite(kernels->solution1->data.F64[ix])) {1082 fprintf (stderr, "WARNING: NAN in vector\n");1083 }1084 }1085 #endif1086 1087 891 } else { 1088 892 // Dual convolution solution 1089 1090 // Accumulation of stamp matrices/vectors 893 // Accumulate the least-squares matrices and vectors. These are generated for the 894 // kernel elements, excluding the background and normalization. 1091 895 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); 1092 896 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); … … 1094 898 psVectorInit(sumVector, 0.0); 1095 899 1096 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations1097 1098 int numStamps = 0; // Number of goodstamps900 int numStamps = 0; // Number of good stamps 901 double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps 902 double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps 1099 903 for (int i = 0; i < stamps->num; i++) { 1100 904 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 1103 907 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 1104 908 1105 psVectorAppend(norms, stamp->norm); 909 normSquare1 += stamp->normSquare1; 910 normSquare2 += stamp->normSquare2; 1106 911 1107 912 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 1108 913 numStamps++; 914 } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { 915 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); 1109 916 } 1110 917 } … … 1117 924 #endif 1118 925 1119 #if 1 1120 // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1121 // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); 1122 1123 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1124 calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0); 1125 #endif 926 calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2); 1126 927 1127 928 psVector *solution = NULL; // Solution to equation! … … 1129 930 psVectorInit(solution, 0); 1130 931 1131 #if 0 1132 // Regular, straight-forward solution 1133 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1134 #else 1135 { 1136 // Solve normalisation and background separately 1137 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1138 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1139 1140 #if 0 1141 psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); 1142 psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); 1143 1144 normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; 1145 normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; 1146 normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; 1147 1148 normVector->data.F64[0] = sumVector->data.F64[normIndex]; 1149 normVector->data.F64[1] = sumVector->data.F64[bgIndex]; 1150 1151 psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); 1152 1153 double normValue = normSolution->data.F64[0]; 1154 double bgValue = normSolution->data.F64[1]; 1155 1156 psFree(normMatrix); 1157 psFree(normVector); 1158 psFree(normSolution); 1159 #endif 1160 1161 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 1162 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 1163 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 1164 psFree(stats); 1165 psFree(sumMatrix); 1166 psFree(sumVector); 1167 psFree(norms); 1168 return false; 1169 } 1170 1171 double normValue = stats->robustMedian; 1172 1173 psFree(stats); 1174 1175 #ifdef TESTING 1176 fprintf(stderr, "Norm: %lf\n", normValue); 1177 #endif 1178 1179 // Solve kernel components 1180 for (int i = 0; i < numSolution2; i++) { 1181 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1182 sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; 1183 1184 sumMatrix->data.F64[i][normIndex] = 0.0; 1185 sumMatrix->data.F64[normIndex][i] = 0.0; 1186 1187 sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; 1188 sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; 1189 } 1190 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1191 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1192 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1193 1194 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1195 1196 sumVector->data.F64[normIndex] = 0.0; 1197 1198 // save the matrix and vector after the NULLs have been set 1199 #if 0 1200 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); 1201 psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); 1202 psVectorWriteFile("sumVector.dat", sumVector); 1203 psFree (save); 1204 #endif 1205 1206 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6); 1207 // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 3e-4); 1208 // psVectorCopy (solution, sumVector, PS_TYPE_F64); 1209 // psMatrixGJSolve(sumMatrix, solution); 1210 solution->data.F64[normIndex] = normValue; 1211 } 1212 #endif 1213 932 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6); 1214 933 1215 934 #if (1) … … 1219 938 #endif 1220 939 1221 psFree(sumMatrix);1222 psFree(sumVector);1223 1224 psFree(norms);1225 1226 940 if (!kernels->solution1) { 1227 kernels->solution1 = psVectorAlloc(numSolution1 , PS_TYPE_F64);941 kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); 1228 942 psVectorInit (kernels->solution1, 0.0); 1229 943 } … … 1233 947 } 1234 948 1235 // only update the solutions that we chose to calculate: 1236 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1237 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1238 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1239 } 1240 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1241 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1242 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1243 } 1244 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1245 int numKernels = kernels->num; 1246 for (int i = 0; i < numKernels * numSpatial; i++) { 1247 // XXX fprintf (stderr, "keep\n"); 1248 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1249 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; 1250 } 1251 } 1252 1253 1254 memcpy(kernels->solution1->data.F64, solution->data.F64, 1255 numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1256 memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], 1257 numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1258 949 int numKernels = kernels->num; 950 for (int i = 0; i < numKernels * numSpatial; i++) { 951 // XXX fprintf (stderr, "keep\n"); 952 kernels->solution1->data.F64[i] = solution->data.F64[i]; 953 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; 954 } 955 956 // Apply the normalisation and background separately 957 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 958 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 959 kernels->solution1->data.F64[normIndex] = stamps->normValue; 960 kernels->solution1->data.F64[bgIndex] = 0.0; 961 962 psFree(sumMatrix); 963 psFree(sumVector); 1259 964 psFree(solution); 1260 1261 965 } 1262 966 … … 1927 1631 return true; 1928 1632 } 1929 1930 1931 # if 01932 1933 #ifdef TESTING1934 psFitsWriteImageSimple("A.fits", sumMatrix, NULL);1935 psVectorWriteFile ("B.dat", sumVector);1936 #endif1937 1938 # define SVD_ANALYSIS 01939 # define COEFF_SIG 0.01940 # define SVD_TOL 0.01941 1942 // Use SVD to determine the kernel coeffs (and validate)1943 if (SVD_ANALYSIS) {1944 1945 // We have sumVector and sumMatrix. we are trying to solve the following equation:1946 // sumMatrix * x = sumVector.1947 1948 // we can use any standard matrix inversion to solve this. However, the basis1949 // functions in general have substantial correlation, so that the solution may be1950 // somewhat poorly determined or unstable. If not numerically ill-conditioned, the1951 // system of equations may be statistically ill-conditioned. Noise in the image1952 // will drive insignificant, but correlated, terms in the solution. To avoid these1953 // problems, we can use SVD to identify numerically unconstrained values and to1954 // avoid statistically badly determined value.1955 1956 // A = sumMatrix, B = sumVector1957 // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T1958 // x = V (1/w) (U^T B)1959 // \sigma_x = sqrt(diag(A^{-1}))1960 // solve for x and A^{-1} to get x & dx1961 // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.01962 // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.01963 1964 // If I use the SVD trick to re-condition the matrix, I need to break out the1965 // kernel and normalization terms from the background term.1966 // XXX is this true? or was this due to an error in the analysis?1967 1968 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background1969 1970 // now pull out the kernel elements into their own square matrix1971 psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);1972 psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);1973 1974 for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {1975 if (ix == bgIndex) continue;1976 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {1977 if (iy == bgIndex) continue;1978 kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];1979 ky++;1980 }1981 kernelVector->data.F64[kx] = sumVector->data.F64[ix];1982 kx++;1983 }1984 1985 psImage *U = NULL;1986 psImage *V = NULL;1987 psVector *w = NULL;1988 if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {1989 psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n");1990 return NULL;1991 }1992 1993 // calculate A_inverse:1994 // Ainv = V * w * U^T1995 psImage *wUt = p_pmSubSolve_wUt (w, U);1996 psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);1997 psImage *Xvar = NULL;1998 psFree (wUt);1999 2000 # ifdef TESTING2001 // kernel terms:2002 for (int i = 0; i < w->n; i++) {2003 fprintf (stderr, "w: %f\n", w->data.F64[i]);2004 }2005 # endif2006 // loop over w adding in more and more of the values until chisquare is no longer2007 // dropping significantly.2008 // XXX this does not seem to work very well: we seem to need all terms even for2009 // simple cases...2010 2011 psVector *Xsvd = NULL;2012 {2013 psVector *Ax = NULL;2014 psVector *UtB = NULL;2015 psVector *wUtB = NULL;2016 2017 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);2018 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);2019 psVectorInit (wMask, 1); // start by masking everything2020 2021 double chiSquareLast = NAN;2022 int maxWeight = 0;2023 2024 double Axx, Bx, y2;2025 2026 // XXX this is an attempt to exclude insignificant modes.2027 // it was not successful with the ISIS kernel set: removing even2028 // the least significant mode leaves additional ringing / noise2029 // because the terms are so coupled.2030 for (int k = 0; false && (k < w->n); k++) {2031 2032 // unmask the k-th weight2033 wMask->data.U8[k] = 0;2034 p_pmSubSolve_SetWeights(wApply, w, wMask);2035 2036 // solve for x:2037 // x = V * w * (U^T * B)2038 p_pmSubSolve_UtB (&UtB, U, kernelVector);2039 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);2040 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);2041 2042 // chi-square for this system of equations:2043 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^22044 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^22045 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);2046 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);2047 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);2048 p_pmSubSolve_y2 (&y2, kernels, stamps);2049 2050 // apparently, this works (compare with the brute force value below2051 double chiSquare = Axx - 2.0*Bx + y2;2052 double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;2053 chiSquareLast = chiSquare;2054 2055 // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);2056 if (k && !maxWeight && (deltaChi < 1.0)) {2057 maxWeight = k;2058 }2059 }2060 2061 // keep all terms or we get extra ringing2062 maxWeight = w->n;2063 psVectorInit (wMask, 1);2064 for (int k = 0; k < maxWeight; k++) {2065 wMask->data.U8[k] = 0;2066 }2067 p_pmSubSolve_SetWeights(wApply, w, wMask);2068 2069 // solve for x:2070 // x = V * w * (U^T * B)2071 p_pmSubSolve_UtB (&UtB, U, kernelVector);2072 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);2073 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);2074 2075 // chi-square for this system of equations:2076 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^22077 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^22078 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);2079 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);2080 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);2081 p_pmSubSolve_y2 (&y2, kernels, stamps);2082 2083 // apparently, this works (compare with the brute force value below2084 double chiSquare = Axx - 2.0*Bx + y2;2085 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);2086 2087 // re-calculate A^{-1} to get new variances:2088 // Ainv = V * w * U^T2089 // XXX since we keep all terms, this is identical to Ainv2090 psImage *wUt = p_pmSubSolve_wUt (wApply, U);2091 Xvar = p_pmSubSolve_VwUt (V, wUt);2092 psFree (wUt);2093 2094 psFree (Ax);2095 psFree (UtB);2096 psFree (wUtB);2097 psFree (wApply);2098 psFree (wMask);2099 }2100 2101 // copy the kernel solutions to the full solution vector:2102 solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);2103 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);2104 2105 for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {2106 if (ix == bgIndex) {2107 solution->data.F64[ix] = 0;2108 solutionErr->data.F64[ix] = 0.001;2109 continue;2110 }2111 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);2112 solution->data.F64[ix] = Xsvd->data.F64[kx];2113 kx++;2114 }2115 2116 psFree (kernelMatrix);2117 psFree (kernelVector);2118 2119 psFree (U);2120 psFree (V);2121 psFree (w);2122 2123 psFree (Ainv);2124 psFree (Xsvd);2125 } else {2126 psVector *permutation = NULL; // Permutation vector, required for LU decomposition2127 psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);2128 if (!luMatrix) {2129 psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n");2130 psFree(solution);2131 psFree(sumVector);2132 psFree(sumMatrix);2133 psFree(luMatrix);2134 psFree(permutation);2135 return NULL;2136 }2137 2138 solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);2139 psFree(luMatrix);2140 psFree(permutation);2141 if (!solution) {2142 psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n");2143 psFree(solution);2144 psFree(sumVector);2145 psFree(sumMatrix);2146 return NULL;2147 }2148 2149 // XXX LUD does not provide A^{-1}? fake the error for now2150 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);2151 for (int ix = 0; ix < sumVector->n; ix++) {2152 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];2153 }2154 }2155 2156 if (!kernels->solution1) {2157 kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);2158 psVectorInit (kernels->solution1, 0.0);2159 }2160 2161 // only update the solutions that we chose to calculate:2162 if (mode & PM_SUBTRACTION_EQUATION_NORM) {2163 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation2164 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];2165 }2166 if (mode & PM_SUBTRACTION_EQUATION_BG) {2167 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background2168 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];2169 }2170 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {2171 int numKernels = kernels->num;2172 int spatialOrder = kernels->spatialOrder; // Order of spatial variation2173 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms2174 for (int i = 0; i < numKernels * numPoly; i++) {2175 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));2176 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {2177 // XXX fprintf (stderr, "drop\n");2178 kernels->solution1->data.F64[i] = 0.0;2179 } else {2180 // XXX fprintf (stderr, "keep\n");2181 kernels->solution1->data.F64[i] = solution->data.F64[i];2182 }2183 }2184 }2185 // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);2186 // fprintf (stderr, "chi square Brute = %f\n", chiSquare);2187 2188 psFree(solution);2189 psFree(sumVector);2190 psFree(sumMatrix);2191 # endif2192 2193 #ifdef TESTING2194 // XXX double-check for NAN in data:2195 for (int iy = 0; iy < stamp->matrix->numRows; iy++) {2196 for (int ix = 0; ix < stamp->matrix->numCols; ix++) {2197 if (!isfinite(stamp->matrix->data.F64[iy][ix])) {2198 fprintf (stderr, "WARNING: NAN in matrix\n");2199 }2200 }2201 }2202 for (int ix = 0; ix < stamp->vector->n; ix++) {2203 if (!isfinite(stamp->vector->data.F64[ix])) {2204 fprintf (stderr, "WARNING: NAN in vector\n");2205 }2206 }2207 #endif2208 2209 #ifdef TESTING2210 for (int ix = 0; ix < sumVector->n; ix++) {2211 if (!isfinite(sumVector->data.F64[ix])) {2212 fprintf (stderr, "WARNING: NAN in vector\n");2213 }2214 }2215 #endif2216 2217 #ifdef TESTING2218 for (int ix = 0; ix < sumVector->n; ix++) {2219 if (!isfinite(sumVector->data.F64[ix])) {2220 fprintf (stderr, "WARNING: NAN in vector\n");2221 }2222 }2223 {2224 psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);2225 psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);2226 psFree(inverse);2227 }2228 {2229 psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);2230 psImage *Xt = psMatrixTranspose(NULL, X);2231 psImage *XtX = psMatrixMultiply(NULL, Xt, X);2232 psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);2233 psFree(X);2234 psFree(Xt);2235 psFree(XtX);2236 }2237 #endif2238 -
branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.h
r29004 r29165 65 65 ); 66 66 67 bool pmSubtractionCalculateNormalization( 68 pmSubtractionStampList *stamps, 69 const pmSubtractionMode mode); 70 71 bool pmSubtractionCalculateNormalizationStamp( 72 pmSubtractionStamp *stamp, // stamp on which to save normalization) 73 const psKernel *input, // Input image (target) 74 const psKernel *reference, // Reference image (convolution source) 75 int footprint, // (Half-)Size of stamp 76 int normWindow1, // Window (half-)size for normalisation measurement 77 int normWindow2 // Window (half-)size for normalisation measurement 78 ); 79 67 80 #endif -
branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c
r29004 r29165 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 … … 672 667 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, 673 668 inner, binning, ringsOrder, penalty, bounds, subMode); 674 // pmSubtractionVisualShowKernels(kernels);675 669 } 676 670 … … 678 672 679 673 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 674 pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); 704 675 switch (newMode) { … … 732 703 } 733 704 734 // XXX step 1: calculate normalization 735 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); 705 // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue 706 psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); 707 if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { 708 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); 709 goto MATCH_ERROR; 710 } 711 712 // step 1: generate the elements of the matrix equation Ax = B 713 psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); 736 714 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { 737 715 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); … … 739 717 } 740 718 741 psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n"); 719 // step 2: solve the matrix equation Ax = B 720 psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); 742 721 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { 743 722 psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); … … 746 725 memCheck(" solve equation"); 747 726 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 727 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations 768 728 if (!deviations) { … … 770 730 goto MATCH_ERROR; 771 731 } 772 773 732 memCheck(" calculate deviations"); 774 733 … … 787 746 // if we hit the max number of iterations and we have rejected stamps, re-solve 788 747 if (numRejected > 0) { 789 // XXX step 1: calculate normalization 748 749 // step 1: generate the elements of the matrix equation Ax = B 790 750 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); 791 751 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { … … 794 754 } 795 755 796 // s olve normalization756 // step 2: solve the matrix equation Ax = B 797 757 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); 798 758 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { … … 801 761 } 802 762 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 763 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations 823 764 if (!deviations) { … … 837 778 goto MATCH_ERROR; 838 779 } 839 840 780 memCheck("diag outputs"); 841 781 … … 1134 1074 } 1135 1075 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 1076 psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); 1154 1077 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations … … 1181 1104 } 1182 1105 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 1106 1202 1107 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations … … 1301 1206 PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); 1302 1207 1303 //float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference1208 // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference 1304 1209 float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor 1305 1210 -
branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c
r29148 r29165 244 244 list->flux = NULL; 245 245 list->normFrac = normFrac; 246 list->normValue = NAN; 246 247 list->window1 = NULL; 247 248 list->window2 = NULL; … … 343 344 stamp->vector = NULL; 344 345 stamp->norm = NAN; 346 stamp->normSquare1 = NAN; 347 stamp->normSquare2 = NAN; 345 348 346 349 return stamp; -
branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h
r29004 r29165 25 25 int footprint; ///< Half-size of stamps 26 26 float normFrac; ///< Fraction of flux in window for normalisation window 27 float normValue; ///< calculated normalization 27 28 psKernel *window1; ///< window function generated from ensemble of stamps (input 1) 28 29 psKernel *window2; ///< window function generated from ensemble of stamps (input 2) 29 float normWindow1; 30 float normWindow2; 30 float normWindow1; ///< Size of window for measuring normalisation 31 float normWindow2; ///< Size of window for measuring normalisation 31 32 float sysErr; ///< Systematic error 32 33 float skyErr; ///< increase effective readnoise … … 85 86 psVector *vector; ///< Least-squares vector, or NULL 86 87 double norm; ///< Normalisation difference 88 double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) 89 double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) 87 90 pmSubtractionStampStatus status; ///< Status of stamp 88 91 } pmSubtractionStamp;
Note:
See TracChangeset
for help on using the changeset viewer.
