Changeset 29543 for trunk/psModules/src/imcombine/pmSubtractionEquation.c
- Timestamp:
- Oct 25, 2010, 2:45:41 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/imcombine/pmSubtractionEquation.c
r29004 r29543 19 19 20 20 //#define USE_WEIGHT // Include weight (1/variance) in equation? 21 //#define USE_WINDOW // Include weight (1/variance) in equation? 22 21 22 // XXX TEST: 23 //# define USE_WINDOW // window to avoid neighbor contamination 24 25 # define PENALTY false 26 # define MOMENTS (!PENALTY) 23 27 24 28 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 27 31 28 32 // Calculate the least-squares matrix and vector 29 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated30 psVector *vector, // Least-squares vector, updated31 double *norm, // Normalisation, updated33 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated 34 psVector *vector, // Least-squares vector, updated 35 double normValue, // Normalisation, supplied 32 36 const psKernel *input, // Input image (target) 33 37 const psKernel *reference, // Reference image (convolution source) … … 37 41 const pmSubtractionKernels *kernels, // Kernels 38 42 const psImage *polyValues, // Spatial polynomial values 39 int footprint, // (Half-)Size of stamp 40 int normWindow1, // Window (half-)size for normalisation measurement 41 int normWindow2, // Window (half-)size for normalisation measurement 42 const pmSubtractionEquationCalculationMode mode 43 int footprint // (Half-)Size of stamp 43 44 ) 44 45 { … … 50 51 51 52 int numKernels = kernels->num; // Number of kernels 52 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation53 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background54 53 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 55 54 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 56 55 double poly[numPoly]; // Polynomial terms 57 56 double poly2[numPoly][numPoly]; // Polynomial-polynomial values 57 int numParams = numKernels * numPoly; 58 59 psAssert(matrix && 60 matrix->type.type == PS_TYPE_F64 && 61 matrix->numCols == numParams && 62 matrix->numRows == numParams, 63 "Least-squares matrix is bad."); 64 psAssert(vector && 65 vector->type.type == PS_TYPE_F64 && 66 vector->n == numParams, 67 "Least-squares vector is bad."); 58 68 59 69 // Evaluate polynomial-polynomial terms 60 // XXX we can skip this if we are not calculating kernel coeffs61 70 for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) { 62 71 for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) { … … 84 93 // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] 85 94 // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] 86 // normalization87 // bg 0, bg 1, bg 2 (only 0 is currently used?)88 95 89 96 for (int i = 0; i < numKernels; i++) { … … 107 114 108 115 // Spatial variation of kernel coeffs 109 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 110 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 111 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 112 double value = sumCC * poly2[iTerm][jTerm]; 113 matrix->data.F64[iIndex][jIndex] = value; 114 matrix->data.F64[jIndex][iIndex] = value; 115 } 116 } 117 } 116 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 117 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 118 double value = sumCC * poly2[iTerm][jTerm]; 119 matrix->data.F64[iIndex][jIndex] = value; 120 matrix->data.F64[jIndex][iIndex] = value; 121 } 122 } 118 123 } 119 124 120 125 double sumRC = 0.0; // Sum of the reference-convolution products 121 126 double sumIC = 0.0; // Sum of the input-convolution products 122 double sumC = 0.0; // Sum of the convolution123 127 for (int y = - footprint; y <= footprint; y++) { 124 128 for (int x = - footprint; x <= footprint; x++) { … … 128 132 double ic = in * conv; 129 133 double rc = ref * conv; 130 double c = conv;131 134 if (weight) { 132 135 float wtVal = weight->kernel[y][x]; 133 136 ic *= wtVal; 134 137 rc *= wtVal; 135 c *= wtVal;136 138 } 137 139 if (window) { … … 139 141 ic *= winVal; 140 142 rc *= winVal; 141 c *= winVal;142 143 } 143 144 sumIC += ic; 144 145 sumRC += rc; 145 sumC += c;146 147 } 146 } 147 } 148 148 149 // Spatial variation 149 150 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 150 double normTerm = sumRC * poly[iTerm]; 151 double bgTerm = sumC * poly[iTerm]; 152 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 153 matrix->data.F64[iIndex][normIndex] = normTerm; 154 matrix->data.F64[normIndex][iIndex] = normTerm; 155 } 156 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 157 matrix->data.F64[iIndex][bgIndex] = bgTerm; 158 matrix->data.F64[bgIndex][iIndex] = bgTerm; 159 } 160 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 161 vector->data.F64[iIndex] = sumIC * poly[iTerm]; 162 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 163 // subtract norm * sumRC * poly[iTerm] 164 psAssert (kernels->solution1, "programming error: define solution first!"); 165 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 166 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 167 vector->data.F64[iIndex] -= norm * normTerm; 168 } 169 } 170 } 171 } 172 173 double sumRR = 0.0; // Sum of the reference product 174 double sumIR = 0.0; // Sum of the input-reference product 175 double sum1 = 0.0; // Sum of the background 176 double sumR = 0.0; // Sum of the reference 177 double sumI = 0.0; // Sum of the input 178 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 179 for (int y = - footprint; y <= footprint; y++) { 180 for (int x = - footprint; x <= footprint; x++) { 181 double in = input->kernel[y][x]; 182 double ref = reference->kernel[y][x]; 183 double ir = in * ref; 184 double rr = PS_SQR(ref); 185 double one = 1.0; 186 187 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 188 normI1 += ref; 189 } 190 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 191 normI2 += in; 192 } 193 194 if (weight) { 195 float wtVal = weight->kernel[y][x]; 196 rr *= wtVal; 197 ir *= wtVal; 198 in *= wtVal; 199 ref *= wtVal; 200 one *= wtVal; 201 } 202 if (window) { 203 float winVal = window->kernel[y][x]; 204 rr *= winVal; 205 ir *= winVal; 206 in *= winVal; 207 ref *= winVal; 208 one *= winVal; 209 } 210 sumRR += rr; 211 sumIR += ir; 212 sumR += ref; 213 sumI += in; 214 sum1 += one; 215 } 216 } 217 218 *norm = normI2 / normI1; 219 220 fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); 221 222 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 223 matrix->data.F64[normIndex][normIndex] = sumRR; 224 vector->data.F64[normIndex] = sumIR; 225 // subtract sum over kernels * kernel solution 226 } 227 if (mode & PM_SUBTRACTION_EQUATION_BG) { 228 matrix->data.F64[bgIndex][bgIndex] = sum1; 229 vector->data.F64[bgIndex] = sumI; 230 } 231 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 232 matrix->data.F64[normIndex][bgIndex] = sumR; 233 matrix->data.F64[bgIndex][normIndex] = sumR; 151 vector->data.F64[iIndex] = (sumIC - normValue*sumRC) * poly[iTerm]; 152 } 234 153 } 235 154 … … 255 174 256 175 // Calculate the least-squares matrix and vector for dual convolution 257 static bool calculateDualMatrixVector(psImage *matrix, // Least-squares matrix, updated 258 psVector *vector, // Least-squares vector, updated 259 double *norm, // Normalisation, updated 260 const psKernel *image1, // Image 1 261 const psKernel *image2, // Image 2 176 // XXX we could avoid calculating these values on successive passes *if* the stamp has not changed. 177 static bool calculateDualMatrixVector(pmSubtractionStamp *stamp, // stamp of interest 178 double normValue, // Normalisation, updated 179 double normValue2, // Normalisation, updated 262 180 const psKernel *weight, // Weight image 263 181 const psKernel *window, // Window image 264 const psArray *convolutions1, // Convolutions of image 1 for each kernel265 const psArray *convolutions2, // Convolutions of image 2 for each kernel266 182 const pmSubtractionKernels *kernels, // Kernels 267 183 const psImage *polyValues, // Spatial polynomial values 268 int footprint, // (Half-)Size of stamp 269 int normWindow1, // Window (half-)size for normalisation measurement 270 int normWindow2, // Window (half-)size for normalisation measurement 271 const pmSubtractionEquationCalculationMode mode 184 int footprint // (Half-)Size of stamp 272 185 ) 273 186 { 274 int numKernels = kernels->num; // Number of kernels 275 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 276 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 277 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 187 int numKernels = kernels->num; // Number of kernels 188 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 278 189 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 279 190 double poly[numPoly]; // Polynomial terms 280 191 double poly2[numPoly][numPoly]; // Polynomial-polynomial values 281 192 282 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 283 int numParams = numKernels * numPoly + 1 + numBackground; // Number of regular parameters 284 int numParams2 = numKernels * numPoly; // Number of additional parameters for dual 285 int numDual = numParams + numParams2; // Total number of parameters for dual 193 int numParams = numKernels * numPoly; // Number of regular parameters 194 int numParams2 = numKernels * numPoly; // Number of additional parameters for dual 195 int numDual = numParams + numParams2; // Total number of parameters for dual 196 197 psImage *matrix = stamp->matrix; // Least-squares matrix, updated 198 psVector *vector = stamp->vector; // Least-squares vector, updated 199 psKernel *image1 = stamp->image1; // Image 1 200 psKernel *image2 = stamp->image2; // Image 2 201 psArray *convolutions1 = stamp->convolutions1; // Convolutions of image 1 for each kernel 202 psArray *convolutions2 = stamp->convolutions2; // Convolutions of image 2 for each kernel 286 203 287 204 psAssert(matrix && … … 290 207 matrix->numRows == numDual, 291 208 "Least-squares matrix is bad."); 209 292 210 psAssert(vector && 293 211 vector->type.type == PS_TYPE_F64 && … … 318 236 } 319 237 238 // the order of the elements in the matrix and vector is: 239 // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0] 240 // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0] 241 // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m] 242 243 // for DUAL convolution analysis, we apply the normalization to I1 as follows: 244 // norm = I2 / I1 245 // 246 // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i 247 // I2c = I2 + \Sum_i b_i I2 \cross k_i 248 249 // we cannot absorb the normalization into a_i until the analysis is complete, or the 250 // second moment terms are incorrectly calculated. 251 320 252 for (int i = 0; i < numKernels; i++) { 321 253 psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i … … 328 260 double sumBB = 0.0; // Sum of convolution products between image 2 329 261 double sumAB = 0.0; // Sum of convolution products across images 1 and 2 262 263 double MxxAA = 0.0; 264 double MyyAA = 0.0; 265 double MxxBB = 0.0; 266 double MyyBB = 0.0; 330 267 for (int y = - footprint; y <= footprint; y++) { 331 268 for (int x = - footprint; x <= footprint; x++) { 332 double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] ;269 double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue); 333 270 double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; 334 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] ;271 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; 335 272 if (weight) { 336 273 float wtVal = weight->kernel[y][x]; … … 348 285 sumBB += bb; 349 286 sumAB += ab; 350 } 351 } 287 288 if (MOMENTS) { 289 MxxAA += x*x*aa; 290 MyyAA += y*y*aa; 291 MxxBB += x*x*bb; 292 MyyBB += y*y*bb; 293 } 294 } 295 } 296 297 if (MOMENTS) { 298 MxxAA /= stamp->normSquare1 * PS_SQR(normValue); 299 MyyAA /= stamp->normSquare1 * PS_SQR(normValue); 300 MxxBB /= stamp->normSquare2; 301 MyyBB /= stamp->normSquare2; 302 } 303 304 // XXX this makes the Chisq portion independent of the normalization and star flux 305 // but may be mis-scaling between stars of different fluxes 306 sumAA /= PS_SQR(stamp->normI2); 307 sumAB /= PS_SQR(stamp->normI2); 308 sumBB /= PS_SQR(stamp->normI2); 309 310 // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB); 352 311 353 312 // Spatial variation of kernel coeffs 354 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 355 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 356 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 357 double aa = sumAA * poly2[iTerm][jTerm]; 358 double bb = sumBB * poly2[iTerm][jTerm]; 359 double ab = sumAB * poly2[iTerm][jTerm]; 360 361 matrix->data.F64[iIndex][jIndex] = aa; 362 matrix->data.F64[jIndex][iIndex] = aa; 363 364 matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; 365 matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; 366 367 matrix->data.F64[iIndex][jIndex + numParams] = ab; 368 matrix->data.F64[jIndex + numParams][iIndex] = ab; 369 } 370 } 371 } 372 } 313 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 314 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 315 double aa = sumAA * poly2[iTerm][jTerm]; 316 double bb = sumBB * poly2[iTerm][jTerm]; 317 double ab = sumAB * poly2[iTerm][jTerm]; 318 319 matrix->data.F64[iIndex][jIndex] = aa; 320 matrix->data.F64[jIndex][iIndex] = aa; 321 322 matrix->data.F64[iIndex + numParams][jIndex + numParams] = bb; 323 matrix->data.F64[jIndex + numParams][iIndex + numParams] = bb; 324 325 // add in second moments 326 if (MOMENTS) { 327 matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA; 328 matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA; 329 330 matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA; 331 matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA; 332 333 matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB; 334 matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB; 335 336 matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB; 337 matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB; 338 339 // XXX this uses the single moments, first try 340 // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; 341 // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; 342 // 343 // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j]; 344 // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j]; 345 // 346 // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; 347 // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; 348 // 349 // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j]; 350 // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j]; 351 } 352 matrix->data.F64[iIndex][jIndex + numParams] = ab; 353 matrix->data.F64[jIndex + numParams][iIndex] = ab; 354 } 355 } 356 } 357 358 // we need to calculate the lower-diagonal AB elements since they are not symmetric for A <-> B 373 359 for (int j = 0; j < i; j++) { 374 360 psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j … … 376 362 for (int y = - footprint; y <= footprint; y++) { 377 363 for (int x = - footprint; x <= footprint; x++) { 378 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] ;364 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; 379 365 if (weight) { 380 366 ab *= weight->kernel[y][x]; … … 387 373 } 388 374 375 // XXX this makes the Chisq portion independent of the normalization and star flux 376 // but may be mis-scaling between stars of different fluxes 377 sumAB /= PS_SQR(stamp->normI2); 378 389 379 // Spatial variation of kernel coeffs 390 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 391 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 392 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 393 double ab = sumAB * poly2[iTerm][jTerm]; 394 matrix->data.F64[iIndex][jIndex + numParams] = ab; 395 matrix->data.F64[jIndex + numParams][iIndex] = ab; 396 } 397 } 380 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 381 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 382 double ab = sumAB * poly2[iTerm][jTerm]; 383 matrix->data.F64[iIndex][jIndex + numParams] = ab; 384 matrix->data.F64[jIndex + numParams][iIndex] = ab; 385 } 398 386 } 399 387 } … … 402 390 double sumBI2 = 0.0; // Sum of B.I_2 products (for vector) 403 391 double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix, normalisation) 404 double sumA = 0.0; // Sum of A (for matrix, background)405 392 double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix, normalisation) 406 double sumB = 0.0; // Sum of B products (for matrix, background) 407 double sumI2 = 0.0; // Sum of I_2 (for vector, background) 393 394 double MxxAI1 = 0.0; 395 double MyyAI1 = 0.0; 396 double MxxBI2 = 0.0; 397 double MyyBI2 = 0.0; 408 398 for (int y = - footprint; y <= footprint; y++) { 409 399 for (int x = - footprint; x <= footprint; x++) { … … 413 403 float i2 = image2->kernel[y][x]; 414 404 415 double ai2 = a * i2 ;405 double ai2 = a * i2 * normValue; 416 406 double bi2 = b * i2; 417 double ai1 = a * i1 ;418 double bi1 = b * i1 ;407 double ai1 = a * i1 * PS_SQR(normValue); 408 double bi1 = b * i1 * normValue; 419 409 420 410 if (weight) { … … 424 414 ai1 *= wtVal; 425 415 bi1 *= wtVal; 426 a *= wtVal;427 b *= wtVal;428 i2 *= wtVal;429 416 } 430 417 if (window) { … … 434 421 ai1 *= wtVal; 435 422 bi1 *= wtVal; 436 a *= wtVal;437 b *= wtVal;438 i2 *= wtVal;439 423 } 440 424 sumAI2 += ai2; 441 425 sumBI2 += bi2; 442 426 sumAI1 += ai1; 443 sumA += a;444 427 sumBI1 += bi1; 445 sumB += b; 446 sumI2 += i2; 447 } 448 } 428 429 if (MOMENTS) { 430 MxxAI1 += x*x*ai1; 431 MyyAI1 += y*y*ai1; 432 MxxBI2 += x*x*bi2; 433 MyyBI2 += y*y*bi2; 434 } 435 } 436 } 437 438 if (MOMENTS) { 439 MxxAI1 /= stamp->normSquare1 * PS_SQR(normValue); 440 MyyAI1 /= stamp->normSquare1 * PS_SQR(normValue); 441 MxxBI2 /= stamp->normSquare2; 442 MyyBI2 /= stamp->normSquare2; 443 } 444 445 // fprintf (stderr, "i : %d : M(xx,yy)(AI1,BI2) : %f %f %f %f\n", i, MxxAI1, MyyAI1, MxxBI2, MyyBI2); 446 447 // XXX this makes the Chisq portion independent of the normalization and star flux 448 // but may be mis-scaling between stars of different fluxes 449 sumAI1 /= PS_SQR(stamp->normI2); 450 sumBI1 /= PS_SQR(stamp->normI2); 451 sumAI2 /= PS_SQR(stamp->normI2); 452 sumBI2 /= PS_SQR(stamp->normI2); 453 449 454 // Spatial variation 450 455 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { … … 452 457 double bi2 = sumBI2 * poly[iTerm]; 453 458 double ai1 = sumAI1 * poly[iTerm]; 454 double a = sumA * poly[iTerm];455 459 double bi1 = sumBI1 * poly[iTerm]; 456 double b = sumB * poly[iTerm]; 457 458 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 459 matrix->data.F64[iIndex][normIndex] = ai1; 460 matrix->data.F64[normIndex][iIndex] = ai1; 461 matrix->data.F64[iIndex + numParams][normIndex] = bi1; 462 matrix->data.F64[normIndex][iIndex + numParams] = bi1; 463 } 464 if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) { 465 matrix->data.F64[iIndex][bgIndex] = a; 466 matrix->data.F64[bgIndex][iIndex] = a; 467 matrix->data.F64[iIndex + numParams][bgIndex] = b; 468 matrix->data.F64[bgIndex][iIndex + numParams] = b; 469 } 470 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 471 vector->data.F64[iIndex] = ai2; 472 vector->data.F64[iIndex + numParams] = bi2; 473 if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) { 474 // subtract norm * sumRC * poly[iTerm] 475 psAssert (kernels->solution1, "programming error: define solution first!"); 476 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 477 double norm = fabs(kernels->solution1->data.F64[normIndex]); // Normalisation 478 vector->data.F64[iIndex] -= norm * ai1; 479 vector->data.F64[iIndex + numParams] -= norm * bi1; 480 } 481 } 482 } 483 } 484 485 double sumI1 = 0.0; // Sum of I_1 (for matrix, background-normalisation) 486 double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix, normalisation-normalisation) 487 double sum1 = 0.0; // Sum of 1 (for matrix, background-background) 488 double sumI2 = 0.0; // Sum of I_2 (for vector, background) 489 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector, normalisation) 490 double normI1 = 0.0, normI2 = 0.0; // Sum of I_1 and I_2 within the normalisation window 491 for (int y = - footprint; y <= footprint; y++) { 492 for (int x = - footprint; x <= footprint; x++) { 493 double i1 = image1->kernel[y][x]; 494 double i2 = image2->kernel[y][x]; 495 496 double i1i1 = i1 * i1; 497 double one = 1.0; 498 double i1i2 = i1 * i2; 499 500 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 501 normI1 += i1; 502 } 503 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 504 normI2 += i2; 505 } 506 507 if (weight) { 508 float wtVal = weight->kernel[y][x]; 509 i1 *= wtVal; 510 i1i1 *= wtVal; 511 one *= wtVal; 512 i2 *= wtVal; 513 i1i2 *= wtVal; 514 } 515 if (window) { 516 float wtVal = window->kernel[y][x]; 517 i1 *= wtVal; 518 i1i1 *= wtVal; 519 one *= wtVal; 520 i2 *= wtVal; 521 i1i2 *= wtVal; 522 } 523 sumI1 += i1; 524 sumI1I1 += i1i1; 525 sum1 += one; 526 sumI2 += i2; 527 sumI1I2 += i1i2; 528 } 529 } 530 531 *norm = normI2 / normI1; 532 fprintf (stderr, "normValue: %f %f %f\n", normI1, normI2, *norm); 533 534 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 535 matrix->data.F64[normIndex][normIndex] = sumI1I1; 536 vector->data.F64[normIndex] = sumI1I2; 537 } 538 if (mode & PM_SUBTRACTION_EQUATION_BG) { 539 matrix->data.F64[bgIndex][bgIndex] = sum1; 540 vector->data.F64[bgIndex] = sumI2; 541 } 542 if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) { 543 matrix->data.F64[bgIndex][normIndex] = sumI1; 544 matrix->data.F64[normIndex][bgIndex] = sumI1; 460 vector->data.F64[iIndex] = ai2 - ai1; 461 vector->data.F64[iIndex + numParams] = bi2 - bi1; 462 463 // fprintf (stderr, "i : %d : V(I1,I2) : %f %f\n", i, vector->data.F64[iIndex], vector->data.F64[iIndex + numParams]); 464 465 // add in second moments 466 if (MOMENTS) { 467 vector->data.F64[iIndex] -= kernels->penalty * MxxAI1; 468 vector->data.F64[iIndex] -= kernels->penalty * MyyAI1; 469 470 vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2; 471 vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2; 472 } 473 } 545 474 } 546 475 … … 565 494 } 566 495 567 #if 1568 496 // Add in penalty term to least-squares vector 569 497 bool calculatePenalty(psImage *matrix, // Matrix to which to add in penalty term 570 498 psVector *vector, // Vector to which to add in penalty term 571 499 const pmSubtractionKernels *kernels, // Kernel parameters 572 float norm // Normalisation 500 float normSquare1, // Normalisation for image 1 501 float normSquare2 // Normalisation for image 2 573 502 ) 574 503 { 504 psAssert (kernels->mode == PM_SUBTRACTION_MODE_DUAL, "only use penalties for dual convolution"); 505 575 506 if (kernels->penalty == 0.0) { 576 507 return true; … … 589 520 // [p_0,x_1,y_0 p_1,x_1,y_0, p_2,x_1,y_0] 590 521 // [p_0,x_0,y_1 p_1,x_0,y_1, p_2,x_0,y_1] 591 // [norm] 592 // [bg] 522 593 523 // [q_0,x_0,y_0 q_1,x_0,y_0, q_2,x_0,y_0] 594 524 // [q_0,x_1,y_0 q_1,x_1,y_0, q_2,x_1,y_0] … … 600 530 // Contribution to chi^2: a_i^2 P_i 601 531 psAssert(isfinite(penalties1->data.F32[i]), "Invalid penalty"); 602 fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], norm * penalties1->data.F32[i], norm, penalties1->data.F32[i]); 603 matrix->data.F64[index][index] += norm * penalties1->data.F32[i]; 604 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 605 fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 606 matrix->data.F64[index + numParams + 2][index + numParams + 2], norm * penalties2->data.F32[i], norm, penalties2->data.F32[i]); 607 matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties2->data.F32[i]; 608 // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) 609 // penalties scale with second moments 610 // 611 } 612 } 613 } 614 } 615 616 return true; 617 } 618 # endif 532 // fprintf (stderr, "penalty: %f + %f (%f * %f)\n", matrix->data.F64[index][index], normSquare1 * penalties1->data.F32[i], normSquare1, penalties1->data.F32[i]); 533 matrix->data.F64[index][index] += normSquare1 * penalties1->data.F32[i]; 534 535 // fprintf (stderr, "penalty: (x^%d y^%d fwhm %f) : %f + %f (%f * %f)\n", kernels->u->data.S32[index], kernels->v->data.S32[index], kernels->widths->data.F32[index], 536 // matrix->data.F64[index + numParams][index + numParams], normSquare2 * penalties2->data.F32[i], normSquare2, penalties2->data.F32[i]); 537 matrix->data.F64[index + numParams][index + numParams] += normSquare2 * penalties2->data.F32[i]; 538 // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1) 539 // penalties scale with second moments 540 // 541 } 542 } 543 } 544 545 return true; 546 } 619 547 620 548 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// … … 647 575 int index, bool wantDual) 648 576 { 649 #if 0650 577 // This is probably in a tight loop, so don't check inputs 651 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);652 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);653 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);654 PS_ASSERT_INT_POSITIVE(index, NAN);655 #endif656 657 578 psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector 658 579 return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, … … 662 583 double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) 663 584 { 664 #if 0665 585 // This is probably in a tight loop, so don't check inputs 666 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);667 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);668 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);669 #endif670 671 586 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 672 587 return kernels->solution1->data.F64[normIndex]; … … 676 591 const psImage *polyValues) 677 592 { 678 #if 0679 593 // This is probably in a tight loop, so don't check inputs 680 PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN);681 PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN);682 PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN);683 #endif684 685 594 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 686 595 return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); … … 690 599 // Public functions 691 600 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 601 602 bool pmSubtractionCalculateMoments( 603 pmSubtractionKernels *kernels, // Kernels 604 pmSubtractionStampList *stamps) 605 { 606 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 607 608 // XXX skip this, right? 609 return true; 610 611 // these are only used by DUAL mode 612 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) return true; 613 614 psTimerStart("pmSubtractionCalculateMoments"); 615 616 int footprint = stamps->footprint; // Half-size of stamps 617 618 // Loop over each stamp and calculate its normalization factor 619 for (int i = 0; i < stamps->num; i++) { 620 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 621 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; 622 if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; 623 624 pmSubtractionCalculateMomentsStamp(kernels, stamp, footprint, stamps->normWindow2, stamps->normWindow1); 625 } 626 627 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate moments: %f sec", psTimerClear("pmSubtractionCalculateMoments")); 628 629 return true; 630 } 631 632 bool pmSubtractionCalculateMomentsStamp( 633 pmSubtractionKernels *kernels, // Kernels 634 pmSubtractionStamp *stamp, // stamp on which to save normalization) 635 int footprint, // (Half-)Size of stamp 636 int normWindow1, // Window (half-)size for normalisation measurement 637 int normWindow2 // Window (half-)size for normalisation measurement 638 ) 639 { 640 double Mxx, Myy; 641 642 int numKernels = kernels->num; 643 644 // Generate convolutions: these are generated once and saved 645 if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { 646 psError(psErrorCodeLast(), false, "Unable to convolve stamp"); 647 return false; 648 } 649 650 if (!stamp->MxxI1) { 651 stamp->MxxI1 = psVectorAlloc (numKernels, PS_TYPE_F32); 652 } 653 if (!stamp->MyyI1) { 654 stamp->MyyI1 = psVectorAlloc (numKernels, PS_TYPE_F32); 655 } 656 if (!stamp->MxxI2) { 657 stamp->MxxI2 = psVectorAlloc (numKernels, PS_TYPE_F32); 658 } 659 if (!stamp->MyyI2) { 660 stamp->MyyI2 = psVectorAlloc (numKernels, PS_TYPE_F32); 661 } 662 663 for (int i = 0; i < numKernels; i++) { 664 pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions1->data[i], footprint, normWindow1); 665 stamp->MxxI1->data.F32[i] = Mxx / stamp->normI1; 666 stamp->MyyI1->data.F32[i] = Myy / stamp->normI1; 667 pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->convolutions2->data[i], footprint, normWindow2); 668 stamp->MxxI2->data.F32[i] = Mxx / stamp->normI2; 669 stamp->MyyI2->data.F32[i] = Myy / stamp->normI2; 670 } 671 672 pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image1, footprint, normWindow1); 673 stamp->MxxI1raw = Mxx / stamp->normI1; 674 stamp->MyyI1raw = Myy / stamp->normI1; 675 676 pmSubtractionCalculateMomentsKernel(&Mxx, &Myy, stamp->image2, footprint, normWindow2); 677 stamp->MxxI2raw = Mxx / stamp->normI2; 678 stamp->MyyI2raw = Myy / stamp->normI2; 679 680 // fprintf (stderr, "Mxx I1: %f, Myy I1: %f, Mxx I2: %f, Myy I2: %f\n", stamp->MxxI1raw, stamp->MyyI1raw, stamp->MxxI2raw, stamp->MyyI2raw); 681 682 return true; 683 } 684 685 bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window) { 686 687 double Sxx = 0.0; 688 double Syy = 0.0; 689 for (int y = - footprint; y <= footprint; y++) { 690 for (int x = - footprint; x <= footprint; x++) { 691 if (PS_SQR(x) + PS_SQR(y) > PS_SQR(window)) continue; 692 693 double flux = image->kernel[y][x]; 694 695 Sxx += PS_SQR(x) * flux; 696 Syy += PS_SQR(y) * flux; 697 } 698 } 699 *Mxx = Sxx; 700 *Myy = Syy; 701 return true; 702 } 703 704 ///--------- 705 706 bool pmSubtractionCalculateNormalization( 707 pmSubtractionStampList *stamps, 708 const pmSubtractionMode mode) 709 { 710 PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); 711 712 psTimerStart("pmSubtractionCalculateNormalization"); 713 714 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 715 psVector *norm2 = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 716 717 int footprint = stamps->footprint; // Half-size of stamps 718 719 // Loop over each stamp and calculate its normalization factor 720 for (int i = 0; i < stamps->num; i++) { 721 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 722 if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) continue; 723 if (stamp->status == PM_SUBTRACTION_STAMP_NONE) continue; 724 725 // XXX skip this if we have already calculated it? (stamp->norm does not change, just the median statistic) 726 // XXX maybe not: the star may have changed for a given stamp -- only if norm is reset to NAN can we do this 727 if (mode == PM_SUBTRACTION_MODE_2) { 728 pmSubtractionCalculateNormalizationStamp(stamp, stamp->image2, stamp->image1, footprint, stamps->normWindow2, stamps->normWindow1); 729 } else { 730 pmSubtractionCalculateNormalizationStamp(stamp, stamp->image1, stamp->image2, footprint, stamps->normWindow1, stamps->normWindow2); 731 } 732 psVectorAppend(norms, stamp->norm); 733 psVectorAppend(norm2, stamp->normSquare2 / stamp->normSquare1); 734 } 735 736 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 737 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 738 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 739 psFree(stats); 740 psFree(norms); 741 psFree(norm2); 742 return false; 743 } 744 stamps->normValue = stats->robustMedian; 745 746 psStatsInit(stats); 747 if (!psVectorStats(stats, norm2, NULL, NULL, 0)) { 748 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 749 psFree(stats); 750 psFree(norms); 751 psFree(norm2); 752 return false; 753 } 754 stamps->normValue2 = stats->robustMedian; 755 756 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2); 757 758 psFree(stats); 759 psFree(norms); 760 psFree(norm2); 761 762 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate normalization: %f sec", psTimerClear("pmSubtractionCalculateNormalization")); 763 764 return true; 765 } 766 767 bool pmSubtractionCalculateNormalizationStamp( 768 pmSubtractionStamp *stamp, // stamp on which to save normalization) 769 const psKernel *image1, // Input image (target) 770 const psKernel *image2, // Reference image (convolution source) 771 int footprint, // (Half-)Size of stamp 772 int normWindow1, // Window (half-)size for normalisation measurement 773 int normWindow2 // Window (half-)size for normalisation measurement 774 ) 775 { 776 double normI1 = 0.0; // Sum of I_1 within the normalisation window (aperture) 777 double normI2 = 0.0; // Sum of I_2 within the normalisation window (aperture) 778 double normSquare1 = 0.0; // Sum of (I_1)^2 within the normalisation window (aperture) 779 double normSquare2 = 0.0; // Sum of (I_2)^2 within the normalisation window (aperture) 780 for (int y = - footprint; y <= footprint; y++) { 781 for (int x = - footprint; x <= footprint; x++) { 782 double im1 = image1->kernel[y][x]; 783 double im2 = image2->kernel[y][x]; 784 785 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow1)) { 786 normI1 += im1; 787 normSquare1 += PS_SQR(im1); 788 } 789 if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow2)) { 790 normI2 += im2; 791 normSquare2 += PS_SQR(im2); 792 } 793 } 794 } 795 stamp->norm = normI2 / normI1; 796 stamp->normI1 = normI1; 797 stamp->normI2 = normI2; 798 stamp->normSquare1 = normSquare1; 799 stamp->normSquare2 = normSquare2; 800 801 psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); 802 803 return true; 804 } 692 805 693 806 bool pmSubtractionCalculateEquationThread(psThreadJob *job) … … 715 828 int numKernels = kernels->num; // Number of kernel basis functions 716 829 int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations 717 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms718 830 719 831 // numKernels is the number of unique kernel images (one for each Gaussian modified by a specific polynomial). … … 722 834 // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the 723 835 // number of coefficients for the spatial polynomial, normalisation and a constant background offset. 724 int numParams = numKernels * numSpatial + 1 + numBackground; 836 // XXX we no longer solve for the normalization and background in the matrix inversion 837 int numParams = numKernels * numSpatial; 725 838 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 726 839 // An additional image is convolved … … 790 903 switch (kernels->mode) { 791 904 case PM_SUBTRACTION_MODE_1: 792 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image2, stamp->image1, 793 weight, window, stamp->convolutions1, kernels, 794 polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); 905 status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image2, stamp->image1, 906 weight, window, stamp->convolutions1, kernels, polyValues, footprint); 795 907 break; 796 908 case PM_SUBTRACTION_MODE_2: 797 status = calculateMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, stamp->image1, stamp->image2, 798 weight, window, stamp->convolutions2, kernels, 799 polyValues, footprint, stamps->normWindow2, stamps->normWindow1, mode); 909 status = calculateMatrixVector(stamp->matrix, stamp->vector, stamps->normValue, stamp->image1, stamp->image2, 910 weight, window, stamp->convolutions2, kernels, polyValues, footprint); 800 911 break; 801 912 case PM_SUBTRACTION_MODE_DUAL: 802 status = calculateDualMatrixVector(stamp->matrix, stamp->vector, &stamp->norm, 803 stamp->image1, stamp->image2, 804 weight, window, stamp->convolutions1, stamp->convolutions2, 805 kernels, polyValues, footprint, stamps->normWindow1, stamps->normWindow2, mode); 913 status = calculateDualMatrixVector(stamp, stamps->normValue, stamps->normValue2, weight, window, kernels, polyValues, footprint); 806 914 break; 807 915 default: … … 928 1036 int numKernels = kernels->num; // Number of kernel basis functions 929 1037 int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations 930 int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 931 int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for 932 int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution 1038 // XXX int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms 1039 // XXX int numParams = numKernels * numSpatial + 1 + numBackground; // Number of parameters being solved for 1040 int numParams = numKernels * numSpatial; // Number of parameters being solved for 1041 int numSolution1 = numParams, numSolution2 = 0; // Number of parameters for each solution 933 1042 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { 934 1043 // An additional image is convolved … … 960 1069 961 1070 if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { 962 // Accumulate the least-squares matricies and vectors 1071 // Accumulate the least-squares matrices and vectors. These are generated for the 1072 // kernel elements, excluding the background and normalization. 963 1073 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix 964 1074 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector … … 966 1076 psImageInit(sumMatrix, 0.0); 967 1077 968 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations969 970 1078 int numStamps = 0; // Number of good stamps 1079 for (int i = 0; i < stamps->num; i++) { 1080 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 1081 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 1082 1083 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); 1084 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 1085 1086 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 1087 numStamps++; 1088 } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) { 1089 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); 1090 } 1091 } 1092 1093 #if 0 1094 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); 1095 psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); 1096 psVectorWriteFile("sumVector.dat", sumVector); 1097 psFree (save); 1098 #endif 1099 1100 psVector *solution = NULL; // Solution to equation! 1101 solution = psVectorAlloc(numParams, PS_TYPE_F64); 1102 psVectorInit(solution, 0); 1103 1104 // XXX TEST: try some constraint on the svd solution 1105 // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1106 // SINGLE solution 1107 if (1) { 1108 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1109 } else { 1110 psVector *PERM = NULL; 1111 psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix); 1112 solution = psMatrixLUSolution(solution, LU, sumVector, PERM); 1113 psFree (LU); 1114 psFree (PERM); 1115 } 1116 1117 # if (0) 1118 for (int i = 0; i < solution->n; i++) { 1119 psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]); 1120 } 1121 # endif 1122 1123 if (!kernels->solution1) { 1124 kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg 1125 psVectorInit(kernels->solution1, 0.0); 1126 } 1127 1128 int numKernels = kernels->num; 1129 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 1130 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 1131 for (int i = 0; i < numKernels * numPoly; i++) { 1132 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1133 } 1134 1135 // Apply the normalisation and background separately 1136 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1137 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1138 kernels->solution1->data.F64[normIndex] = stamps->normValue; 1139 kernels->solution1->data.F64[bgIndex] = 0.0; 1140 1141 psFree(solution); 1142 psFree(sumVector); 1143 psFree(sumMatrix); 1144 1145 } else { 1146 // Dual convolution solution 1147 // Accumulate the least-squares matrices and vectors. These are generated for the 1148 // kernel elements, excluding the background and normalization. 1149 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); 1150 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); 1151 psImageInit(sumMatrix, 0.0); 1152 psVectorInit(sumVector, 0.0); 1153 1154 int numStamps = 0; // Number of good stamps 1155 double normSquare1 = 0.0; // Sum of (I_1)^2 over stamps 1156 double normSquare2 = 0.0; // Sum of (I_2)^2 over stamps 971 1157 for (int i = 0; i < stamps->num; i++) { 972 1158 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest … … 974 1160 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); 975 1161 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 976 psVectorAppend(norms, stamp->norm); 1162 1163 normSquare1 += stamp->normSquare1; 1164 normSquare2 += stamp->normSquare2; 1165 977 1166 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 978 1167 numStamps++; … … 983 1172 984 1173 #if 0 985 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 986 calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); 1174 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); 1175 psFitsWriteImageSimple ("sumMatrix.fits", save, NULL); 1176 psVectorWriteFile("sumVector.dat", sumVector); 1177 psFree (save); 987 1178 #endif 1179 1180 if (PENALTY) { 1181 calculatePenalty(sumMatrix, sumVector, kernels, normSquare1, normSquare2); 1182 } 988 1183 989 1184 psVector *solution = NULL; // Solution to equation! … … 991 1186 psVectorInit(solution, 0); 992 1187 993 #if 0 994 // Regular, straight-forward solution 995 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 996 #else 997 { 998 // Solve normalisation and background separately 999 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1000 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1001 1002 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 1003 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 1004 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 1005 psFree(stats); 1006 psFree(sumMatrix); 1007 psFree(sumVector); 1008 psFree(norms); 1009 return false; 1010 } 1011 1012 // double normValue = 1.0; 1013 double normValue = stats->robustMedian; 1014 // double bgValue = 0.0; 1015 1016 psFree(stats); 1017 1018 #ifdef TESTING 1019 fprintf(stderr, "Norm: %lf\n", normValue); 1020 #endif 1021 // Solve kernel components 1022 for (int i = 0; i < numSolution1; i++) { 1023 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1024 1025 sumMatrix->data.F64[i][normIndex] = 0.0; 1026 sumMatrix->data.F64[normIndex][i] = 0.0; 1027 } 1028 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1029 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1030 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1031 1032 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1033 sumVector->data.F64[normIndex] = 0.0; 1034 1035 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1036 1037 solution->data.F64[normIndex] = normValue; 1038 } 1039 # endif 1040 1041 #if (1) 1042 for (int i = 0; i < solution->n; i++) { 1043 fprintf(stderr, "Single solution %d: %lf\n", i, solution->data.F64[i]); 1044 } 1045 #endif 1046 1047 if (!kernels->solution1) { 1048 kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64); 1049 psVectorInit(kernels->solution1, 0.0); 1050 } 1051 1052 // only update the solutions that we chose to calculate: 1053 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1054 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1055 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1056 } 1057 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1058 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1059 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1060 } 1061 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1062 int numKernels = kernels->num; 1063 int spatialOrder = kernels->spatialOrder; // Order of spatial variation 1064 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms 1065 for (int i = 0; i < numKernels * numPoly; i++) { 1066 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1067 } 1068 } 1069 1070 psFree(norms); 1071 psFree(solution); 1072 psFree(sumVector); 1073 psFree(sumMatrix); 1074 1075 #ifdef TESTING 1076 // XXX double-check for NAN in data: 1077 for (int ix = 0; ix < kernels->solution1->n; ix++) { 1078 if (!isfinite(kernels->solution1->data.F64[ix])) { 1079 fprintf (stderr, "WARNING: NAN in vector\n"); 1080 } 1081 } 1082 #endif 1083 1084 } else { 1085 // Dual convolution solution 1086 1087 // Accumulation of stamp matrices/vectors 1088 psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); 1089 psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); 1090 psImageInit(sumMatrix, 0.0); 1091 psVectorInit(sumVector, 0.0); 1092 1093 psVector *norms = psVectorAllocEmpty(stamps->num, PS_TYPE_F64); // Normalisations 1094 1095 int numStamps = 0; // Number of good stamps 1096 for (int i = 0; i < stamps->num; i++) { 1097 pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest 1098 if (stamp->status == PM_SUBTRACTION_STAMP_USED) { 1099 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); 1100 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); 1101 1102 psVectorAppend(norms, stamp->norm); 1103 1104 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); 1105 numStamps++; 1106 } 1107 } 1108 1109 #ifdef TESTING 1110 psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL); 1111 psVectorWriteFile("sumVector.dat", sumVector); 1112 #endif 1113 1114 #if 1 1115 // int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1116 // calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]); 1117 1118 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1119 calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[normIndex][normIndex] / 100.0); 1120 #endif 1121 1122 psVector *solution = NULL; // Solution to equation! 1123 solution = psVectorAlloc(numParams, PS_TYPE_F64); 1124 psVectorInit(solution, 0); 1125 1126 #if 0 1127 // Regular, straight-forward solution 1128 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1129 #else 1130 { 1131 // Solve normalisation and background separately 1132 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1133 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1134 1135 #if 0 1136 psImage *normMatrix = psImageAlloc(2, 2, PS_TYPE_F64); 1137 psVector *normVector = psVectorAlloc(2, PS_TYPE_F64); 1138 1139 normMatrix->data.F64[0][0] = sumMatrix->data.F64[normIndex][normIndex]; 1140 normMatrix->data.F64[1][1] = sumMatrix->data.F64[bgIndex][bgIndex]; 1141 normMatrix->data.F64[0][1] = normMatrix->data.F64[1][0] = sumMatrix->data.F64[normIndex][bgIndex]; 1142 1143 normVector->data.F64[0] = sumVector->data.F64[normIndex]; 1144 normVector->data.F64[1] = sumVector->data.F64[bgIndex]; 1145 1146 psVector *normSolution = psMatrixSolveSVD(NULL, normMatrix, normVector, NAN); 1147 1148 double normValue = normSolution->data.F64[0]; 1149 double bgValue = normSolution->data.F64[1]; 1150 1151 psFree(normMatrix); 1152 psFree(normVector); 1153 psFree(normSolution); 1154 #endif 1155 1156 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm 1157 if (!psVectorStats(stats, norms, NULL, NULL, 0)) { 1158 psError(PM_ERR_DATA, false, "Unable to determine median normalisation"); 1159 psFree(stats); 1160 psFree(sumMatrix); 1161 psFree(sumVector); 1162 psFree(norms); 1163 return false; 1164 } 1165 1166 double normValue = stats->robustMedian; 1167 1168 psFree(stats); 1169 1170 #ifdef TESTING 1171 fprintf(stderr, "Norm: %lf\n", normValue); 1172 #endif 1173 1174 // Solve kernel components 1175 for (int i = 0; i < numSolution2; i++) { 1176 sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; 1177 sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; 1178 1179 sumMatrix->data.F64[i][normIndex] = 0.0; 1180 sumMatrix->data.F64[normIndex][i] = 0.0; 1181 1182 sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0; 1183 sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0; 1184 } 1185 sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex]; 1186 sumMatrix->data.F64[bgIndex][normIndex] = 0.0; 1187 sumMatrix->data.F64[normIndex][bgIndex] = 0.0; 1188 1189 sumMatrix->data.F64[normIndex][normIndex] = 1.0; 1190 1191 sumVector->data.F64[normIndex] = 0.0; 1192 1193 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1194 1195 solution->data.F64[normIndex] = normValue; 1196 } 1197 #endif 1198 1199 1200 #if (1) 1188 // DUAL solution 1189 solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); 1190 1191 #if (0) 1201 1192 for (int i = 0; i < solution->n; i++) { 1202 1193 fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); … … 1204 1195 #endif 1205 1196 1206 psFree(sumMatrix);1207 psFree(sumVector);1208 1209 psFree(norms);1210 1211 1197 if (!kernels->solution1) { 1212 kernels->solution1 = psVectorAlloc(numSolution1 , PS_TYPE_F64);1198 kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); 1213 1199 psVectorInit (kernels->solution1, 0.0); 1214 1200 } … … 1218 1204 } 1219 1205 1220 // only update the solutions that we chose to calculate: 1221 if (mode & PM_SUBTRACTION_EQUATION_NORM) { 1222 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1223 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex]; 1224 } 1225 if (mode & PM_SUBTRACTION_EQUATION_BG) { 1226 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background 1227 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex]; 1228 } 1229 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) { 1230 int numKernels = kernels->num; 1231 for (int i = 0; i < numKernels * numSpatial; i++) { 1232 // XXX fprintf (stderr, "keep\n"); 1233 kernels->solution1->data.F64[i] = solution->data.F64[i]; 1234 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; 1235 } 1236 } 1237 1238 1239 memcpy(kernels->solution1->data.F64, solution->data.F64, 1240 numSolution1 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1241 memcpy(kernels->solution2->data.F64, &solution->data.F64[numSolution1], 1242 numSolution2 * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); 1243 1206 // for DUAL convolution analysis, we apply the normalization to I1 as follows: 1207 // I1c = norm I1 + \Sum_i a_i norm I1 \cross k_i 1208 // I2c = I2 + \Sum_i b_i I2 \cross k_i 1209 1210 // We absorb the normalization into a_i after the analysis is complete to be consistent 1211 // with the SINGLE definitions of the convolutions 1212 1213 int numKernels = kernels->num; 1214 for (int i = 0; i < numKernels * numSpatial; i++) { 1215 // we solve for coefficients 1216 kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue; 1217 kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; 1218 } 1219 1220 // Apply the normalisation and background separately 1221 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation 1222 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background 1223 kernels->solution1->data.F64[normIndex] = stamps->normValue; 1224 kernels->solution1->data.F64[bgIndex] = 0.0; 1225 1226 psFree(sumMatrix); 1227 psFree(sumVector); 1244 1228 psFree(solution); 1245 1246 1229 } 1247 1230 … … 1265 1248 } 1266 1249 1267 bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { 1268 1269 // XXX measure some useful stats on the residuals 1250 // measure some useful stats on the stamp residuals: 1251 // fResSigma : the residual stdev / total flux 1252 // fResOuter : the residual fabs / total flux for R > 2 pix 1253 // fResTotal : the residual fabs / total flux for R > 0 pix 1254 bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) { 1255 1270 1256 float sum = 0.0; 1271 1257 float peak = 0.0; … … 1277 1263 } 1278 1264 1279 // only count pixels with more than X% of the source flux1280 // calculate stdev(dflux)1265 // init counters 1266 int npix = 0; 1281 1267 float dflux1 = 0.0; 1282 1268 float dflux2 = 0.0; 1283 int npix = 0; 1284 1285 float dmax = 0.0; 1286 float dmin = 0.0; 1269 float dOuter = 0.0; 1270 float dTotal = 0.0; 1287 1271 1288 1272 for (int y = - footprint; y <= footprint; y++) { 1289 1273 for (int x = - footprint; x <= footprint; x++) { 1290 float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);1291 if (dflux < 0.02*sum) continue;1292 1274 dflux1 += residual->kernel[y][x]; 1293 1275 dflux2 += PS_SQR(residual->kernel[y][x]); 1294 dmax = PS_MAX(residual->kernel[y][x], dmax); 1295 dmin = PS_MIN(residual->kernel[y][x], dmin); 1276 dTotal += fabs(residual->kernel[y][x]); 1277 if (hypot(x,y) > 2.0) { 1278 dOuter += fabs(residual->kernel[y][x]); 1279 } 1296 1280 npix ++; 1297 1281 } … … 1299 1283 float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix)); 1300 1284 if (!isfinite(sum)) return false; 1301 if (!isfinite(dmax)) return false;1302 if (!isfinite(dmin)) return false;1303 1285 if (!isfinite(peak)) return false; 1304 1305 // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/peak, dmin/peak); 1306 psVectorAppend(fSigRes, sigma/sum); 1307 psVectorAppend(fMaxRes, dmax/peak); 1308 psVectorAppend(fMinRes, dmin/peak); 1286 if (!isfinite(dOuter)) return false; 1287 if (!isfinite(dTotal)) return false; 1288 1289 // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum); 1290 psVectorAppend(fResSigma, sigma/sum); 1291 psVectorAppend(fResOuter, dOuter/sum); 1292 psVectorAppend(fResTotal, dTotal/sum); 1309 1293 return true; 1310 1294 } … … 1329 1313 pmSubtractionVisualShowFitInit (stamps); 1330 1314 1331 psVector *f SigRes= psVectorAllocEmpty(stamps->num, PS_TYPE_F32);1332 psVector *f MinRes= psVectorAllocEmpty(stamps->num, PS_TYPE_F32);1333 psVector *f MaxRes= psVectorAllocEmpty(stamps->num, PS_TYPE_F32);1315 psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1316 psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1317 psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); 1334 1318 1335 1319 // we want to save the residual images for the 9 brightest stamps. … … 1443 1427 } 1444 1428 1445 pmSubtractionResidualStats(f SigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);1429 pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint); 1446 1430 1447 1431 } else { … … 1479 1463 } 1480 1464 1481 pmSubtractionResidualStats(f SigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);1465 pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint); 1482 1466 } 1483 1467 … … 1558 1542 1559 1543 psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); 1560 psVectorStats (stats, f SigRes, NULL, NULL, 0);1561 kernels->f SigResMean= stats->robustMedian;1562 kernels->f SigResStdev = stats->robustStdev;1544 psVectorStats (stats, fResSigma, NULL, NULL, 0); 1545 kernels->fResSigmaMean = stats->robustMedian; 1546 kernels->fResSigmaStdev = stats->robustStdev; 1563 1547 1564 1548 psStatsInit (stats); 1565 psVectorStats (stats, f MaxRes, NULL, NULL, 0);1566 kernels->f MaxResMean= stats->robustMedian;1567 kernels->f MaxResStdev = stats->robustStdev;1549 psVectorStats (stats, fResOuter, NULL, NULL, 0); 1550 kernels->fResOuterMean = stats->robustMedian; 1551 kernels->fResOuterStdev = stats->robustStdev; 1568 1552 1569 1553 psStatsInit (stats); 1570 psVectorStats (stats, f MinRes, NULL, NULL, 0);1571 kernels->f MinResMean= stats->robustMedian;1572 kernels->f MinResStdev = stats->robustStdev;1554 psVectorStats (stats, fResTotal, NULL, NULL, 0); 1555 kernels->fResTotalMean = stats->robustMedian; 1556 kernels->fResTotalStdev = stats->robustStdev; 1573 1557 1574 1558 // XXX save these values somewhere 1575 psLogMsg("psModules.imcombine", PS_LOG_INFO, "f Sigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",1576 kernels->f SigResMean, kernels->fSigResStdev,1577 kernels->f MaxResMean, kernels->fMaxResStdev,1578 kernels->f MinResMean, kernels->fMinResStdev);1579 1580 psFree (f SigRes);1581 psFree (f MaxRes);1582 psFree (f MinRes);1559 psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f", 1560 kernels->fResSigmaMean, kernels->fResSigmaStdev, 1561 kernels->fResOuterMean, kernels->fResOuterStdev, 1562 kernels->fResTotalMean, kernels->fResTotalStdev); 1563 1564 psFree (fResSigma); 1565 psFree (fResOuter); 1566 psFree (fResTotal); 1583 1567 psFree (stats); 1584 1568 } … … 1586 1570 psFree(residual); 1587 1571 psFree(polyValues); 1588 1589 1572 1590 1573 return deviations; … … 1912 1895 return true; 1913 1896 } 1914 1915 1916 # if 01917 1918 #ifdef TESTING1919 psFitsWriteImageSimple("A.fits", sumMatrix, NULL);1920 psVectorWriteFile ("B.dat", sumVector);1921 #endif1922 1923 # define SVD_ANALYSIS 01924 # define COEFF_SIG 0.01925 # define SVD_TOL 0.01926 1927 // Use SVD to determine the kernel coeffs (and validate)1928 if (SVD_ANALYSIS) {1929 1930 // We have sumVector and sumMatrix. we are trying to solve the following equation:1931 // sumMatrix * x = sumVector.1932 1933 // we can use any standard matrix inversion to solve this. However, the basis1934 // functions in general have substantial correlation, so that the solution may be1935 // somewhat poorly determined or unstable. If not numerically ill-conditioned, the1936 // system of equations may be statistically ill-conditioned. Noise in the image1937 // will drive insignificant, but correlated, terms in the solution. To avoid these1938 // problems, we can use SVD to identify numerically unconstrained values and to1939 // avoid statistically badly determined value.1940 1941 // A = sumMatrix, B = sumVector1942 // SVD: A = U w V^T -> A^{-1} = V (1/w) U^T1943 // x = V (1/w) (U^T B)1944 // \sigma_x = sqrt(diag(A^{-1}))1945 // solve for x and A^{-1} to get x & dx1946 // identify the elements of (1/w) that are nan (1/0.0) -> set to 0.01947 // identify the elements of x that are insignificant (x / dx < 1.0? < 0.5?) -> set to 0.01948 1949 // If I use the SVD trick to re-condition the matrix, I need to break out the1950 // kernel and normalization terms from the background term.1951 // XXX is this true? or was this due to an error in the analysis?1952 1953 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background1954 1955 // now pull out the kernel elements into their own square matrix1956 psImage *kernelMatrix = psImageAlloc (sumMatrix->numCols - 1, sumMatrix->numRows - 1, PS_TYPE_F64);1957 psVector *kernelVector = psVectorAlloc (sumMatrix->numCols - 1, PS_TYPE_F64);1958 1959 for (int ix = 0, kx = 0; ix < sumMatrix->numCols; ix++) {1960 if (ix == bgIndex) continue;1961 for (int iy = 0, ky = 0; iy < sumMatrix->numRows; iy++) {1962 if (iy == bgIndex) continue;1963 kernelMatrix->data.F64[ky][kx] = sumMatrix->data.F64[iy][ix];1964 ky++;1965 }1966 kernelVector->data.F64[kx] = sumVector->data.F64[ix];1967 kx++;1968 }1969 1970 psImage *U = NULL;1971 psImage *V = NULL;1972 psVector *w = NULL;1973 if (!psMatrixSVD (&U, &w, &V, kernelMatrix)) {1974 psError(psErrorCodeLast(), false, "failed to perform SVD on sumMatrix\n");1975 return NULL;1976 }1977 1978 // calculate A_inverse:1979 // Ainv = V * w * U^T1980 psImage *wUt = p_pmSubSolve_wUt (w, U);1981 psImage *Ainv = p_pmSubSolve_VwUt (V, wUt);1982 psImage *Xvar = NULL;1983 psFree (wUt);1984 1985 # ifdef TESTING1986 // kernel terms:1987 for (int i = 0; i < w->n; i++) {1988 fprintf (stderr, "w: %f\n", w->data.F64[i]);1989 }1990 # endif1991 // loop over w adding in more and more of the values until chisquare is no longer1992 // dropping significantly.1993 // XXX this does not seem to work very well: we seem to need all terms even for1994 // simple cases...1995 1996 psVector *Xsvd = NULL;1997 {1998 psVector *Ax = NULL;1999 psVector *UtB = NULL;2000 psVector *wUtB = NULL;2001 2002 psVector *wApply = psVectorAlloc(w->n, PS_TYPE_F64);2003 psVector *wMask = psVectorAlloc(w->n, PS_TYPE_U8);2004 psVectorInit (wMask, 1); // start by masking everything2005 2006 double chiSquareLast = NAN;2007 int maxWeight = 0;2008 2009 double Axx, Bx, y2;2010 2011 // XXX this is an attempt to exclude insignificant modes.2012 // it was not successful with the ISIS kernel set: removing even2013 // the least significant mode leaves additional ringing / noise2014 // because the terms are so coupled.2015 for (int k = 0; false && (k < w->n); k++) {2016 2017 // unmask the k-th weight2018 wMask->data.U8[k] = 0;2019 p_pmSubSolve_SetWeights(wApply, w, wMask);2020 2021 // solve for x:2022 // x = V * w * (U^T * B)2023 p_pmSubSolve_UtB (&UtB, U, kernelVector);2024 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);2025 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);2026 2027 // chi-square for this system of equations:2028 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^22029 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^22030 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);2031 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);2032 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);2033 p_pmSubSolve_y2 (&y2, kernels, stamps);2034 2035 // apparently, this works (compare with the brute force value below2036 double chiSquare = Axx - 2.0*Bx + y2;2037 double deltaChi = (k == 0) ? chiSquare : chiSquareLast - chiSquare;2038 chiSquareLast = chiSquare;2039 2040 // fprintf (stderr, "chi square = %f, delta: %f\n", chiSquare, deltaChi);2041 if (k && !maxWeight && (deltaChi < 1.0)) {2042 maxWeight = k;2043 }2044 }2045 2046 // keep all terms or we get extra ringing2047 maxWeight = w->n;2048 psVectorInit (wMask, 1);2049 for (int k = 0; k < maxWeight; k++) {2050 wMask->data.U8[k] = 0;2051 }2052 p_pmSubSolve_SetWeights(wApply, w, wMask);2053 2054 // solve for x:2055 // x = V * w * (U^T * B)2056 p_pmSubSolve_UtB (&UtB, U, kernelVector);2057 p_pmSubSolve_wUtB (&wUtB, wApply, UtB);2058 p_pmSubSolve_VwUtB (&Xsvd, V, wUtB);2059 2060 // chi-square for this system of equations:2061 // chi-square = sum over terms of: (Ax - B)*x - b*x - y^22062 // y^2 = \sum_stamps \sum_pixels input->kernel[y][x]^22063 p_pmSubSolve_Ax (&Ax, kernelMatrix, Xsvd);2064 p_pmSubSolve_VdV (&Axx, Ax, Xsvd);2065 p_pmSubSolve_VdV (&Bx, kernelVector, Xsvd);2066 p_pmSubSolve_y2 (&y2, kernels, stamps);2067 2068 // apparently, this works (compare with the brute force value below2069 double chiSquare = Axx - 2.0*Bx + y2;2070 psLogMsg ("psModules.imcombine", PS_LOG_INFO, "model kernel with %d terms; chi square = %f\n", maxWeight, chiSquare);2071 2072 // re-calculate A^{-1} to get new variances:2073 // Ainv = V * w * U^T2074 // XXX since we keep all terms, this is identical to Ainv2075 psImage *wUt = p_pmSubSolve_wUt (wApply, U);2076 Xvar = p_pmSubSolve_VwUt (V, wUt);2077 psFree (wUt);2078 2079 psFree (Ax);2080 psFree (UtB);2081 psFree (wUtB);2082 psFree (wApply);2083 psFree (wMask);2084 }2085 2086 // copy the kernel solutions to the full solution vector:2087 solution = psVectorAlloc(sumVector->n, PS_TYPE_F64);2088 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);2089 2090 for (int ix = 0, kx = 0; ix < sumVector->n; ix++) {2091 if (ix == bgIndex) {2092 solution->data.F64[ix] = 0;2093 solutionErr->data.F64[ix] = 0.001;2094 continue;2095 }2096 solutionErr->data.F64[ix] = sqrt(Ainv->data.F64[kx][kx]);2097 solution->data.F64[ix] = Xsvd->data.F64[kx];2098 kx++;2099 }2100 2101 psFree (kernelMatrix);2102 psFree (kernelVector);2103 2104 psFree (U);2105 psFree (V);2106 psFree (w);2107 2108 psFree (Ainv);2109 psFree (Xsvd);2110 } else {2111 psVector *permutation = NULL; // Permutation vector, required for LU decomposition2112 psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);2113 if (!luMatrix) {2114 psError(PM_ERR_DATA, true, "LU Decomposition of least-squares matrix failed.\n");2115 psFree(solution);2116 psFree(sumVector);2117 psFree(sumMatrix);2118 psFree(luMatrix);2119 psFree(permutation);2120 return NULL;2121 }2122 2123 solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);2124 psFree(luMatrix);2125 psFree(permutation);2126 if (!solution) {2127 psError(PM_ERR_DATA, true, "Failed to solve the least-squares system.\n");2128 psFree(solution);2129 psFree(sumVector);2130 psFree(sumMatrix);2131 return NULL;2132 }2133 2134 // XXX LUD does not provide A^{-1}? fake the error for now2135 solutionErr = psVectorAlloc(sumVector->n, PS_TYPE_F64);2136 for (int ix = 0; ix < sumVector->n; ix++) {2137 solutionErr->data.F64[ix] = 0.1*solution->data.F64[ix];2138 }2139 }2140 2141 if (!kernels->solution1) {2142 kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);2143 psVectorInit (kernels->solution1, 0.0);2144 }2145 2146 // only update the solutions that we chose to calculate:2147 if (mode & PM_SUBTRACTION_EQUATION_NORM) {2148 int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation2149 kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];2150 }2151 if (mode & PM_SUBTRACTION_EQUATION_BG) {2152 int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background2153 kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];2154 }2155 if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {2156 int numKernels = kernels->num;2157 int spatialOrder = kernels->spatialOrder; // Order of spatial variation2158 int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms2159 for (int i = 0; i < numKernels * numPoly; i++) {2160 // XXX fprintf (stderr, "%f +/- %f (%f) -> ", solution->data.F64[i], solutionErr->data.F64[i], fabs(solution->data.F64[i]/solutionErr->data.F64[i]));2161 if (fabs(solution->data.F64[i] / solutionErr->data.F64[i]) < COEFF_SIG) {2162 // XXX fprintf (stderr, "drop\n");2163 kernels->solution1->data.F64[i] = 0.0;2164 } else {2165 // XXX fprintf (stderr, "keep\n");2166 kernels->solution1->data.F64[i] = solution->data.F64[i];2167 }2168 }2169 }2170 // double chiSquare = p_pmSubSolve_ChiSquare (kernels, stamps);2171 // fprintf (stderr, "chi square Brute = %f\n", chiSquare);2172 2173 psFree(solution);2174 psFree(sumVector);2175 psFree(sumMatrix);2176 # endif2177 2178 #ifdef TESTING2179 // XXX double-check for NAN in data:2180 for (int iy = 0; iy < stamp->matrix->numRows; iy++) {2181 for (int ix = 0; ix < stamp->matrix->numCols; ix++) {2182 if (!isfinite(stamp->matrix->data.F64[iy][ix])) {2183 fprintf (stderr, "WARNING: NAN in matrix\n");2184 }2185 }2186 }2187 for (int ix = 0; ix < stamp->vector->n; ix++) {2188 if (!isfinite(stamp->vector->data.F64[ix])) {2189 fprintf (stderr, "WARNING: NAN in vector\n");2190 }2191 }2192 #endif2193 2194 #ifdef TESTING2195 for (int ix = 0; ix < sumVector->n; ix++) {2196 if (!isfinite(sumVector->data.F64[ix])) {2197 fprintf (stderr, "WARNING: NAN in vector\n");2198 }2199 }2200 #endif2201 2202 #ifdef TESTING2203 for (int ix = 0; ix < sumVector->n; ix++) {2204 if (!isfinite(sumVector->data.F64[ix])) {2205 fprintf (stderr, "WARNING: NAN in vector\n");2206 }2207 }2208 {2209 psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);2210 psFitsWriteImageSimple("matrixInv.fits", inverse, NULL);2211 psFree(inverse);2212 }2213 {2214 psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);2215 psImage *Xt = psMatrixTranspose(NULL, X);2216 psImage *XtX = psMatrixMultiply(NULL, Xt, X);2217 psFitsWriteImageSimple("matrixErr.fits", XtX, NULL);2218 psFree(X);2219 psFree(Xt);2220 psFree(XtX);2221 }2222 #endif2223
Note:
See TracChangeset
for help on using the changeset viewer.
