Changeset 15809 for trunk/psModules/src/imcombine/pmSubtractionEquation.c
- Timestamp:
- Dec 13, 2007, 11:32:44 AM (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r15762 r15809 58 58 for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { 59 59 double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder]; 60 // XXX Not sure I've got this the right way around 60 61 assert(iIndex < matrix->numRows && j < matrix->numCols); 62 61 63 matrix->data.F64[iIndex][j] = convPoly; 62 64 if (symmetric) { 65 66 assert(iIndex < matrix->numCols && j < matrix->numRows); 67 63 68 matrix->data.F64[j][iIndex] = convPoly; 64 69 } … … 93 98 for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) { 94 99 double convPoly = sum * iPoly * polyValues->data.F64[jyOrder][jxOrder]; 100 101 assert(iIndex < matrix->numRows && jIndex < matrix->numCols); 102 95 103 matrix->data.F64[iIndex][jIndex] = convPoly; 96 104 if (symmetric) { 105 106 assert(iIndex < matrix->numCols && jIndex < matrix->numRows); 107 97 108 matrix->data.F64[jIndex][iIndex] = convPoly; 98 109 } … … 139 150 const psArray *convolutions, // Convolutions of source with kernels 140 151 const psKernel *input, // Input stamp, or NULL 141 const psKernel *target, // Target stamp, or NULL142 152 const psKernel *weight, // Weight stamp 143 153 const psImage *polyValues, // Spatial polynomial values … … 156 166 assert(matrix->numCols == numTerms); 157 167 assert(convolutions && convolutions->n == numKernels); 158 assert(target);159 168 assert(polyValues); 160 169 assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image … … 168 177 // XXX To support higher-order background model than simply constant, the below code needs to be updated. 169 178 if (normAndBG) { 170 // Normalisation and background terms 171 int normIndex, bgIndex; // Indices in matrix for normalisation and background 172 PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); 179 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 180 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 173 181 174 182 for (int i = 0; i < numKernels; i++) { 175 183 psKernel *conv = convolutions->data[i]; // Convolution for i-th element 176 // Normalisation 184 185 // Normalisation-convolution terms 177 186 if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, 178 187 footprint, spatialOrder, true)) { … … 181 190 } 182 191 183 // Background 192 // Background-convolution terms 184 193 double sumC = 0.0; // Sum of the convolution 185 194 for (int y = - footprint; y <= footprint; y++) { … … 192 201 return false; 193 202 } 203 194 204 for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { 195 205 for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { 196 matrix->data.F64[index][bgIndex] = matrix->data.F64[bgIndex][index] = 197 sumC * polyValues->data.F64[yOrder][xOrder]; 206 double value = sumC * polyValues->data.F64[yOrder][xOrder]; 207 matrix->data.F64[index][bgIndex] = value; 208 matrix->data.F64[bgIndex][index] = value; 198 209 } 199 210 } 200 211 } 201 212 202 // Normalisation only and background only terms in the matrix213 // Background only, normalisation only, and background-normalisation terms 203 214 double sum1 = 0.0; // Sum of the weighting 204 215 double sumI = 0.0; // Sum of the input … … 206 217 for (int y = - footprint; y <= footprint; y++) { 207 218 for (int x = - footprint; x <= footprint; x++) { 208 floatinvNoise2 = 1.0 / weight->kernel[y][x];209 floatvalue = input->kernel[y][x] * invNoise2;219 double invNoise2 = 1.0 / weight->kernel[y][x]; 220 double value = input->kernel[y][x] * invNoise2; 210 221 sumI += value; 211 sumII += input->kernel[y][x] * value;222 sumII += value * input->kernel[y][x]; 212 223 sum1 += invNoise2; 213 224 } 225 } 226 if (!isfinite(sumI)) { 227 psTrace("psModules.imcombine", 2, "Bad sumI detected"); 228 return false; 214 229 } 215 230 if (!isfinite(sumII)) { 216 231 psTrace("psModules.imcombine", 2, "Bad sumII detected"); 217 return false;218 }219 if (!isfinite(sumI)) {220 psTrace("psModules.imcombine", 2, "Bad sumI detected");221 232 return false; 222 233 } … … 239 250 const pmSubtractionKernels *kernels, // Kernel components 240 251 const psArray *convolutions, // Convolutions of source with kernels 241 const psKernel *input, // Input stamp, or NULL 242 const psKernel *target, // Target stamp , or NULL252 const psKernel *input, // Input stamp, or NULL if !normAndBG 253 const psKernel *target, // Target stamp 243 254 const psKernel *weight, // Weight stamp 244 255 const psImage *polyValues, // Spatial polynomial values 245 int footprint // (Half-)Size of stamp 256 int footprint, // (Half-)Size of stamp 257 bool normAndBG // Calculate normalisation and background terms? 246 258 ) 247 259 { … … 250 262 int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms 251 263 int bgOrder = kernels->bgOrder; // Maximum order of background fit 252 int numBackground = PM_SUBTRACTION_POLYTERMS(bgOrder); // Number of background terms253 int numTerms = numKernels * numSpatial + 1 + numBackground; // Total number of terms264 int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms 265 int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms 254 266 assert(vector && vector->n == numTerms); 255 267 assert(convolutions && convolutions->n == numKernels); 256 assert(input);257 268 assert(target); 258 269 assert(polyValues); 270 assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image 259 271 260 272 // Convolution terms … … 278 290 } 279 291 280 // Normalisation and background terms 281 double sumIT = 0.0; // Sum of the target and input 282 double sumT = 0.0; // Sum of the target 283 for (int y = - footprint; y <= footprint; y++) { 284 for (int x = - footprint; x <= footprint; x++) { 285 float value = target->kernel[y][x] / weight->kernel[y][x]; 286 sumIT += value * input->kernel[y][x]; 287 sumT += value; 288 } 289 } 290 if (!isfinite(sumIT)) { 291 psTrace("psModules.imcombine", 2, "Bad sumIT detected"); 292 return false; 293 } 294 if (!isfinite(sumT)) { 295 psTrace("psModules.imcombine", 2, "Bad sumI detected"); 296 return false; 297 } 298 299 int normIndex, bgIndex; // Indices in vector for normalisation and background terms 300 PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); 301 vector->data.F64[normIndex] = sumIT; 302 vector->data.F64[bgIndex] = sumT; 292 if (normAndBG) { 293 // Background terms 294 double sumT = 0.0; // Sum of the target 295 double sumIT = 0.0; // Sum of the input-target product 296 for (int y = - footprint; y <= footprint; y++) { 297 for (int x = - footprint; x <= footprint; x++) { 298 float value = target->kernel[y][x] / weight->kernel[y][x]; 299 sumIT += value * input->kernel[y][x]; 300 sumT += value; 301 } 302 } 303 if (!isfinite(sumT)) { 304 psTrace("psModules.imcombine", 2, "Bad sumI detected"); 305 return false; 306 } 307 if (!isfinite(sumIT)) { 308 psTrace("psModules.imcombine", 2, "Bad sumIT detected"); 309 return false; 310 } 311 312 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation term 313 vector->data.F64[normIndex] = sumIT; 314 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background term 315 vector->data.F64[bgIndex] = sumT; 316 } 303 317 304 318 return true; … … 383 397 for (int yOrder = 0; yOrder <= order; yOrder++) { 384 398 for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) { 399 400 assert(index < coeff->n); 401 385 402 sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; 386 403 } … … 468 485 } 469 486 487 #ifdef TESTING 488 for (int j = 0; j < numKernels; j++) { 489 if (stamp->convolutions1) { 490 psString convName = NULL; 491 psStringAppend(&convName, "conv1_%03d_%03d.fits", i, j); 492 psFits *fits = psFitsOpen(convName, "w"); 493 psFree(convName); 494 psKernel *conv = stamp->convolutions1->data[j]; 495 psFitsWriteImage(fits, NULL, conv->image, 0, NULL); 496 psFitsClose(fits); 497 } 498 499 if (stamp->convolutions2) { 500 psString convName = NULL; 501 psStringAppend(&convName, "conv2_%03d_%03d.fits", i, j); 502 psFits *fits = psFitsOpen(convName, "w"); 503 psFree(convName); 504 psKernel *conv = stamp->convolutions2->data[j]; 505 psFitsWriteImage(fits, NULL, conv->image, 0, NULL); 506 psFitsClose(fits); 507 } 508 } 509 #endif 510 470 511 polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm); 471 512 472 stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); 513 514 stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); 515 stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64); 473 516 #ifdef TESTING 474 psVectorInit(stamp->vector, NAN); 517 psImageInit(stamp->matrix1, NAN); 518 psVectorInit(stamp->vector1, NAN); 475 519 #endif 476 520 … … 479 523 case PM_SUBTRACTION_MODE_TARGET: 480 524 case PM_SUBTRACTION_MODE_1: 481 stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);482 #ifdef TESTING483 psImageInit(stamp->matrix1, NAN);484 #endif485 525 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 486 stamp-> image2, stamp->weight, polyValues, footprint, true);487 status &= calculateVector(stamp->vector , kernels, stamp->convolutions1, stamp->image1,488 stamp->image2, stamp->weight, polyValues, footprint );526 stamp->weight, polyValues, footprint, true); 527 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 528 stamp->image2, stamp->weight, polyValues, footprint, true); 489 529 break; 490 530 case PM_SUBTRACTION_MODE_2: 491 stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);492 #ifdef TESTING493 psImageInit(stamp->matrix1, NAN);494 #endif495 531 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, 496 stamp-> image1, stamp->weight, polyValues, footprint, true);497 status &= calculateVector(stamp->vector , kernels, stamp->convolutions2, stamp->image2,498 stamp->image1, stamp->weight, polyValues, footprint );532 stamp->weight, polyValues, footprint, true); 533 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, 534 stamp->image1, stamp->weight, polyValues, footprint, true); 499 535 break; 500 536 case PM_SUBTRACTION_MODE_DUAL: 501 stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);502 537 stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); 503 538 stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64); 539 stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); 504 540 #ifdef TESTING 505 psImageInit(stamp->matrix1, NAN);506 541 psImageInit(stamp->matrix2, NAN); 507 542 psImageInit(stamp->matrixX, NAN); 508 #endif 509 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 510 stamp->image2, stamp->weight, polyValues, footprint, true);511 status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, stamp->image2,512 stamp->image1, stamp->weight, polyValues, footprint, false);513 // XXX Not sure I'm passing the correct stamp->image to calculateMatrixCross:543 psVectorInit(stamp->vector2, NAN); 544 #endif 545 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 546 stamp->weight, polyValues, footprint, true); 547 status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, 548 stamp->weight, polyValues, footprint, false); 514 549 status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, 515 stamp->convolutions2, stamp->image 2, stamp->weight, polyValues,550 stamp->convolutions2, stamp->image1, stamp->weight, polyValues, 516 551 footprint); 517 status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1, 518 stamp->image2, stamp->weight, polyValues, footprint); 552 status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, 553 stamp->image2, stamp->weight, polyValues, footprint, true); 554 status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, 555 stamp->image2, stamp->weight, polyValues, footprint, false); 519 556 break; 520 557 default: … … 540 577 541 578 matrixName = NULL; 542 psStringAppend(&matrixName, "vector _%d.fits", i);543 psImage *dummy = psImageAlloc(stamp->vector ->n, 1, PS_TYPE_F64);544 memcpy(dummy->data.F64[0], stamp->vector ->data.F64,545 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector ->n);579 psStringAppend(&matrixName, "vector1_%d.fits", i); 580 psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64); 581 memcpy(dummy->data.F64[0], stamp->vector1->data.F64, 582 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n); 546 583 matrixFile = psFitsOpen(matrixName, "w"); 547 584 psFree(matrixName); … … 550 587 psFitsClose(matrixFile); 551 588 589 if (stamp->vector2) { 590 matrixName = NULL; 591 psStringAppend(&matrixName, "vector2_%d.fits", i); 592 dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64); 593 memcpy(dummy->data.F64[0], stamp->vector2->data.F64, 594 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n); 595 matrixFile = psFitsOpen(matrixName, "w"); 596 psFree(matrixName); 597 psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); 598 psFree(dummy); 599 psFitsClose(matrixFile); 600 } 552 601 553 602 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { … … 575 624 } 576 625 577 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps, 578 float tolerance) 626 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) 579 627 { 580 628 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); 581 629 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 582 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {583 PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, false);584 }585 630 586 631 // Check inputs … … 594 639 } 595 640 641 PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false); 642 if (numParams == -1) { 643 numParams = stamp->vector1->n; 644 } 645 PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false); 646 PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false); 647 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); 648 PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); 649 PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); 650 596 651 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 597 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);598 652 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false); 599 653 PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); 600 if (numParams == -1) { 601 numParams = stamp->matrix1->numCols; 654 if (numParams2 == 0) { 602 655 numParams2 = stamp->matrix2->numCols; 603 656 } 604 PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);605 657 PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); 606 658 PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false); 607 PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);608 659 PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false); 609 660 PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false); 610 PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); 611 PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); 612 PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); 613 } else { 614 PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); 615 if (numParams == -1) { 616 numParams = stamp->vector->n; 617 } 618 PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); 619 PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); 620 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); 621 PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); 622 PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); 661 PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false); 662 PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false); 663 PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false); 623 664 } 624 665 } … … 638 679 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 639 680 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); 640 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector );681 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); 641 682 } 642 683 } … … 667 708 psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); 668 709 psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); 669 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); 710 psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); 711 psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64); 670 712 psImageInit(sumMatrix1, 0.0); 671 713 psImageInit(sumMatrix2, 0.0); 672 714 psImageInit(sumMatrixX, 0.0); 673 psVectorInit(sumVector, 0.0); 715 psVectorInit(sumVector1, 0.0); 716 psVectorInit(sumVector2, 0.0); 674 717 675 718 for (int i = 0; i < stamps->num; i++) { … … 679 722 (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2); 680 723 (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); 681 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 682 } 683 } 684 685 psVector *permutation1 = NULL, *permutation2 = NULL; // Permutation vectors for LUD 686 psImage *lu1 = psMatrixLUD(NULL, &permutation1, sumMatrix1); // LUD of matrix 1 687 psImage *lu2 = psMatrixLUD(NULL, &permutation2, sumMatrix2); // LUD of matrix 2 688 psFree(sumMatrix1); 689 psFree(sumMatrix2); 690 if (!lu1 || !lu2) { 691 psError(PS_ERR_UNKNOWN, false, "Unable to generate LU decomposed matrices."); 692 psFree(permutation1); 693 psFree(permutation2); 694 psFree(lu1); 695 psFree(lu2); 696 psFree(sumMatrixX); 697 psFree(sumVector); 698 return NULL; 699 } 700 701 psVector *lastSol2 = psVectorAlloc(numParams2, PS_TYPE_F64); // Value of sol2 in last iteration 702 psVectorInit(lastSol2, 0.0); 703 704 double diff = NAN; // Difference between iterations 705 psVector *vector = NULL; // RHS for equations 706 int iter = 0; // Iteration number 707 do { 708 if (!vector) { 709 // Start with traditional Alard-Lupton solution 710 vector = psVectorCopy(NULL, sumVector, PS_TYPE_F64); 711 } else { 712 // Generate vector for RHS of equation 1 713 vector->n = numParams; 714 psVector *sol2 = kernels->solution2; 715 assert(sol2); 716 for (int i = 0; i < numParams; i++) { 717 double value = 0.0; // Value of vector element 718 for (int j = 0; j < numParams2; j++) { 719 // XXX Not sure I've got this the right way around 720 value += sol2->data.F64[j] * sumMatrixX->data.F64[j][i]; 721 } 722 vector->data.F64[i] = value; 723 } 724 } 725 726 kernels->solution1 = psMatrixLUSolve(kernels->solution1, lu1, vector, permutation1); 727 if (!kernels->solution1) { 728 psError(PS_ERR_UNKNOWN, false, "Unable to solve initial matrix equation."); 729 psFree(permutation1); 730 psFree(permutation2); 731 psFree(lu1); 732 psFree(lu2); 733 psFree(sumMatrixX); 734 psFree(sumVector); 735 psFree(vector); 736 psFree(lastSol2); 737 return NULL; 738 } 739 740 // Generate vector for RHS of equation 2 741 vector->n = numParams2; 742 psVector *sol1 = kernels->solution1; // Solution of first equation 743 for (int i = 0; i < numParams2; i++) { 744 double value = 0.0; // Value of vector element 745 for (int j = 0; j < numParams; j++) { 746 // XXX Not sure I've got this the right way around 747 value += sol1->data.F64[j] * sumMatrixX->data.F64[i][j]; 748 } 749 vector->data.F64[i] = value; 750 } 751 752 kernels->solution2 = psMatrixLUSolve(kernels->solution2, lu2, vector, permutation2); 753 if (!kernels->solution2) { 754 psError(PS_ERR_UNKNOWN, false, "Unable to solve matrix equation 2."); 755 psFree(permutation1); 756 psFree(permutation2); 757 psFree(lu1); 758 psFree(lu2); 759 psFree(sumMatrixX); 760 psFree(sumVector); 761 psFree(vector); 762 psFree(lastSol2); 763 return NULL; 764 } 765 766 // Get difference produced by iteration 767 diff = 0.0; 724 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); 725 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); 726 } 727 } 728 768 729 #if 0 769 psVector *sol2 = kernels->solution2; // Solution of second equation 770 for (int i = 0; i < numParams2; i++) { 771 diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]); 772 } 773 lastSol2 = psVectorCopy(lastSol2, sol2, PS_TYPE_F64); 774 psTrace("psModules.imcombine", 4, "Solution iteration %d: %lf\n", iter, diff); 775 #endif 776 777 psVectorInit(kernels->solution2, 0.0); 778 779 iter++; 780 } while (diff > tolerance); 781 782 psFree(permutation1); 783 psFree(permutation2); 784 psFree(lu1); 785 psFree(lu2); 786 psFree(sumMatrixX); 787 psFree(sumVector); 788 psFree(vector); 789 psFree(lastSol2); 730 // Apply weighting to maximise the difference between solution coefficients for the same functions 731 double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor 732 for (int i = 0; i < kernels->num; i++) { 733 sumMatrix1->data.F64[i][i] -= fudge; 734 sumMatrix2->data.F64[i][i] += fudge; 735 } 736 #endif 737 738 // Pure matrix operations 739 740 // A * a = Ct * b + d 741 // C * a = B * b + e 742 // 743 // a = (Ct * Bi * C - A)i (Ct * Bi * e - d) 744 // b = Bi * (C * a - e) 745 psVector *a = psVectorRecycle(kernels->solution1, numParams, PS_TYPE_F64); 746 psVector *b = psVectorRecycle(kernels->solution2, numParams2, PS_TYPE_F64); 747 #ifdef TESTING 748 psVectorInit(a, NAN); 749 psVectorInit(b, NAN); 750 #endif 751 psImage *A = sumMatrix1; 752 psImage *B = sumMatrix2; 753 psImage *C = sumMatrixX; 754 psVector *d = sumVector1; 755 psVector *e = sumVector2; 756 757 assert(a->n == numParams); 758 assert(b->n == numParams2); 759 assert(A->numRows == numParams && A->numCols == numParams); 760 assert(B->numRows == numParams2 && B->numCols == numParams2); 761 assert(C->numRows == numParams2 && C->numCols == numParams); 762 assert(d->n == numParams); 763 assert(e->n == numParams2); 764 765 psImage *Bi = psMatrixInvert(NULL, B, NULL); 766 assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); 767 psImage *Ct = psMatrixTranspose(NULL, C); 768 assert(Ct->numRows == numParams && Ct->numCols == numParams2); 769 770 psImage *BiC = psMatrixMultiply(NULL, Bi, C); 771 assert(BiC->numRows == numParams2 && BiC->numCols == numParams); 772 psImage *CtBi = psMatrixMultiply(NULL, Ct, Bi); 773 assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); 774 775 psImage *CtBiC = psMatrixMultiply(NULL, Ct, BiC); 776 assert(CtBiC->numRows == numParams && CtBiC->numCols == numParams); 777 778 psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A); 779 assert(F->numRows == numParams && F->numCols == numParams); 780 float det = NAN; 781 psImage *Fi = psMatrixInvert(NULL, F, &det); 782 assert(Fi->numRows == numParams && Fi->numCols == numParams); 783 psTrace("psModules.imcombine", 4, "Determinant of F: %f\n", det); 784 785 psVector *g = psVectorAlloc(numParams, PS_TYPE_F64); 786 #ifdef TESTING 787 psVectorInit(g, NAN); 788 #endif 789 assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); 790 assert(e->n == numParams2); 791 assert(d->n == numParams); 792 for (int i = 0; i < numParams; i++) { 793 double value = 0.0; 794 for (int j = 0; j < numParams2; j++) { 795 value += CtBi->data.F64[i][j] * e->data.F64[j]; 796 } 797 g->data.F64[i] = value - d->data.F64[i]; 798 } 799 800 assert(Fi->numRows == numParams && Fi->numCols == numParams); 801 assert(g->n == numParams); 802 for (int i = 0; i < numParams; i++) { 803 double value = 0.0; 804 for (int j = 0; j < numParams; j++) { 805 value += Fi->data.F64[i][j] * g->data.F64[j]; 806 } 807 a->data.F64[i] = value; 808 } 809 810 psVector *h = psVectorAlloc(numParams2, PS_TYPE_F64); 811 #ifdef TESTING 812 psVectorInit(h, NAN); 813 #endif 814 assert(C->numRows == numParams2 && C->numCols == numParams); 815 assert(a->n == numParams); 816 assert(e->n == numParams2); 817 for (int i = 0; i < numParams2; i++) { 818 double value = 0.0; 819 for (int j = 0; j < numParams; j++) { 820 value += C->data.F64[i][j] * a->data.F64[j]; 821 } 822 h->data.F64[i] = value - e->data.F64[i]; 823 } 824 825 assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); 826 assert(h->n == numParams2); 827 for (int i = 0; i < numParams2; i++) { 828 double value = 0.0; 829 for (int j = 0; j < numParams2; j++) { 830 value += Bi->data.F64[i][j] * h->data.F64[j]; 831 } 832 b->data.F64[i] = value; 833 } 834 835 836 #if 0 837 for (int i = 0; i < numParams; i++) { 838 double aVal1 = 0.0, bVal1 = 0.0; 839 for (int j = 0; j < numParams2; j++) { 840 aVal1 += A->data.F64[i][j] * a->data.F64[j]; 841 bVal1 += Ct->data.F64[i][j] * b->data.F64[j]; 842 } 843 bVal1 += d->data.F64[i]; 844 for (int j = numParams2; j < numParams; j++) { 845 aVal1 += A->data.F64[i][j] * a->data.F64[j]; 846 } 847 printf("%d: %lf\n", i, aVal1 - bVal1); 848 } 849 850 for (int i = 0; i < numParams2; i++) { 851 double aVal2 = 0.0, bVal2 = 0.0; 852 for (int j = 0; j < numParams2; j++) { 853 aVal2 += C->data.F64[i][j] * a->data.F64[j]; 854 bVal2 += B->data.F64[i][j] * b->data.F64[j]; 855 } 856 bVal2 += e->data.F64[i]; 857 for (int j = numParams2; j < numParams; j++) { 858 aVal2 += C->data.F64[i][j] * a->data.F64[j]; 859 } 860 printf("%d: %lf\n", i, aVal2 - bVal2); 861 } 862 #endif 863 864 { 865 psFits *fits = psFitsOpen("sumMatrix1.fits", "w"); 866 psFitsWriteImage(fits, NULL, sumMatrix1, 0, NULL); 867 psFitsClose(fits); 868 } 869 { 870 psFits *fits = psFitsOpen("sumMatrix2.fits", "w"); 871 psFitsWriteImage(fits, NULL, sumMatrix2, 0, NULL); 872 psFitsClose(fits); 873 } 874 { 875 psFits *fits = psFitsOpen("sumMatrixX.fits", "w"); 876 psFitsWriteImage(fits, NULL, sumMatrixX, 0, NULL); 877 psFitsClose(fits); 878 } 879 { 880 psFits *fits = psFitsOpen("sumFinverse.fits", "w"); 881 psFitsWriteImage(fits, NULL, Fi, 0, NULL); 882 psFitsClose(fits); 883 } 884 885 886 kernels->solution1 = a; 887 kernels->solution2 = b; 888 889 // XXXXX Free temporary matrices and vectors 890 790 891 } 791 892 … … 817 918 int numKernels = kernels->num; // Number of kernels 818 919 819 psVector *coeff1 = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for first image820 psVector *coeff2 = kernels->mode == PM_SUBTRACTION_MODE_DUAL ?821 psVectorAlloc(numKernels, PS_TYPE_F64) : NULL; // Coefficients for second image822 920 psImage *polyValues = NULL; // Polynomial values 823 921 psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image … … 832 930 // Calculate coefficients of the kernel basis functions 833 931 polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); 834 for (int j = 0; j < numKernels; j++) {835 coeff1->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false);836 }837 932 double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation 838 933 double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background … … 863 958 for (int j = 0; j < numKernels; j++) { 864 959 psKernel *convolution = convolutions->data[j]; // Convolution 865 double coefficient = coeff1->data.F64[j]; // Coefficient 960 double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, 961 false); // Coefficient 866 962 for (int y = - footprint; y <= footprint; y++) { 867 963 for (int x = - footprint; x <= footprint; x++) { … … 880 976 psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image 881 977 psKernel *image1 = stamp->image1; // The first image 882 883 for (int j = 0; j < numKernels; j++) { 884 coeff2->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); 885 } 978 psKernel *image2 = stamp->image2; // The second image 886 979 887 980 for (int j = 0; j < numKernels; j++) { 888 981 psKernel *conv1 = convolutions1->data[j]; // Convolution of first image 889 982 psKernel *conv2 = convolutions2->data[j]; // Convolution of second image 890 891 double coefficient1 = coeff1->data.F64[j]; // Coefficient for convolution 1 892 double coefficient2 = coeff2->data.F64[j]; // Coefficient for convolution 2 983 double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 984 double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 893 985 894 986 for (int y = - footprint; y <= footprint; y++) { 895 987 for (int x = - footprint; x <= footprint; x++) { 896 residual->kernel[y][x] += conv1->kernel[y][x] * coefficient1 + 897 conv2->kernel[y][x] * coefficient2; 988 residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1; 898 989 } 899 990 } … … 901 992 for (int y = - footprint; y <= footprint; y++) { 902 993 for (int x = - footprint; x <= footprint; x++) { 903 residual->kernel[y][x] += image 1->kernel[y][x] * norm - background;994 residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm; 904 995 } 905 996 } … … 909 1000 for (int y = - footprint; y <= footprint; y++) { 910 1001 for (int x = - footprint; x <= footprint; x++) { 911 deviation += PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; 1002 double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; 1003 deviation += dev; 1004 #ifdef TESTING 1005 residual->kernel[y][x] = dev; 1006 #endif 912 1007 } 913 1008 } … … 932 1027 psFitsClose(fits); 933 1028 } 1029 { 1030 psString filename = NULL; 1031 psStringAppend(&filename, "stamp_image_%03d.fits", i); 1032 psFits *fits = psFitsOpen(filename, "w"); 1033 psFree(filename); 1034 psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); 1035 psFitsClose(fits); 1036 } 1037 { 1038 psString filename = NULL; 1039 psStringAppend(&filename, "stamp_weight_%03d.fits", i); 1040 psFits *fits = psFitsOpen(filename, "w"); 1041 psFree(filename); 1042 psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); 1043 psFitsClose(fits); 1044 } 934 1045 #endif 935 1046 … … 937 1048 psFree(residual); 938 1049 psFree(polyValues); 939 psFree(coeff1);940 1050 941 1051 return deviations;
Note:
See TracChangeset
for help on using the changeset viewer.
