Changeset 26893 for trunk/psModules/src/imcombine/pmSubtractionEquation.c
- Timestamp:
- Feb 10, 2010, 7:34:39 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r26035 r26893 17 17 //#define TESTING // TESTING output for debugging; may not work with threads! 18 18 19 #define USE_WEIGHT // Include weight (1/variance) in equation? 19 //#define USE_WEIGHT // Include weight (1/variance) in equation? 20 //#define USE_WINDOW // Include weight (1/variance) in equation? 20 21 21 22 … … 27 28 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated 28 29 psVector *vector, // Least-squares vector, updated 30 double *norm, // Normalisation, updated 29 31 const psKernel *input, // Input image (target) 30 32 const psKernel *reference, // Reference image (convolution source) 31 33 const psKernel *weight, // Weight image 34 const psKernel *window, // Window image 32 35 const psArray *convolutions, // Convolutions for each kernel 33 36 const pmSubtractionKernels *kernels, // Kernels 34 37 const psImage *polyValues, // Spatial polynomial values 35 int footprint // (Half-)Size of stamp 38 int footprint, // (Half-)Size of stamp 39 int normWindow, // Window (half-)size for normalisation measurement 40 const pmSubtractionEquationCalculationMode mode 36 41 ) 37 42 { … … 51 56 52 57 // Evaluate polynomial-polynomial terms 58 // XXX we can skip this if we are not calculating kernel coeffs 53 59 for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { 54 60 for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { … … 64 70 } 65 71 72 // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we 73 // choose to calculate 74 psImageInit(matrix, 0.0); 75 psVectorInit(vector, 1.0); 76 for (int i = 0; i < matrix->numCols; i++) { 77 matrix->data.F64[i][i] = 1.0; 78 } 79 80 // the order of the elements in the matrix and vector is: 81 // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] 82 // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] 83 // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] 84 // normalization 85 // bg 0, bg 1, bg 2 (only 0 is currently used?) 66 86 67 87 for (int i = 0; i < numKernels; i++) { … … 74 94 for (int x = - footprint; x <= footprint; x++) { 75 95 double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; 76 #ifdef USE_WEIGHT 77 cc *= weight->kernel[y][x]; 78 #endif 96 if (weight) { 97 cc *= weight->kernel[y][x]; 98 } 99 if (window) { 100 cc *= window->kernel[y][x]; 101 } 79 102 sumCC += cc; 80 103 } 81 104 } 82 105 83 // Spatial variation 84 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 85 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 86 double value = sumCC * poly2[iTerm][jTerm]; 87 matrix->data.F64[iIndex][jIndex] = value; 88 matrix->data.F64[jIndex][iIndex] = value; 106 // Spatial variation of kernel coeffs 107 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 108 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 109 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 110 double value = sumCC * poly2[iTerm][jTerm]; 111 matrix->data.F64[iIndex][jIndex] = value; 112 matrix->data.F64[jIndex][iIndex] = value; 113 } 89 114 } 90 115 } … … 102 127 double rc = ref * conv; 103 128 double c = conv; 104 #ifdef USE_WEIGHT 105 float wtVal = weight->kernel[y][x]; 106 ic *= wtVal; 107 rc *= wtVal; 108 c *= wtVal; 109 #endif 129 if (weight) { 130 float wtVal = weight->kernel[y][x]; 131 ic *= wtVal; 132 rc *= wtVal; 133 c *= wtVal; 134 } 135 if (window) { 136 float winVal = window->kernel[y][x]; 137 ic *= winVal; 138 rc *= winVal; 139 c *= winVal; 140 } 110 141 sumIC += ic; 111 142 sumRC += rc; … … 117 148 double normTerm = sumRC * poly[iTerm]; 118 149 double bgTerm = sumC * poly[iTerm]; 119 matrix->data.F64[iIndex][normIndex] = normTerm; 120 matrix->data.F64[normIndex][iIndex] = normTerm; 121 matrix->data.F64[iIndex][bgIndex] = bgTerm; 122 matrix->data.F64[bgIndex][iIndex] = bgTerm; 123 vector->data.F64[iIndex] = sumIC * poly[iTerm]; 150 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 151 matrix->data.F64[iIndex][normIndex] = normTerm; 152 matrix->data.F64[normIndex][iIndex] = normTerm; 153 } 154 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 155 matrix->data.F64[iIndex][bgIndex] = bgTerm; 156 matrix->data.F64[bgIndex][iIndex] = bgTerm; 157 } 158 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 159 vector->data.F64[iIndex] = sumIC * poly[iTerm]; 160 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 161 // subtract norm * sumRC * poly[iTerm] 162 psAssert (kernels->solution1, "programming error: define solution first!"); 163 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 164 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 165 vector->data.F64[iIndex] -= norm * normTerm; 166 } 167 } 124 168 } 125 169 } … … 130 174 double sumR = 0.0; // Sum of the reference 131 175 double sumI = 0.0; // Sum of the input 176 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 132 177 for (int y = - footprint; y <= footprint; y++) { 133 178 for (int x = - footprint; x <= footprint; x++) { … … 137 182 double rr = PS_SQR(ref); 138 183 double one = 1.0; 139 #ifdef USE_WEIGHT 140 float wtVal = weight->kernel[y][x]; 141 rr *= wtVal; 142 ir *= wtVal; 143 in *= wtVal; 144 ref *= wtVal; 145 one *= wtVal; 146 #endif 184 185 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { 186 normI1 += ref; 187 normI2 += in; 188 } 189 190 if (weight) { 191 float wtVal = weight->kernel[y][x]; 192 rr *= wtVal; 193 ir *= wtVal; 194 in *= wtVal; 195 ref *= wtVal; 196 one *= wtVal; 197 } 198 if (window) { 199 float winVal = window->kernel[y][x]; 200 rr *= winVal; 201 ir *= winVal; 202 in *= winVal; 203 ref *= winVal; 204 one *= winVal; 205 } 147 206 sumRR += rr; 148 207 sumIR += ir; … … 152 211 } 153 212 } 154 matrix->data.F64[normIndex][normIndex] = sumRR; 155 matrix->data.F64[bgIndex][bgIndex] = sum1; 156 matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR; 157 vector->data.F64[normIndex] = sumIR; 158 vector->data.F64[bgIndex] = sumI; 213 214 *norm = normI2 / normI1; 215 216 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 217 matrix->data.F64[normIndex][normIndex] = sumRR; 218 vector->data.F64[normIndex] = sumIR; 219 // subtract sum over kernels * kernel solution 220 } 221 if (mode & PM_SUBTRACTION_EQUATION_BG) { 222 matrix->data.F64[bgIndex][bgIndex] = sum1; 223 vector->data.F64[bgIndex] = sumI; 224 } 225 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 226 matrix->data.F64[normIndex][bgIndex] = sumR; 227 matrix->data.F64[bgIndex][normIndex] = sumR; 228 } 229 230 // check for any NAN values in the result, skip if found: 231 for (int iy = 0; iy < matrix->numRows; iy++) { 232 for (int ix = 0; ix < matrix->numCols; ix++) { 233 if (!isfinite(matrix->data.F64[iy][ix])) { 234 fprintf (stderr, "WARNING: NAN in matrix\n"); 235 return false; 236 } 237 } 238 } 239 for (int ix = 0; ix < vector->n; ix++) { 240 if (!isfinite(vector->data.F64[ix])) { 241 fprintf (stderr, "WARNING: NAN in vector\n"); 242 return false; 243 } 244 } 159 245 160 246 return true; 161 247 } 162 248 249 163 250 // Calculate the least-squares matrix and vector for dual convolution 164 static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated 165 psVector *vector1, // Least-squares vector, updated 166 psImage *matrix2, // Least-squares matrix, updated 167 psVector *vector2, // Least-squares vector, updated 168 psImage *matrixX, // Cross-matrix 251 static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated 252 psVector *vector, // Least-squares vector, updated 253 double *norm, // Normalisation, updated 169 254 const psKernel *image1, // Image 1 170 255 const psKernel *image2, // Image 2 171 256 const psKernel *weight, // Weight image 257 const psKernel *window, // Window image 172 258 const psArray *convolutions1, // Convolutions of image 1 for each kernel 173 259 const psArray *convolutions2, // Convolutions of image 2 for each kernel 174 260 const pmSubtractionKernels *kernels, // Kernels 175 261 const psImage *polyValues, // Spatial polynomial values 176 int footprint // (Half-)Size of stamp 262 int footprint, // (Half-)Size of stamp 263 int normWindow, // Window (half-)size for normalisation measurement 264 const pmSubtractionEquationCalculationMode mode 177 265 ) 178 266 { 179 // A_ij = A_i A_j180 // B_ij = B_i B_j181 // C_ij = A_i B_j182 // d_i = A_i I_2183 // e_i = B_i I_2184 185 // A_i = I_1 * k_i186 // B_i = I_2 * k_i187 188 // Background: A_i = 1.0189 // Normalisation: A_i = I_1190 191 267 int numKernels = kernels->num; // Number of kernels 192 268 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation … … 196 272 double poly[numPoly]; // Polynomial terms 197 273 double poly2[numPoly][numPoly]; // Polynomial-polynomial values 274 275 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 276 int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters 277 int numParams2 = numKernels * numPoly; // Number of additional parameters for dual 278 int numDual = numParams + numParams2; // Total number of parameters for dual 279 280 psAssert(matrix && 281 matrix->type.type == PS_TYPE_F64 && 282 matrix->numCols == numDual && 283 matrix->numRows == numDual, 284 "Least-squares matrix is bad."); 285 psAssert(vector && 286 vector->type.type == PS_TYPE_F64 && 287 vector->n == numDual, 288 "Least-squares vector is bad."); 198 289 199 290 // Evaluate polynomial-polynomial terms … … 212 303 213 304 305 // initialize the matrix and vector for NOP on all coeffs. we only fill in the coeffs we 306 // choose to calculate 307 psImageInit(matrix, 0.0); 308 psVectorInit(vector, 1.0); 309 for (int i = 0; i < matrix->numCols; i++) { 310 matrix->data.F64[i][i] = 1.0; 311 } 312 214 313 for (int i = 0; i < numKernels; i++) { 215 314 psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i … … 219 318 psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j 220 319 221 double sumAA = 0.0; // Sum of convolution products for matrix A222 double sumBB = 0.0; // Sum of convolution products for matrix B223 double sumAB = 0.0; // Sum of convolution products for matrix C320 double sumAA = 0.0; // Sum of convolution products between image 1 321 double sumBB = 0.0; // Sum of convolution products between image 2 322 double sumAB = 0.0; // Sum of convolution products across images 1 and 2 224 323 for (int y = - footprint; y <= footprint; y++) { 225 324 for (int x = - footprint; x <= footprint; x++) { … … 227 326 double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; 228 327 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 229 #ifdef USE_WEIGHT 230 float wtVal = weight->kernel[y][x]; 231 aa *= wtVal; 232 bb *= wtVal; 233 ab *= wtVal; 234 #endif 328 if (weight) { 329 float wtVal = weight->kernel[y][x]; 330 aa *= wtVal; 331 bb *= wtVal; 332 ab *= wtVal; 333 } 334 if (window) { 335 float wtVal = window->kernel[y][x]; 336 aa *= wtVal; 337 bb *= wtVal; 338 ab *= wtVal; 339 } 235 340 sumAA += aa; 236 341 sumBB += bb; … … 239 344 } 240 345 241 // Spatial variation 242 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 243 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 244 double aa = sumAA * poly2[iTerm][jTerm]; 245 double bb = sumBB * poly2[iTerm][jTerm]; 246 double ab = sumAB * poly2[iTerm][jTerm]; 247 matrix1->data.F64[iIndex][jIndex] = aa; 248 matrix1->data.F64[jIndex][iIndex] = aa; 249 matrix2->data.F64[iIndex][jIndex] = bb; 250 matrix2->data.F64[jIndex][iIndex] = bb; 251 matrixX->data.F64[iIndex][jIndex] = ab; 346 // Spatial variation of kernel coeffs 347 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 348 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 349 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 350 double aa = sumAA * poly2[iTerm][jTerm]; 351 double bb = sumBB * poly2[iTerm][jTerm]; 352 double ab = sumAB * poly2[iTerm][jTerm]; 353 354 matrix->data.F64[iIndex][jIndex] = aa; 355 matrix->data.F64[jIndex][iIndex] = aa; 356 357 matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; 358 matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; 359 360 matrix->data.F64[iIndex][jIndex + numParams] = ab; 361 matrix->data.F64[jIndex + numParams][iIndex] = ab; 362 } 252 363 } 253 364 } … … 259 370 for (int x = - footprint; x <= footprint; x++) { 260 371 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 261 #ifdef USE_WEIGHT 262 ab *= weight->kernel[y][x]; 263 #endif 372 if (weight) { 373 ab *= weight->kernel[y][x]; 374 } 375 if (window) { 376 ab *= window->kernel[y][x]; 377 } 264 378 sumAB += ab; 265 379 } 266 380 } 267 381 268 // Spatial variation 269 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 270 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 271 double ab = sumAB * poly2[iTerm][jTerm]; 272 matrixX->data.F64[iIndex][jIndex] = ab; 273 } 274 } 275 } 276 277 double sumAI2 = 0.0; // Sum of A.I_2 products (for vector 1) 278 double sumBI2 = 0.0; // Sum of B.I_2 products (for vector 2) 279 double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix 1, normalisation) 280 double sumA = 0.0; // Sum of A (for matrix 1, background) 281 double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix X, normalisation) 282 double sumB = 0.0; // Sum of B products (for matrix X, background) 283 double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) 284 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) 382 // Spatial variation of kernel coeffs 383 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 384 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 385 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 386 double ab = sumAB * poly2[iTerm][jTerm]; 387 matrix->data.F64[iIndex][jIndex + numParams] = ab; 388 matrix->data.F64[jIndex + numParams][iIndex] = ab; 389 } 390 } 391 } 392 } 393 394 double sumAI2 = 0.0; // Sum of A.I_2 products (for vector) 395 double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) 396 double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) 397 double sumA = 0.0; // Sum of A (for matrix, background) 398 double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) 399 double sumB = 0.0; // Sum of B products (for matrix, background) 400 double sumI2 = 0.0; // Sum of I_2 (for vector, background) 285 401 for (int y = - footprint; y <= footprint; y++) { 286 402 for (int x = - footprint; x <= footprint; x++) { 287 floata = iConv1->kernel[y][x];288 floatb = iConv2->kernel[y][x];403 double a = iConv1->kernel[y][x]; 404 double b = iConv2->kernel[y][x]; 289 405 float i1 = image1->kernel[y][x]; 290 406 float i2 = image2->kernel[y][x]; … … 294 410 double ai1 = a * i1; 295 411 double bi1 = b * i1; 296 double i1i2 = i1 * i2; 297 298 #ifdef USE_WEIGHT 299 float wtVal = weight->kernel[y][x]; 300 ai2 *= wtVal; 301 bi2 *= wtVal; 302 ai1 *= wtVal; 303 bi1 *= wtVal; 304 i1i2 *= wtVal; 305 a *= wtVal; 306 b *= wtVal; 307 i2 *= wtVal; 308 #endif 309 412 413 if (weight) { 414 float wtVal = weight->kernel[y][x]; 415 ai2 *= wtVal; 416 bi2 *= wtVal; 417 ai1 *= wtVal; 418 bi1 *= wtVal; 419 a *= wtVal; 420 b *= wtVal; 421 i2 *= wtVal; 422 } 423 if (window) { 424 float wtVal = window->kernel[y][x]; 425 ai2 *= wtVal; 426 bi2 *= wtVal; 427 ai1 *= wtVal; 428 bi1 *= wtVal; 429 a *= wtVal; 430 b *= wtVal; 431 i2 *= wtVal; 432 } 310 433 sumAI2 += ai2; 311 434 sumBI2 += bi2; … … 315 438 sumB += b; 316 439 sumI2 += i2; 317 sumI1I2 += i1i2;318 440 } 319 441 } … … 323 445 double bi2 = sumBI2 * poly[iTerm]; 324 446 double ai1 = sumAI1 * poly[iTerm]; 325 double a = sumA * poly[iTerm];447 double a = sumA * poly[iTerm]; 326 448 double bi1 = sumBI1 * poly[iTerm]; 327 double b = sumB * poly[iTerm]; 328 329 matrix1->data.F64[iIndex][normIndex] = ai1; 330 matrix1->data.F64[normIndex][iIndex] = ai1; 331 matrix1->data.F64[iIndex][bgIndex] = a; 332 matrix1->data.F64[bgIndex][iIndex] = a; 333 vector1->data.F64[iIndex] = ai2; 334 vector2->data.F64[iIndex] = bi2; 335 matrixX->data.F64[iIndex][normIndex] = bi1; 336 matrixX->data.F64[iIndex][bgIndex] = b; 337 } 338 } 339 340 double sumI1 = 0.0; // Sum of I_1 (for matrix 1, background-normalisation) 341 double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix 1, normalisation-normalisation) 342 double sum1 = 0.0; // Sum of 1 (for matrix 1, background-background) 343 double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) 344 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) 449 double b = sumB * poly[iTerm]; 450 451 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 452 matrix->data.F64[iIndex][normIndex] = ai1; 453 matrix->data.F64[normIndex][iIndex] = ai1; 454 matrix->data.F64[iIndex + numParams][normIndex] = bi1; 455 matrix->data.F64[normIndex][iIndex + numParams] = bi1; 456 } 457 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 458 matrix->data.F64[iIndex][bgIndex] = a; 459 matrix->data.F64[bgIndex][iIndex] = a; 460 matrix->data.F64[iIndex + numParams][bgIndex] = b; 461 matrix->data.F64[bgIndex][iIndex + numParams] = b; 462 } 463 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 464 vector->data.F64[iIndex] = ai2; 465 vector->data.F64[iIndex + numParams] = bi2; 466 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 467 // subtract norm * sumRC * poly[iTerm] 468 psAssert (kernels->solution1, "programming error: define solution first!"); 469 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 470 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 471 vector->data.F64[iIndex] -= norm * ai1; 472 vector->data.F64[iIndex + numParams] -= norm * bi1; 473 } 474 } 475 } 476 } 477 478 double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) 479 double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) 480 double sum1 = 0.0; // Sum of 1 (for matrix, background-background) 481 double sumI2 = 0.0; // Sum of I_2 (for vector, background) 482 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) 483 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 345 484 for (int y = - footprint; y <= footprint; y++) { 346 485 for (int x = - footprint; x <= footprint; x++) { 347 floati1 = image1->kernel[y][x];348 floati2 = image2->kernel[y][x];486 double i1 = image1->kernel[y][x]; 487 double i2 = image2->kernel[y][x]; 349 488 350 489 double i1i1 = i1 * i1; … … 352 491 double i1i2 = i1 * i2; 353 492 354 #ifdef USE_WEIGHT 355 float wtVal = weight->kernel[y][x]; 356 i1 *= wtVal; 357 i1i1 *= wtVal; 358 one *= wtVal; 359 i2 *= wtVal; 360 i1i2 *= wtVal; 361 #endif 362 493 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) { 494 normI1 += i1; 495 normI2 += i2; 496 } 497 498 if (weight) { 499 float wtVal = weight->kernel[y][x]; 500 i1 *= wtVal; 501 i1i1 *= wtVal; 502 one *= wtVal; 503 i2 *= wtVal; 504 i1i2 *= wtVal; 505 } 506 if (window) { 507 float wtVal = window->kernel[y][x]; 508 i1 *= wtVal; 509 i1i1 *= wtVal; 510 one *= wtVal; 511 i2 *= wtVal; 512 i1i2 *= wtVal; 513 } 363 514 sumI1 += i1; 364 515 sumI1I1 += i1i1; … … 368 519 } 369 520 } 370 matrix1->data.F64[bgIndex][normIndex] = sumI1; 371 matrix1->data.F64[normIndex][bgIndex] = sumI1; 372 matrix1->data.F64[normIndex][normIndex] = sumI1I1; 373 matrix1->data.F64[bgIndex][bgIndex] = sum1; 374 vector1->data.F64[bgIndex] = sumI2; 375 vector1->data.F64[normIndex] = sumI1I2; 521 522 *norm = normI2 / normI1; 523 524 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 525 matrix->data.F64[normIndex][normIndex] = sumI1I1; 526 vector->data.F64[normIndex] = sumI1I2; 527 } 528 if (mode & PM_SUBTRACTION_EQUATION_BG) { 529 matrix->data.F64[bgIndex][bgIndex] = sum1; 530 vector->data.F64[bgIndex] = sumI2; 531 } 532 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 533 matrix->data.F64[bgIndex][normIndex] = sumI1; 534 matrix->data.F64[normIndex][bgIndex] = sumI1; 535 } 536 537 // check for any NAN values in the result, skip if found: 538 for (int iy = 0; iy < matrix->numRows; iy++) { 539 for (int ix = 0; ix < matrix->numCols; ix++) { 540 if (!isfinite(matrix->data.F64[iy][ix])) { 541 fprintf (stderr, "WARNING: NAN in matrix\n"); 542 return false; 543 } 544 } 545 } 546 for (int ix = 0; ix < vector->n; ix++) { 547 if (!isfinite(vector->data.F64[ix])) { 548 fprintf (stderr, "WARNING: NAN in vector\n"); 549 return false; 550 } 551 } 552 376 553 377 554 return true; 378 555 } 379 556 380 // Merge dual matrices and vectors into single matrix equation 381 // Have: Aa = Ct.b + d 382 // Have: Ca = Bb + e 383 // Set: F = ( A -Ct ; C -B ) 384 // Set: g = ( a ; b ) 385 // Set: h = ( d ; e ) 386 // So that we combine the above two equations: Fg = h 387 static bool calculateEquationDual(psImage **outMatrix, 388 psVector **outVector, 389 const psImage *sumMatrix1, 390 const psImage *sumMatrix2, 391 const psImage *sumMatrixX, 392 const psVector *sumVector1, 393 const psVector *sumVector2 394 ) 395 { 396 psAssert(sumMatrix1 && sumMatrix2 && sumMatrixX, "Require input matrices"); 397 psAssert(sumVector1 && sumVector2, "Require input vectors"); 398 int num1 = sumVector1->n; // Number of parameters in first set 399 int num2 = sumVector2->n; // Number of parameters in second set 400 int num = num1 + num2; // Number of parameters in new set 401 402 psAssert(sumMatrix1->type.type == PS_TYPE_F64 && 403 sumMatrix2->type.type == PS_TYPE_F64 && 404 sumMatrixX->type.type == PS_TYPE_F64 && 405 sumVector1->type.type == PS_TYPE_F64 && 406 sumVector2->type.type == PS_TYPE_F64, 407 "Require input type is F64"); 408 409 psAssert(outMatrix, "Require output matrix"); 410 psAssert(outVector, "Require output vector"); 411 if (!*outMatrix) { 412 *outMatrix = psImageAlloc(num, num, PS_TYPE_F64); 413 } 414 if (!*outVector) { 415 *outVector = psVectorAlloc(num, PS_TYPE_F64); 416 } 417 psImage *matrix = *outMatrix; 418 psVector *vector = *outVector; 419 420 psAssert(sumMatrix1->numCols == num1 && sumMatrix1->numRows == num1, "Require size NxN"); 421 psAssert(sumMatrix2->numCols == num2 && sumMatrix2->numRows == num2, "Require size MxM"); 422 psAssert(sumMatrixX->numCols == num1 && sumMatrixX->numRows == num2, "Require size MxN"); 423 424 memcpy(vector->data.F64, sumVector1->data.F64, num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 425 memcpy(&vector->data.F64[num1], sumVector2->data.F64, num2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 426 427 for (int i = 0; i < num1; i++) { 428 memcpy(matrix->data.F64[i], sumMatrix1->data.F64[i], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 429 for (int j = 0, k = num1; j < num2; j++, k++) { 430 matrix->data.F64[i][k] = - sumMatrixX->data.F64[j][i]; 431 } 432 } 433 for (int i1 = 0, i2 = num1; i1 < num2; i1++, i2++) { 434 memcpy(matrix->data.F64[i2], sumMatrixX->data.F64[i1], num1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 435 for (int j = 0, k = num1; j < num2; j++, k++) { 436 matrix->data.F64[i2][k] = - sumMatrix2->data.F64[i1][j]; 437 } 438 } 439 440 return true; 441 } 442 443 557 #if 1 444 558 // Add in penalty term to least-squares vector 445 staticbool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term559 bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term 446 560 psVector *vector, // Vector to which to add in penalty term 447 561 const pmSubtractionKernels *kernels, // Kernel parameters … … 456 570 int spatialOrder = kernels->spatialOrder; // Order of spatial variations 457 571 int numKernels = kernels->num; // Number of kernel components 572 int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations 573 int numParams = numKernels * numSpatial; // Number of kernel parameters 574 575 // order is : 576 // [p_0,x_0,y_0 p_1,x_0,y_0, p_2,x_0,y_0] 577 // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] 578 // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] 579 // [norm] 580 // [bg] 581 // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] 582 // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] 583 // [q_0,x_0,y_1 q_1,x_0,y_1, q_2,x_0,y_1] 584 458 585 for (int i = 0; i < numKernels; i++) { 459 586 for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { 460 587 for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { 461 588 // Contribution to chi^2: a_i^2 P_i 462 matrix->data.F64[index][index] -= norm * penalties->data.F32[i]; 589 psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty"); 590 matrix->data.F64[index][index] += norm * penalties->data.F32[i]; 591 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 592 matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i]; 593 // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) 594 // penalties scale with second moments 595 // 596 } 463 597 } 464 598 } … … 467 601 return true; 468 602 } 603 # endif 469 604 470 605 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 476 611 // Calculate the value of a polynomial, specified by coefficients and polynomial values 477 612 double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients 478 479 480 481 482 613 const psImage *polyValues, // Polynomial values 614 int order, // Order of polynomials 615 int index, // Index at which to begin 616 int step // Step between subsequent indices 617 ) 483 618 { 484 619 double sum = 0.0; // Value of the polynomial sum … … 495 630 496 631 double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues, 497 632 int index, bool wantDual) 498 633 { 499 634 #if 0 … … 548 683 const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels 549 684 int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index 550 551 return pmSubtractionCalculateEquationStamp(stamps, kernels, index); 685 pmSubtractionEquationCalculationMode mode = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model 686 687 return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode); 552 688 } 553 689 554 690 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, 555 int index )691 int index, const pmSubtractionEquationCalculationMode mode) 556 692 { 557 693 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); … … 566 702 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 567 703 704 // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). 705 // = \sum_i^N_Gaussians [(order + 1) * (order + 2) / 2], eg for 1 Gauss and 1st order, numKernels = 3 706 568 707 // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the 569 708 // number of coefficients for the spatial polynomial, normalisation and a constant background offset. 570 709 int numParams = numKernels * numSpatial + 1 + numBackground; 710 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 711 // An additional image is convolved 712 numParams += numKernels * numSpatial; 713 } 571 714 572 715 pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest 573 716 psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state."); 574 717 575 // Generate convolutions 718 // Generate convolutions: these are generated once and saved 576 719 if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { 577 720 psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", index); … … 603 746 #endif 604 747 748 // XXX visualize the set of convolved stamps 749 605 750 psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, 606 751 stamp->xNorm, stamp->yNorm); // Polynomial terms 607 752 608 bool new = stamp->vector 1? false : true; // Is this a new run?753 bool new = stamp->vector ? false : true; // Is this a new run? 609 754 if (new) { 610 stamp->matrix 1= psImageAlloc(numParams, numParams, PS_TYPE_F64);611 stamp->vector 1= psVectorAlloc(numParams, PS_TYPE_F64);755 stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); 756 stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); 612 757 } 613 758 #ifdef TESTING 614 psImageInit(stamp->matrix 1, NAN);615 psVectorInit(stamp->vector 1, NAN);759 psImageInit(stamp->matrix, NAN); 760 psVectorInit(stamp->vector, NAN); 616 761 #endif 617 762 618 763 bool status; // Status of least-squares matrix/vector calculation 764 765 psKernel *weight = NULL; 766 psKernel *window = NULL; 767 768 #ifdef USE_WEIGHT 769 weight = stamp->weight; 770 #endif 771 #ifdef USE_WINDOW 772 window = stamps->window; 773 #endif 774 619 775 switch (kernels->mode) { 620 776 case PM_SUBTRACTION_MODE_1: 621 status = calculateMatrixVector(stamp->matrix 1, stamp->vector1, stamp->image2, stamp->image1,622 stamp->weight, stamp->convolutions1, kernels, polyValues,623 footprint);777 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, 778 weight, window, stamp->convolutions1, kernels, 779 polyValues, footprint, stamps->normWindow, mode); 624 780 break; 625 781 case PM_SUBTRACTION_MODE_2: 626 status = calculateMatrixVector(stamp->matrix 1, stamp->vector1, stamp->image1, stamp->image2,627 stamp->weight, stamp->convolutions2, kernels, polyValues,628 footprint);782 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, 783 weight, window, stamp->convolutions2, kernels, 784 polyValues, footprint, stamps->normWindow, mode); 629 785 break; 630 786 case PM_SUBTRACTION_MODE_DUAL: 631 if (new) { 632 stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); 633 stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64); 634 stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); 635 } 636 #ifdef TESTING 637 psImageInit(stamp->matrix2, NAN); 638 psImageInit(stamp->matrixX, NAN); 639 psVectorInit(stamp->vector2, NAN); 640 #endif 641 status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, 642 stamp->matrixX, stamp->image1, stamp->image2, stamp->weight, 643 stamp->convolutions1, stamp->convolutions2, kernels, polyValues, 644 footprint); 787 status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, 788 stamp->image1, stamp->image2, 789 weight, window, stamp->convolutions1, stamp->convolutions2, 790 kernels, polyValues, footprint, stamps->normWindow, mode); 645 791 break; 646 792 default: … … 651 797 stamp->status = PM_SUBTRACTION_STAMP_REJECTED; 652 798 psWarning("Rejecting stamp %d (%d,%d) because of bad equation", 653 index, (int)(stamp->x + 0.5), (int)(stamp->y +0.5));799 index, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); 654 800 } else { 655 801 stamp->status = PM_SUBTRACTION_STAMP_USED; … … 659 805 { 660 806 psString matrixName = NULL; 661 psStringAppend(&matrixName, "matrix 1_%d.fits", index);807 psStringAppend(&matrixName, "matrix_%d.fits", index); 662 808 psFits *matrixFile = psFitsOpen(matrixName, "w"); 663 809 psFree(matrixName); 664 psFitsWriteImage(matrixFile, NULL, stamp->matrix 1, 0, NULL);810 psFitsWriteImage(matrixFile, NULL, stamp->matrix, 0, NULL); 665 811 psFitsClose(matrixFile); 666 812 667 813 matrixName = NULL; 668 psStringAppend(&matrixName, "vector 1_%d.fits", index);669 psImage *dummy = psImageAlloc(stamp->vector 1->n, 1, PS_TYPE_F64);670 memcpy(dummy->data.F64[0], stamp->vector 1->data.F64,671 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector 1->n);814 psStringAppend(&matrixName, "vector_%d.fits", index); 815 psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64); 816 memcpy(dummy->data.F64[0], stamp->vector->data.F64, 817 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n); 672 818 matrixFile = psFitsOpen(matrixName, "w"); 673 819 psFree(matrixName); … … 675 821 psFree(dummy); 676 822 psFitsClose(matrixFile); 677 678 if (stamp->vector2) {679 matrixName = NULL;680 psStringAppend(&matrixName, "vector2_%d.fits", index);681 dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);682 memcpy(dummy->data.F64[0], stamp->vector2->data.F64,683 PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);684 matrixFile = psFitsOpen(matrixName, "w");685 psFree(matrixName);686 psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);687 psFree(dummy);688 psFitsClose(matrixFile);689 }690 691 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {692 matrixName = NULL;693 psStringAppend(&matrixName, "matrix2_%d.fits", index);694 matrixFile = psFitsOpen(matrixName, "w");695 psFree(matrixName);696 psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL);697 psFitsClose(matrixFile);698 699 matrixName = NULL;700 psStringAppend(&matrixName, "matrixX_%d.fits", index);701 matrixFile = psFitsOpen(matrixName, "w");702 psFree(matrixName);703 psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL);704 psFitsClose(matrixFile);705 }706 823 } 707 824 #endif … … 712 829 } 713 830 714 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels) 831 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, 832 const pmSubtractionEquationCalculationMode mode) 715 833 { 716 834 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); … … 727 845 } 728 846 847 if ((stamp->x <= 0.0) && (stamp->y <= 0.0)) { 848 psAbort ("bad stamp"); 849 } 850 if (!isfinite(stamp->x) && !isfinite(stamp->y)) { 851 psAbort ("bad stamp"); 852 } 853 729 854 if (pmSubtractionThreaded()) { 730 855 psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION"); … … 732 857 psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array 733 858 PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); 859 PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32); 734 860 if (!psThreadJobAddPending(job)) { 735 861 psFree(job); … … 738 864 psFree(job); 739 865 } else { 740 pmSubtractionCalculateEquationStamp(stamps, kernels, i );866 pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode); 741 867 } 742 868 } … … 748 874 749 875 pmSubtractionVisualPlotLeastSquares(stamps); 876 pmSubtractionVisualShowKernels((pmSubtractionKernels *)kernels); 877 pmSubtractionVisualShowBasis(stamps); 750 878 751 879 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", … … 756 884 } 757 885 758 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) 886 // private functions used on pmSubtractionSolveEquation 887 bool psVectorWriteFile (char *filename, const psVector *vector); 888 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header); 889 890 psImage *p_pmSubSolve_wUt (psVector *w, psImage *U); 891 psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt); 892 893 bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask); 894 895 bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B); 896 bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB); 897 bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB); 898 899 bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x); 900 bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y); 901 bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); 902 903 psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w); 904 905 double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); 906 907 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 908 const pmSubtractionStampList *stamps, 909 const pmSubtractionEquationCalculationMode mode) 759 910 { 760 911 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); … … 762 913 763 914 // Check inputs 764 int numParams = -1; // Number of parameters 765 int numParams2 = 0; // Number of parameters for part solution (DUAL mode) 915 int numKernels = kernels->num; // Number of kernel basis functions 916 int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations 917 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 918 int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for 919 int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution 920 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 921 // An additional image is convolved 922 numSolution2 = numKernels * numSpatial; 923 numParams += numSolution2; 924 } 925 766 926 for (int i = 0; i < stamps->num; i++) { 767 927 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 771 931 } 772 932 773 PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false); 774 if (numParams == -1) { 775 numParams = stamp->vector1->n; 776 } 777 PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false); 778 PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false); 779 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); 780 PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); 781 PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); 782 783 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 784 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false); 785 PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); 786 if (numParams2 == 0) { 787 numParams2 = stamp->matrix2->numCols; 788 } 789 PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); 790 PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false); 791 PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false); 792 PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false); 793 PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false); 794 PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false); 795 PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false); 796 } 797 } 798 if (numParams == -1) { 799 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); 800 return NULL; 933 PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); 934 PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); 935 PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); 936 PS_ASSERT_IMAGE_NON_NULL(stamp->matrix, false); 937 PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, false); 938 PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, false); 801 939 } 802 940 … … 814 952 psVectorInit(sumVector, 0.0); 815 953 psImageInit(sumMatrix, 0.0); 954 955 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 956 816 957 int numStamps = 0; // Number of good stamps 817 958 for (int i = 0; i < stamps->num; i++) { 818 959 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 819 820 960 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 821 822 #ifdef TESTING 823 // XXX double-check for NAN in data: 824 for (int iy = 0; iy < stamp->matrix1->numRows; iy++) { 825 for (int ix = 0; ix < stamp->matrix1->numCols; ix++) { 826 if (!isfinite(stamp->matrix1->data.F64[iy][ix])) { 827 fprintf (stderr, "WARNING: NAN in matrix1\n"); 828 } 829 } 830 } 831 for (int ix = 0; ix < stamp->vector1->n; ix++) { 832 if (!isfinite(stamp->vector1->data.F64[ix])) { 833 fprintf (stderr, "WARNING: NAN in vector1\n"); 834 } 835 } 836 #endif 837 838 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); 839 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); 961 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); 962 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 963 psVectorAppend(norms, stamp->norm); 840 964 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 841 965 numStamps++; … … 845 969 } 846 970 847 #ifdef TESTING848 for (int ix = 0; ix < sumVector->n; ix++) {849 if (!isfinite(sumVector->data.F64[ix])) {850 fprintf (stderr, "WARNING: NAN in vector1\n");851 }852 }853 #endif854 855 971 #if 0 856 972 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background … … 858 974 #endif 859 975 860 #ifdef TESTING 861 for (int ix = 0; ix < sumVector->n; ix++) { 862 if (!isfinite(sumVector->data.F64[ix])) { 863 fprintf (stderr, "WARNING: NAN in vector1\n"); 864 } 865 } 976 psVector *solution = NULL; // Solution to equation! 977 solution = psVectorAlloc(numParams, PS_TYPE_F64); 978 psVectorInit(solution, 0); 979 980 #if 0 981 // Regular, straight-forward solution 982 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 983 #else 866 984 { 867 psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); 868 psFits *fits = psFitsOpen("matrixInv.fits", "w"); 869 psFitsWriteImage(fits, NULL, inverse, 0, NULL); 870 psFitsClose(fits); 871 psFree(inverse); 872 } 873 { 874 psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); 875 psImage *Xt = psMatrixTranspose(NULL, X); 876 psImage *XtX = psMatrixMultiply(NULL, Xt, X); 877 psFits *fits = psFitsOpen("matrixErr.fits", "w"); 878 psFitsWriteImage(fits, NULL, XtX, 0, NULL); 879 psFitsClose(fits); 880 psFree(X); 881 psFree(Xt); 882 psFree(XtX); 883 } 884 #endif 885 886 psVector *permutation = NULL; // Permutation vector, required for LU decomposition 887 psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); 985 // Solve normalisation and background separately 986 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 987 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 988 989 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 990 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 991 psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); 992 psFree(stats); 993 psFree(sumMatrix); 994 psFree(sumVector); 995 psFree(norms); 996 return false; 997 } 998 999 double normValue = stats->robustMedian; 1000 // double bgValue = 0.0; 1001 1002 psFree(stats); 1003 1004 fprintf(stderr, "Norm: %lf\n", normValue); 1005 1006 // Solve kernel components 1007 for (int i = 0; i < numSolution1; i++) { 1008 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1009 1010 sumMatrix->data.F64[i][normIndex] = 0.0; 1011 sumMatrix->data.F64[normIndex][i] = 0.0; 1012 } 1013 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1014 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1015 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1016 1017 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1018 sumVector->data.F64[normIndex] = 0.0; 1019 1020 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1021 1022 solution->data.F64[normIndex] = normValue; 1023 } 1024 # endif 1025 1026 if (!kernels->solution1) { 1027 kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); 1028 psVectorInit(kernels->solution1, 0.0); 1029 } 1030 1031 // only update the solutions that we chose to calculate: 1032 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1033 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1034 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1035 } 1036 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1037 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1038 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1039 } 1040 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1041 int numKernels = kernels->num; 1042 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 1043 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 1044 for (int i = 0; i < numKernels * numPoly; i++) { 1045 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1046 } 1047 } 1048 1049 psFree(solution); 1050 psFree(sumVector); 888 1051 psFree(sumMatrix); 889 if (!luMatrix) {890 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");891 psFree(sumVector);892 psFree(luMatrix);893 psFree(permutation);894 return NULL;895 }896 kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);897 1052 898 1053 #ifdef TESTING … … 900 1055 for (int ix = 0; ix < kernels->solution1->n; ix++) { 901 1056 if (!isfinite(kernels->solution1->data.F64[ix])) { 902 fprintf (stderr, "WARNING: NAN in vector1\n"); 903 } 904 } 905 #endif 906 907 psFree(sumVector); 908 psFree(luMatrix); 909 psFree(permutation); 910 if (!kernels->solution1) { 911 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); 912 return NULL; 913 } 1057 fprintf (stderr, "WARNING: NAN in vector\n"); 1058 } 1059 } 1060 #endif 1061 914 1062 } else { 915 1063 // Dual convolution solution 916 1064 917 1065 // Accumulation of stamp matrices/vectors 918 psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); 919 psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); 920 psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); 921 psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); 922 psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64); 923 psImageInit(sumMatrix1, 0.0); 924 psImageInit(sumMatrix2, 0.0); 925 psImageInit(sumMatrixX, 0.0); 926 psVectorInit(sumVector1, 0.0); 927 psVectorInit(sumVector2, 0.0); 1066 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); 1067 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); 1068 psImageInit(sumMatrix, 0.0); 1069 psVectorInit(sumVector, 0.0); 1070 1071 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 928 1072 929 1073 int numStamps = 0; // Number of good stamps … … 931 1075 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 932 1076 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 933 (void)psBinaryOp(sumMatrix 1, sumMatrix1, "+", stamp->matrix1);934 (void)psBinaryOp(sum Matrix2, sumMatrix2, "+", stamp->matrix2);935 (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); 936 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);937 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); 1077 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); 1078 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 1079 1080 psVectorAppend(norms, stamp->norm); 1081 938 1082 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 939 1083 numStamps++; … … 941 1085 } 942 1086 943 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background944 calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);945 calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]);946 947 psImage *sumMatrix = NULL; // Combined matrix948 psVector *sumVector = NULL; // Combined vector949 calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,950 sumMatrixX, sumVector1, sumVector2);951 952 1087 #ifdef TESTING 1088 psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); 1089 psVectorWriteFile("sumVector.dat", sumVector); 1090 #endif 1091 1092 #if 1 1093 // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1094 // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); 1095 1096 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1097 calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 1000.0); 1098 #endif 1099 1100 psVector *solution = NULL; // Solution to equation! 1101 solution = psVectorAlloc(numParams, PS_TYPE_F64); 1102 psVectorInit(solution, 0); 1103 1104 #if 0 1105 // Regular, straight-forward solution 1106 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1107 #else 953 1108 { 954 psFits *fits = psFitsOpen("sumMatrix.fits", "w"); 955 psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); 956 psFitsClose(fits); 957 } 958 { 959 psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); 960 psFits *fits = psFitsOpen("sumVector.fits", "w"); 961 for (int i = 0; i < numParams + numParams2; i++) { 962 image->data.F64[0][i] = sumVector->data.F64[i]; 963 } 964 psFitsWriteImage(fits, NULL, image, 0, NULL); 965 psFree(image); 966 psFitsClose(fits); 967 } 968 #endif 969 970 psVector *solution = NULL; // Solution to equation! 971 { 972 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); 973 if (!solution) { 974 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); 1109 // Solve normalisation and background separately 1110 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1111 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1112 1113 #if 0 1114 psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); 1115 psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); 1116 1117 normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; 1118 normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; 1119 normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; 1120 1121 normVector->data.F64[0] = sumVector->data.F64[normIndex]; 1122 normVector->data.F64[1] = sumVector->data.F64[bgIndex]; 1123 1124 psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); 1125 1126 double normValue = normSolution->data.F64[0]; 1127 double bgValue = normSolution->data.F64[1]; 1128 1129 psFree(normMatrix); 1130 psFree(normVector); 1131 psFree(normSolution); 1132 #endif 1133 1134 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 1135 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 1136 psError(PS_ERR_UNKNOWN, false, "Unable to determine median normalisation"); 1137 psFree(stats); 975 1138 psFree(sumMatrix); 976 1139 psFree(sumVector); 977 return NULL; 978 } 979 980 int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations 981 int numKernels = kernels->num; // Number of kernel basis functions 982 983 // Remove a kernel basis for image 1 from the equation 984 #define MASK_BASIS_1(INDEX) \ 985 { \ 986 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ 987 for (int k = 0; k < numParams2; k++) { \ 988 sumMatrix1->data.F64[k][index] = 0.0; \ 989 sumMatrix1->data.F64[index][k] = 0.0; \ 990 sumMatrixX->data.F64[k][index] = 0.0; \ 991 } \ 992 sumMatrix1->data.F64[bgIndex][index] = 0.0; \ 993 sumMatrix1->data.F64[index][bgIndex] = 0.0; \ 994 sumMatrix1->data.F64[normIndex][index] = 0.0; \ 995 sumMatrix1->data.F64[index][normIndex] = 0.0; \ 996 sumMatrix1->data.F64[index][index] = 1.0; \ 997 sumVector1->data.F64[index] = 0.0; \ 998 } \ 999 } 1000 1001 // Remove a kernel basis for image 2 from the equation 1002 #define MASK_BASIS_2(INDEX) \ 1003 { \ 1004 for (int j = 0, index = INDEX; j < numSpatial; j++, index += numKernels) { \ 1005 for (int k = 0; k < numParams2; k++) { \ 1006 sumMatrix2->data.F64[k][index] = 0.0; \ 1007 sumMatrix2->data.F64[index][k] = 0.0; \ 1008 sumMatrixX->data.F64[index][k] = 0.0; \ 1009 } \ 1010 sumMatrix2->data.F64[index][index] = 1.0; \ 1011 sumMatrixX->data.F64[index][normIndex] = 0.0; \ 1012 sumMatrixX->data.F64[index][bgIndex] = 0.0; \ 1013 sumVector2->data.F64[index] = 0.0; \ 1014 } \ 1015 } 1016 1017 #define TOL 1.0e-5 1018 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1019 double norm = fabs(solution->data.F64[normIndex]); // Normalisation 1020 double thresh = norm * TOL; // Threshold for low parameters 1021 for (int i = 0; i < numKernels; i++) { 1022 // Getting 0th order parameter value. In the presence of spatial variation, the actual value 1023 // of the parameter will vary over the image. We are in effect getting the value in the 1024 // centre of the image. If we use different polynomial functions (e.g., Chebyshev), we may 1025 // have to change this to properly determine the value of the parameter at the centre. 1026 double param1 = solution->data.F64[i], 1027 param2 = solution->data.F64[numParams + i]; // Parameters of interest 1028 bool mask1 = false, mask2 = false; // Masked the parameter? 1029 if (fabs(param1) < thresh) { 1030 psTrace("psModules.imcombine", 7, "Parameter %d: 1 below threshold\n", i); 1031 MASK_BASIS_1(i); 1032 mask1 = true; 1033 } 1034 if (fabs(param2) < thresh) { 1035 psTrace("psModules.imcombine", 7, "Parameter %d: 2 below threshold\n", i); 1036 MASK_BASIS_2(i); 1037 mask2 = true; 1038 } 1039 1040 if (!mask1 && !mask2) { 1041 if (fabs(param1) < fabs(param2)) { 1042 psTrace("psModules.imcombine", 7, "Parameter %d: 1 < 2\n", i); 1043 MASK_BASIS_1(i); 1044 } else { 1045 psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", i); 1046 MASK_BASIS_2(i); 1047 } 1048 } 1049 } 1050 } 1051 1052 calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2, 1053 sumMatrixX, sumVector1, sumVector2); 1140 psFree(norms); 1141 return false; 1142 } 1143 1144 double normValue = stats->robustMedian; 1145 1146 psFree(stats); 1147 1148 fprintf(stderr, "Norm: %lf\n", normValue); 1149 1150 // Solve kernel components 1151 for (int i = 0; i < numSolution2; i++) { 1152 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1153 sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; 1154 1155 sumMatrix->data.F64[i][normIndex] = 0.0; 1156 sumMatrix->data.F64[normIndex][i] = 0.0; 1157 1158 sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; 1159 sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; 1160 } 1161 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1162 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1163 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1164 1165 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1166 1167 sumVector->data.F64[normIndex] = 0.0; 1168 1169 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1170 1171 solution->data.F64[normIndex] = normValue; 1172 } 1173 #endif 1174 1054 1175 1055 1176 #ifdef TESTING 1056 { 1057 psFits *fits = psFitsOpen("sumMatrixFix.fits", "w"); 1058 psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL); 1059 psFitsClose(fits); 1060 } 1061 { 1062 psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); 1063 psFits *fits = psFitsOpen("sumVectorFix.fits", "w"); 1064 for (int i = 0; i < numParams + numParams2; i++) { 1065 image->data.F64[0][i] = sumVector->data.F64[i]; 1066 } 1067 psFitsWriteImage(fits, NULL, image, 0, NULL); 1068 psFree(image); 1069 psFitsClose(fits); 1070 } 1071 #endif 1072 1073 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector); 1074 if (!solution) { 1075 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n"); 1076 psFree(sumMatrix); 1077 psFree(sumVector); 1078 return NULL; 1079 } 1080 1081 psFree(sumMatrix1); 1082 psFree(sumMatrix2); 1083 psFree(sumMatrixX); 1084 psFree(sumVector1); 1085 psFree(sumVector2); 1177 for (int i = 0; i < solution->n; i++) { 1178 fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); 1179 } 1180 #endif 1086 1181 1087 1182 psFree(sumMatrix); 1088 1183 psFree(sumVector); 1089 1184 1090 #ifdef TESTING 1091 { 1092 psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64); 1093 psFits *fits = psFitsOpen("solnVector.fits", "w"); 1094 for (int i = 0; i < numParams + numParams2; i++) { 1095 image->data.F64[0][i] = solution->data.F64[i]; 1096 } 1097 psFitsWriteImage(fits, NULL, image, 0, NULL); 1098 psFree(image); 1099 psFitsClose(fits); 1100 } 1101 #endif 1185 psFree(norms); 1102 1186 1103 1187 if (!kernels->solution1) { 1104 kernels->solution1 = psVectorAlloc(numParams, PS_TYPE_F64); 1188 kernels->solution1 = psVectorAlloc(numSolution1, PS_TYPE_F64); 1189 psVectorInit (kernels->solution1, 0.0); 1105 1190 } 1106 1191 if (!kernels->solution2) { 1107 kernels->solution2 = psVectorAlloc(numParams2, PS_TYPE_F64); 1108 } 1109 1110 memcpy(kernels->solution1->data.F64, solution->data.F64, numParams * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1111 memcpy(kernels->solution2->data.F64, &solution->data.F64[numParams], 1112 numParams2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1192 kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); 1193 psVectorInit (kernels->solution2, 0.0); 1194 } 1195 1196 // only update the solutions that we chose to calculate: 1197 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1198 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1199 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1200 } 1201 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1202 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1203 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1204 } 1205 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1206 int numKernels = kernels->num; 1207 for (int i = 0; i < numKernels * numSpatial; i++) { 1208 // XXX fprintf (stderr, "keep\n"); 1209 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1210 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; 1211 } 1212 } 1213 1214 1215 memcpy(kernels->solution1->data.F64, solution->data.F64, 1216 numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1217 memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], 1218 numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1113 1219 1114 1220 psFree(solution); … … 1131 1237 } 1132 1238 1133 pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const1239 // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const 1134 1240 return true; 1135 1241 } 1136 1242 1243 bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { 1244 1245 // XXX measure some useful stats on the residuals 1246 float sum = 0.0; 1247 float peak = 0.0; 1248 for (int y = - footprint; y <= footprint; y++) { 1249 for (int x = - footprint; x <= footprint; x++) { 1250 sum += 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); 1251 peak = PS_MAX(peak, 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm)); 1252 } 1253 } 1254 1255 // only count pixels with more than X% of the source flux 1256 // calculate stdev(dflux) 1257 float dflux1 = 0.0; 1258 float dflux2 = 0.0; 1259 int npix = 0; 1260 1261 float dmax = 0.0; 1262 float dmin = 0.0; 1263 1264 for (int y = - footprint; y <= footprint; y++) { 1265 for (int x = - footprint; x <= footprint; x++) { 1266 float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm); 1267 if (dflux < 0.02*sum) continue; 1268 dflux1 += residual->kernel[y][x]; 1269 dflux2 += PS_SQR(residual->kernel[y][x]); 1270 dmax = PS_MAX(residual->kernel[y][x], dmax); 1271 dmin = PS_MIN(residual->kernel[y][x], dmin); 1272 npix ++; 1273 } 1274 } 1275 float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); 1276 if (!isfinite(sum)) return false; 1277 if (!isfinite(dmax)) return false; 1278 if (!isfinite(dmin)) return false; 1279 if (!isfinite(peak)) return false; 1280 1281 // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak); 1282 psVectorAppend(fSigRes, sigma/sum); 1283 psVectorAppend(fMaxRes, dmax/peak); 1284 psVectorAppend(fMinRes, dmin/peak); 1285 return true; 1286 } 1287 1137 1288 psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, 1138 constpmSubtractionKernels *kernels)1289 pmSubtractionKernels *kernels) 1139 1290 { 1140 1291 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); … … 1151 1302 psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image 1152 1303 1304 // set up holding images for the visualization 1305 pmSubtractionVisualShowFitInit (stamps); 1306 1307 psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1308 psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1309 psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1310 1311 // we want to save the residual images for the 9 brightest stamps. 1312 // identify the 9 brightest stamps 1313 psVector *keepStamps = psVectorAlloc(stamps->num, PS_TYPE_S32); 1314 psVectorInit (keepStamps, 0); 1315 { 1316 psVector *flux = psVectorAlloc(stamps->num, PS_TYPE_F32); 1317 psVectorInit (flux, 0.0); 1318 1319 for (int i = 0; i < stamps->num; i++) { 1320 pmSubtractionStamp *stamp = stamps->stamps->data[i]; 1321 if (!isfinite(stamp->flux)) continue; 1322 flux->data.F32[i] = stamp->flux; 1323 } 1324 1325 psVector *index = psVectorSortIndex(NULL, flux); 1326 for (int i = 0; (i < stamps->num) && (i < 9); i++) { 1327 int n = stamps->num - i - 1; 1328 keepStamps->data.S32[index->data.S32[n]] = 1; 1329 fprintf (stderr, "keeping %d (%d of %d)\n", index->data.S32[n], n, 9); 1330 } 1331 psFree (flux); 1332 psFree (index); 1333 1334 // this function is called multiple times in the iteration, but 1335 // we only know after the interation is done if we will try again. 1336 // therefore we must save the sample each time, and blow away the old one 1337 // if it exists. 1338 psFree (kernels->sampleStamps); 1339 kernels->sampleStamps = psArrayAllocEmpty(9); 1340 } 1341 1342 psString log = psStringCopy("Deviations:\n"); // Log message with deviations 1153 1343 for (int i = 0; i < stamps->num; i++) { 1154 1344 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest … … 1210 1400 for (int y = - footprint; y <= footprint; y++) { 1211 1401 for (int x = - footprint; x <= footprint; x++) { 1212 residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;1402 residual->kernel[y][x] += convolution->kernel[y][x] * coefficient; 1213 1403 } 1214 1404 } 1215 1405 } 1406 1407 // XXX visualize the target, source, convolution and residual 1408 pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); 1409 1216 1410 for (int y = - footprint; y <= footprint; y++) { 1217 1411 for (int x = - footprint; x <= footprint; x++) { 1218 residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm; 1219 } 1220 } 1412 residual->kernel[y][x] += background + source->kernel[y][x] * norm - target->kernel[y][x]; 1413 } 1414 } 1415 1416 if (keepStamps->data.S32[i]) { 1417 psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32); 1418 psArrayAdd (kernels->sampleStamps, 9, sample); 1419 psFree (sample); 1420 } 1421 1422 pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint); 1423 1221 1424 } else { 1222 1425 // Dual convolution … … 1234 1437 for (int y = - footprint; y <= footprint; y++) { 1235 1438 for (int x = - footprint; x <= footprint; x++) { 1236 residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 -conv1->kernel[y][x] * coeff1;1439 residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1; 1237 1440 } 1238 1441 } 1239 1442 } 1443 1444 // XXX visualize the target, source, convolution and residual 1445 pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); 1446 1240 1447 for (int y = - footprint; y <= footprint; y++) { 1241 1448 for (int x = - footprint; x <= footprint; x++) { 1242 residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm; 1243 } 1244 } 1449 residual->kernel[y][x] += background + image1->kernel[y][x] * norm - image2->kernel[y][x]; 1450 } 1451 } 1452 if (keepStamps->data.S32[i]) { 1453 psImage *sample = psImageCopy(NULL, residual->image, PS_TYPE_F32); 1454 psArrayAdd (kernels->sampleStamps, 9, sample); 1455 psFree (sample); 1456 } 1457 1458 pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint); 1245 1459 } 1246 1460 … … 1257 1471 deviations->data.F32[i] = devNorm * deviation; 1258 1472 psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", 1259 i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); 1473 i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]); 1474 psStringAppend(&log, "Stamp %d (%d,%d): %f\n", 1475 i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]); 1260 1476 if (!isfinite(deviations->data.F32[i])) { 1261 1477 stamp->status = PM_SUBTRACTION_STAMP_REJECTED; 1262 1478 psTrace("psModules.imcombine", 5, 1263 1479 "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", 1264 i, (int)(stamp->x + 0.5), (int)(stamp->y +0.5));1480 i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); 1265 1481 continue; 1266 1482 } … … 1302 1518 1303 1519 } 1520 1521 psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); 1522 psFree(log); 1523 1524 // calculate and report the normalization and background for the image center 1525 { 1526 polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0); 1527 double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation 1528 double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background 1529 psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); 1530 1531 pmSubtractionVisualShowFit(norm); 1532 pmSubtractionVisualPlotFit(kernels); 1533 1534 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); 1535 psVectorStats (stats, fSigRes, NULL, NULL, 0); 1536 kernels->fSigResMean = stats->robustMedian; 1537 kernels->fSigResStdev = stats->robustStdev; 1538 1539 psStatsInit (stats); 1540 psVectorStats (stats, fMaxRes, NULL, NULL, 0); 1541 kernels->fMaxResMean = stats->robustMedian; 1542 kernels->fMaxResStdev = stats->robustStdev; 1543 1544 psStatsInit (stats); 1545 psVectorStats (stats, fMinRes, NULL, NULL, 0); 1546 kernels->fMinResMean = stats->robustMedian; 1547 kernels->fMinResStdev = stats->robustStdev; 1548 1549 // XXX save these values somewhere 1550 psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f", 1551 kernels->fSigResMean, kernels->fSigResStdev, 1552 kernels->fMaxResMean, kernels->fMaxResStdev, 1553 kernels->fMinResMean, kernels->fMinResStdev); 1554 1555 psFree (fSigRes); 1556 psFree (fMaxRes); 1557 psFree (fMinRes); 1558 psFree (stats); 1559 } 1560 1304 1561 psFree(residual); 1305 1562 psFree(polyValues); 1306 1563 1564 1307 1565 return deviations; 1308 1566 } 1567 1568 // we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w) 1569 psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) { 1570 1571 psAssert (w->n == U->numCols, "w and U dimensions do not match"); 1572 1573 // wUt has dimensions transposed relative to Ut. 1574 psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64); 1575 psImageInit (wUt, 0.0); 1576 1577 for (int i = 0; i < wUt->numCols; i++) { 1578 for (int j = 0; j < wUt->numRows; j++) { 1579 if (!isfinite(w->data.F64[j])) continue; 1580 if (w->data.F64[j] == 0.0) continue; 1581 wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j]; 1582 } 1583 } 1584 return wUt; 1585 } 1586 1587 // XXX this is just standard matrix multiplication: use psMatrixMultiply? 1588 psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) { 1589 1590 psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match"); 1591 1592 psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64); 1593 1594 for (int i = 0; i < Ainv->numCols; i++) { 1595 for (int j = 0; j < Ainv->numRows; j++) { 1596 double sum = 0.0; 1597 for (int k = 0; k < V->numCols; k++) { 1598 sum += V->data.F64[j][k] * wUt->data.F64[k][i]; 1599 } 1600 Ainv->data.F64[j][i] = sum; 1601 } 1602 } 1603 return Ainv; 1604 } 1605 1606 // we are supplied U, not Ut 1607 bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) { 1608 1609 psAssert (U->numRows == B->n, "U and B dimensions do not match"); 1610 1611 UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64); 1612 1613 for (int i = 0; i < U->numCols; i++) { 1614 double sum = 0.0; 1615 for (int j = 0; j < U->numRows; j++) { 1616 sum += B->data.F64[j] * U->data.F64[j][i]; 1617 } 1618 UtB[0]->data.F64[i] = sum; 1619 } 1620 return true; 1621 } 1622 1623 // w is diagonal 1624 bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) { 1625 1626 psAssert (w->n == UtB->n, "w and UtB dimensions do not match"); 1627 1628 // wUt has dimensions transposed relative to Ut. 1629 wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64); 1630 psVectorInit (wUtB[0], 0.0); 1631 1632 for (int i = 0; i < w->n; i++) { 1633 if (!isfinite(w->data.F64[i])) continue; 1634 if (w->data.F64[i] == 0.0) continue; 1635 wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; 1636 } 1637 return true; 1638 } 1639 1640 // this is basically matrix * vector 1641 bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) { 1642 1643 psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match"); 1644 1645 VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64); 1646 1647 for (int j = 0; j < V->numRows; j++) { 1648 double sum = 0.0; 1649 for (int i = 0; i < V->numCols; i++) { 1650 sum += V->data.F64[j][i] * wUtB->data.F64[i]; 1651 } 1652 VwUtB[0]->data.F64[j] = sum; 1653 } 1654 return true; 1655 } 1656 1657 // this is basically matrix * vector 1658 bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) { 1659 1660 psAssert (A->numCols == x->n, "A and x dimensions do not match"); 1661 1662 B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64); 1663 1664 for (int j = 0; j < A->numRows; j++) { 1665 double sum = 0.0; 1666 for (int i = 0; i < A->numCols; i++) { 1667 sum += A->data.F64[j][i] * x->data.F64[i]; 1668 } 1669 B[0]->data.F64[j] = sum; 1670 } 1671 return true; 1672 } 1673 1674 // this is basically Vector * vector 1675 bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) { 1676 1677 psAssert (x->n == y->n, "x and y dimensions do not match"); 1678 1679 double sum = 0.0; 1680 for (int i = 0; i < x->n; i++) { 1681 sum += x->data.F64[i] * y->data.F64[i]; 1682 } 1683 *value = sum; 1684 return true; 1685 } 1686 1687 bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { 1688 1689 int footprint = stamps->footprint; // Half-size of stamps 1690 1691 double sum = 0.0; 1692 for (int i = 0; i < stamps->num; i++) { 1693 1694 pmSubtractionStamp *stamp = stamps->stamps->data[i]; 1695 if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; 1696 1697 psKernel *weight = NULL; 1698 psKernel *window = NULL; 1699 psKernel *input = NULL; 1700 1701 #ifdef USE_WEIGHT 1702 weight = stamp->weight; 1703 #endif 1704 #ifdef USE_WINDOW 1705 window = stamps->window; 1706 #endif 1707 1708 switch (kernels->mode) { 1709 // MODE_1 : convolve image 1 to match image 2 (and vice versa) 1710 case PM_SUBTRACTION_MODE_1: 1711 input = stamp->image2; 1712 break; 1713 case PM_SUBTRACTION_MODE_2: 1714 input = stamp->image1; 1715 break; 1716 default: 1717 psAbort ("programming error"); 1718 } 1719 1720 for (int y = - footprint; y <= footprint; y++) { 1721 for (int x = - footprint; x <= footprint; x++) { 1722 double in = input->kernel[y][x]; 1723 double value = in*in; 1724 if (weight) { 1725 float wtVal = weight->kernel[y][x]; 1726 value *= wtVal; 1727 } 1728 if (window) { 1729 float winVal = window->kernel[y][x]; 1730 value *= winVal; 1731 } 1732 sum += value; 1733 } 1734 } 1735 } 1736 *y2 = sum; 1737 return true; 1738 } 1739 1740 double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { 1741 1742 int footprint = stamps->footprint; // Half-size of stamps 1743 int numKernels = kernels->num; // Number of kernels 1744 1745 double sum = 0.0; 1746 1747 psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image 1748 psImageInit(residual->image, 0.0); 1749 1750 psImage *polyValues = NULL; // Polynomial values 1751 1752 for (int i = 0; i < stamps->num; i++) { 1753 1754 pmSubtractionStamp *stamp = stamps->stamps->data[i]; 1755 if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; 1756 1757 psKernel *weight = NULL; 1758 psKernel *window = NULL; 1759 psKernel *target = NULL; 1760 psKernel *source = NULL; 1761 1762 psArray *convolutions = NULL; 1763 1764 #ifdef USE_WEIGHT 1765 weight = stamp->weight; 1766 #endif 1767 #ifdef USE_WINDOW 1768 window = stamps->window; 1769 #endif 1770 1771 switch (kernels->mode) { 1772 // MODE_1 : convolve image 1 to match image 2 (and vice versa) 1773 case PM_SUBTRACTION_MODE_1: 1774 target = stamp->image2; 1775 source = stamp->image1; 1776 convolutions = stamp->convolutions1; 1777 break; 1778 case PM_SUBTRACTION_MODE_2: 1779 target = stamp->image1; 1780 source = stamp->image2; 1781 convolutions = stamp->convolutions2; 1782 break; 1783 default: 1784 psAbort ("programming error"); 1785 } 1786 1787 // Calculate coefficients of the kernel basis functions 1788 polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); 1789 double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation 1790 double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background 1791 1792 psImageInit(residual->image, 0.0); 1793 for (int j = 0; j < numKernels; j++) { 1794 psKernel *convolution = convolutions->data[j]; // Convolution 1795 double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1796 for (int y = - footprint; y <= footprint; y++) { 1797 for (int x = - footprint; x <= footprint; x++) { 1798 residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; 1799 } 1800 } 1801 } 1802 1803 for (int y = - footprint; y <= footprint; y++) { 1804 for (int x = - footprint; x <= footprint; x++) { 1805 double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x]; 1806 double value = PS_SQR(resid); 1807 if (weight) { 1808 float wtVal = weight->kernel[y][x]; 1809 value *= wtVal; 1810 } 1811 if (window) { 1812 float winVal = window->kernel[y][x]; 1813 value *= winVal; 1814 } 1815 sum += value; 1816 } 1817 } 1818 } 1819 psFree (polyValues); 1820 psFree (residual); 1821 1822 return sum; 1823 } 1824 1825 bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) { 1826 1827 for (int i = 0; i < w->n; i++) { 1828 wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i]; 1829 } 1830 return true; 1831 } 1832 1833 // we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w) 1834 psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) { 1835 1836 psAssert (w->n == V->numCols, "w and U dimensions do not match"); 1837 1838 psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); 1839 psImageInit (Vn, 0.0); 1840 1841 // generate Vn = V * w^{-1} 1842 for (int j = 0; j < Vn->numRows; j++) { 1843 for (int i = 0; i < Vn->numCols; i++) { 1844 if (!isfinite(w->data.F64[i])) continue; 1845 if (w->data.F64[i] == 0.0) continue; 1846 Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i]; 1847 } 1848 } 1849 1850 psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); 1851 psImageInit (Xvar, 0.0); 1852 1853 // generate Xvar = Vn * Vn^T 1854 for (int j = 0; j < Vn->numRows; j++) { 1855 for (int i = 0; i < Vn->numCols; i++) { 1856 double sum = 0.0; 1857 for (int k = 0; k < Vn->numCols; k++) { 1858 sum += Vn->data.F64[k][i]*Vn->data.F64[k][j]; 1859 } 1860 Xvar->data.F64[j][i] = sum; 1861 } 1862 } 1863 return Xvar; 1864 } 1865 1866 // I get confused by the index values between the image vs matrix usage: In terms 1867 // of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix 1868 // multiplication is: A_k,j * B_i,k = C_i,j 1869 1870 1871 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) { 1872 1873 psFits *fits = psFitsOpen(filename, "w"); 1874 psFitsWriteImage(fits, header, image, 0, NULL); 1875 psFitsClose(fits); 1876 1877 return true; 1878 } 1879 1880 bool psVectorWriteFile (char *filename, const psVector *vector) { 1881 1882 FILE *f = fopen (filename, "w"); 1883 int fd = fileno(f); 1884 p_psVectorPrint (fd, vector, "unnamed"); 1885 fclose (f); 1886 1887 return true; 1888 } 1889 1890 1891 # if 0 1892 1893 #ifdef TESTING 1894 psFitsWriteImageSimple("A.fits", sumMatrix, NULL); 1895 psVectorWriteFile ("B.dat", sumVector); 1896 #endif 1897 1898 # define SVD_ANALYSIS 0 1899 # define COEFF_SIG 0.0 1900 # define SVD_TOL 0.0 1901 1902 // Use SVD to determine the kernel coeffs (and validate) 1903 if (SVD_ANALYSIS) { 1904 1905 // We have sumVector and sumMatrix. we are trying to solve the following equation: 1906 // sumMatrix * x = sumVector. 1907 1908 // we can use any standard matrix inversion to solve this. However, the basis 1909 // functions in general have substantial correlation, so that the solution may be 1910 // somewhat poorly determined or unstable. If not numerically ill-conditioned, the 1911 // system of equations may be statistically ill-conditioned. Noise in the image 1912 // will drive insignificant, but correlated, terms in the solution. To avoid these 1913 // problems, we can use SVD to identify numerically unconstrained values and to 1914 // avoid statistically badly determined value. 1915 1916 // A = sumMatrix, B = sumVector 1917 // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T 1918 // x = V (1/w) (U^T B) 1919 // \sigma_x = sqrt(diag(A^{-1})) 1920 // solve for x and A^{-1} to get x & dx 1921 // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.0 1922 // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.0 1923 1924 // If I use the SVD trick to re-condition the matrix, I need to break out the 1925 // kernel and normalization terms from the background term. 1926 // XXX is this true? or was this due to an error in the analysis? 1927 1928 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1929 1930 // now pull out the kernel elements into their own square matrix 1931 psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64); 1932 psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64); 1933 1934 for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) { 1935 if (ix == bgIndex) continue; 1936 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) { 1937 if (iy == bgIndex) continue; 1938 kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix]; 1939 ky++; 1940 } 1941 kernelVector->data.F64[kx] = sumVector->data.F64[ix]; 1942 kx++; 1943 } 1944 1945 psImage *U = NULL; 1946 psImage *V = NULL; 1947 psVector *w = NULL; 1948 if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) { 1949 psError(PS_ERR_UNKNOWN, false, "failed to perform SVD on sumMatrix\n"); 1950 return NULL; 1951 } 1952 1953 // calculate A_inverse: 1954 // Ainv = V * w * U^T 1955 psImage *wUt = p_pmSubSolve_wUt (w, U); 1956 psImage *Ainv = p_pmSubSolve_VwUt (V, wUt); 1957 psImage *Xvar = NULL; 1958 psFree (wUt); 1959 1960 # ifdef TESTING 1961 // kernel terms: 1962 for (int i = 0; i < w->n; i++) { 1963 fprintf (stderr, "w: %f\n", w->data.F64[i]); 1964 } 1965 # endif 1966 // loop over w adding in more and more of the values until chisquare is no longer 1967 // dropping significantly. 1968 // XXX this does not seem to work very well: we seem to need all terms even for 1969 // simple cases... 1970 1971 psVector *Xsvd = NULL; 1972 { 1973 psVector *Ax = NULL; 1974 psVector *UtB = NULL; 1975 psVector *wUtB = NULL; 1976 1977 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64); 1978 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8); 1979 psVectorInit (wMask, 1); // start by masking everything 1980 1981 double chiSquareLast = NAN; 1982 int maxWeight = 0; 1983 1984 double Axx, Bx, y2; 1985 1986 // XXX this is an attempt to exclude insignificant modes. 1987 // it was not successful with the ISIS kernel set: removing even 1988 // the least significant mode leaves additional ringing / noise 1989 // because the terms are so coupled. 1990 for (int k = 0; false && (k < w->n); k++) { 1991 1992 // unmask the k-th weight 1993 wMask->data.U8[k] = 0; 1994 p_pmSubSolve_SetWeights(wApply, w, wMask); 1995 1996 // solve for x: 1997 // x = V * w * (U^T * B) 1998 p_pmSubSolve_UtB (&UtB, U, kernelVector); 1999 p_pmSubSolve_wUtB (&wUtB, wApply, UtB); 2000 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); 2001 2002 // chi-square for this system of equations: 2003 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 2004 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 2005 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); 2006 p_pmSubSolve_VdV (&Axx, Ax, Xsvd); 2007 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); 2008 p_pmSubSolve_y2 (&y2, kernels, stamps); 2009 2010 // apparently, this works (compare with the brute force value below 2011 double chiSquare = Axx - 2.0*Bx + y2; 2012 double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare; 2013 chiSquareLast = chiSquare; 2014 2015 // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi); 2016 if (k && !maxWeight && (deltaChi < 1.0)) { 2017 maxWeight = k; 2018 } 2019 } 2020 2021 // keep all terms or we get extra ringing 2022 maxWeight = w->n; 2023 psVectorInit (wMask, 1); 2024 for (int k = 0; k < maxWeight; k++) { 2025 wMask->data.U8[k] = 0; 2026 } 2027 p_pmSubSolve_SetWeights(wApply, w, wMask); 2028 2029 // solve for x: 2030 // x = V * w * (U^T * B) 2031 p_pmSubSolve_UtB (&UtB, U, kernelVector); 2032 p_pmSubSolve_wUtB (&wUtB, wApply, UtB); 2033 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB); 2034 2035 // chi-square for this system of equations: 2036 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^2 2037 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^2 2038 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd); 2039 p_pmSubSolve_VdV (&Axx, Ax, Xsvd); 2040 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd); 2041 p_pmSubSolve_y2 (&y2, kernels, stamps); 2042 2043 // apparently, this works (compare with the brute force value below 2044 double chiSquare = Axx - 2.0*Bx + y2; 2045 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare); 2046 2047 // re-calculate A^{-1} to get new variances: 2048 // Ainv = V * w * U^T 2049 // XXX since we keep all terms, this is identical to Ainv 2050 psImage *wUt = p_pmSubSolve_wUt (wApply, U); 2051 Xvar = p_pmSubSolve_VwUt (V, wUt); 2052 psFree (wUt); 2053 2054 psFree (Ax); 2055 psFree (UtB); 2056 psFree (wUtB); 2057 psFree (wApply); 2058 psFree (wMask); 2059 } 2060 2061 // copy the kernel solutions to the full solution vector: 2062 solution = psVectorAlloc(sumVector->n, PS_TYPE_F64); 2063 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); 2064 2065 for (int ix = 0, kx = 0; ix < sumVector->n; ix++) { 2066 if (ix == bgIndex) { 2067 solution->data.F64[ix] = 0; 2068 solutionErr->data.F64[ix] = 0.001; 2069 continue; 2070 } 2071 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]); 2072 solution->data.F64[ix] = Xsvd->data.F64[kx]; 2073 kx++; 2074 } 2075 2076 psFree (kernelMatrix); 2077 psFree (kernelVector); 2078 2079 psFree (U); 2080 psFree (V); 2081 psFree (w); 2082 2083 psFree (Ainv); 2084 psFree (Xsvd); 2085 } else { 2086 psVector *permutation = NULL; // Permutation vector, required for LU decomposition 2087 psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix); 2088 if (!luMatrix) { 2089 psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); 2090 psFree(solution); 2091 psFree(sumVector); 2092 psFree(sumMatrix); 2093 psFree(luMatrix); 2094 psFree(permutation); 2095 return NULL; 2096 } 2097 2098 solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation); 2099 psFree(luMatrix); 2100 psFree(permutation); 2101 if (!solution) { 2102 psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); 2103 psFree(solution); 2104 psFree(sumVector); 2105 psFree(sumMatrix); 2106 return NULL; 2107 } 2108 2109 // XXX LUD does not provide A^{-1}? fake the error for now 2110 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64); 2111 for (int ix = 0; ix < sumVector->n; ix++) { 2112 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix]; 2113 } 2114 } 2115 2116 if (!kernels->solution1) { 2117 kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64); 2118 psVectorInit (kernels->solution1, 0.0); 2119 } 2120 2121 // only update the solutions that we chose to calculate: 2122 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 2123 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 2124 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 2125 } 2126 if (mode & PM_SUBTRACTION_EQUATION_BG) { 2127 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 2128 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 2129 } 2130 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 2131 int numKernels = kernels->num; 2132 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 2133 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 2134 for (int i = 0; i < numKernels * numPoly; i++) { 2135 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i])); 2136 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) { 2137 // XXX fprintf (stderr, "drop\n"); 2138 kernels->solution1->data.F64[i] = 0.0; 2139 } else { 2140 // XXX fprintf (stderr, "keep\n"); 2141 kernels->solution1->data.F64[i] = solution->data.F64[i]; 2142 } 2143 } 2144 } 2145 // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps); 2146 // fprintf (stderr, "chi square Brute = %f\n", chiSquare); 2147 2148 psFree(solution); 2149 psFree(sumVector); 2150 psFree(sumMatrix); 2151 # endif 2152 2153 #ifdef TESTING 2154 // XXX double-check for NAN in data: 2155 for (int iy = 0; iy < stamp->matrix->numRows; iy++) { 2156 for (int ix = 0; ix < stamp->matrix->numCols; ix++) { 2157 if (!isfinite(stamp->matrix->data.F64[iy][ix])) { 2158 fprintf (stderr, "WARNING: NAN in matrix\n"); 2159 } 2160 } 2161 } 2162 for (int ix = 0; ix < stamp->vector->n; ix++) { 2163 if (!isfinite(stamp->vector->data.F64[ix])) { 2164 fprintf (stderr, "WARNING: NAN in vector\n"); 2165 } 2166 } 2167 #endif 2168 2169 #ifdef TESTING 2170 for (int ix = 0; ix < sumVector->n; ix++) { 2171 if (!isfinite(sumVector->data.F64[ix])) { 2172 fprintf (stderr, "WARNING: NAN in vector\n"); 2173 } 2174 } 2175 #endif 2176 2177 #ifdef TESTING 2178 for (int ix = 0; ix < sumVector->n; ix++) { 2179 if (!isfinite(sumVector->data.F64[ix])) { 2180 fprintf (stderr, "WARNING: NAN in vector\n"); 2181 } 2182 } 2183 { 2184 psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL); 2185 psFitsWriteImageSimple("matrixInv.fits", inverse, NULL); 2186 psFree(inverse); 2187 } 2188 { 2189 psImage *X = psMatrixInvert(NULL, sumMatrix, NULL); 2190 psImage *Xt = psMatrixTranspose(NULL, X); 2191 psImage *XtX = psMatrixMultiply(NULL, Xt, X); 2192 psFitsWriteImageSimple("matrixErr.fits", XtX, NULL); 2193 psFree(X); 2194 psFree(Xt); 2195 psFree(XtX); 2196 } 2197 #endif 2198
Note:
See TracChangeset
for help on using the changeset viewer.
