Changeset 25834
- Timestamp:
- Oct 13, 2009, 11:30:52 PM (17 years ago)
- Location:
- branches/pap/psModules/src/imcombine
- Files:
-
- 3 edited
-
pmSubtraction.c (modified) (2 diffs)
-
pmSubtractionAnalysis.c (modified) (3 diffs)
-
pmSubtractionEquation.c (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/pap/psModules/src/imcombine/pmSubtraction.c
r25833 r25834 964 964 } 965 965 966 #if 0966 #if 1 967 967 psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) 968 968 { … … 972 972 PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); 973 973 974 psArray *images = psArrayAlloc(kernels->solution1->n - 1); // Images of each kernel to return 975 psVector *fakeSolution = psVectorAlloc(kernels->solution1->n, PS_TYPE_F64); // Fake solution vector 976 psVectorInit(fakeSolution, 0.0); 977 978 for (int i = 0; i < kernels->solution1->n - 1; i++) { 979 fakeSolution->data.F64[i] = kernels->solution1->data.F64[i]; 974 psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution of interest 975 psVector *backup = psVectorCopy(NULL, solution, PS_TYPE_F64); // Backup version 976 977 int num = wantDual ? solution->n - 1 : solution->n; // Number of kernel basis functions 978 979 psArray *images = psArrayAlloc(num); // Images of each kernel to return 980 psVectorInit(solution, 0.0); 981 982 for (int i = 0; i < num; i++) { 983 solution->data.F64[i] = backup->data.F64[i]; 980 984 images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); 981 fakeSolution->data.F64[i] = 0.0;982 } 983 984 psFree( fakeSolution);985 solution->data.F64[i] = 0.0; 986 } 987 psVectorCopy(solution, backup, PS_TYPE_F64); 988 psFree(backup); 985 989 986 990 return images; -
branches/pap/psModules/src/imcombine/pmSubtractionAnalysis.c
r25279 r25834 16 16 #define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image 17 17 18 //#define TESTING 18 19 19 20 bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, … … 117 118 118 119 119 #if 0120 #ifdef TESTING 120 121 // Generate images of the kernel components 121 122 { … … 128 129 } 129 130 psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, false); 130 psFits *kernelFile = psFitsOpen("kernels.fits", "w"); 131 psFits *kernelFile = psFitsOpen("kernels1.fits", "w"); 132 (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); 133 psFitsClose(kernelFile); 134 psFree(kernelImages); 135 psFree(header); 136 } 137 if (kernels->solution2) { 138 psMetadata *header = psMetadataAlloc(); // Header 139 for (int i = 0; i < kernels->solution2->n; i++) { 140 psString name = NULL; // Header keyword 141 psStringAppend(&name, "SOLN%04d", i); 142 psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, kernels->solution2->data.F64[i]); 143 psFree(name); 144 } 145 psArray *kernelImages = pmSubtractionKernelSolutions(kernels, 0.0, 0.0, true); 146 psFits *kernelFile = psFitsOpen("kernels2.fits", "w"); 131 147 (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); 132 148 psFitsClose(kernelFile); -
branches/pap/psModules/src/imcombine/pmSubtractionEquation.c
r25833 r25834 19 19 #define USE_VARIANCE // Include variance in equation? 20 20 21 //#define OLD_FUNCTIONS // Use old functions 22 21 23 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 22 24 // Private (file-static) functions 23 25 ////////////////////////////////////////////////////////////////////////////////////////////////////////////// 24 26 25 #if 127 #ifndef OLD_FUNCTIONS 26 28 // Calculate the least-squares matrix and vector 27 29 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated … … 136 138 double ir = in * ref; 137 139 double rr = PS_SQR(ref); 140 double one = 1.0; 138 141 #ifdef USE_VARIANCE 139 142 float varVal = 1.0 / variance->kernel[y][x]; … … 142 145 in *= varVal; 143 146 ref *= varVal; 147 one *= varVal; 144 148 #endif 145 149 sumRR += rr; … … 147 151 sumR += ref; 148 152 sumI += in; 149 sum1 += 1.0;153 sum1 += one; 150 154 } 151 155 } … … 159 163 } 160 164 161 #if 0162 165 // Calculate the least-squares matrix and vector for dual convolution 163 166 static bool calculateDualMatrixVector(psImage *matrix1, // Least-squares matrix, updated … … 166 169 psVector *vector2, // Least-squares vector, updated 167 170 psImage *matrixX, // Cross-matrix 168 const psKernel *i nput, // Input image (target)169 const psKernel * reference, // Reference image (convolution source)171 const psKernel *image1, // Image 1 172 const psKernel *image2, // Image 2 170 173 const psKernel *variance, // Variance image 171 const psArray *convolutions, // Convolutions for each kernel 174 const psArray *convolutions1, // Convolutions of image 1 for each kernel 175 const psArray *convolutions2, // Convolutions of image 2 for each kernel 172 176 const pmSubtractionKernels *kernels, // Kernels 173 177 const psImage *polyValues, // Spatial polynomial values … … 175 179 ) 176 180 { 177 // (I - R * sum_i a_i k_i - g) (R * k_j) = 0 178 // I C_j = sum_i C_i C_j 179 180 // Background: C_i = 1.0 181 // Normalisation: C_i = R 181 // A_ij = A_i A_j 182 // B_ij = B_i B_j 183 // C_ij = A_i B_j 184 // d_i = A_i I_2 185 // e_i = B_i I_2 186 187 // A_i = I_1 * k_i 188 // B_i = I_2 * k_i 189 190 // Background: A_i = 1.0 191 // Normalisation: A_i = I_1 182 192 183 193 int numKernels = kernels->num; // Number of kernels … … 205 215 206 216 for (int i = 0; i < numKernels; i++) { 207 psKernel *iConv = convolutions->data[i]; // Convolution for index i 217 psKernel *iConv1 = convolutions1->data[i]; // Convolution 1 for index i 218 psKernel *iConv2 = convolutions2->data[i]; // Convolution 2 for index i 208 219 for (int j = i; j < numKernels; j++) { 209 psKernel *jConv = convolutions->data[j]; // Convolution for index j 210 211 double sumCC = 0.0; // Sum of convolution products 220 psKernel *jConv1 = convolutions1->data[j]; // Convolution 1 for index j 221 psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j 222 223 double sumAA = 0.0; // Sum of convolution products for matrix A 224 double sumBB = 0.0; // Sum of convolution products for matrix B 225 double sumAB = 0.0; // Sum of convolution products for matrix C 212 226 for (int y = - footprint; y <= footprint; y++) { 213 227 for (int x = - footprint; x <= footprint; x++) { 214 double cc = iConv->kernel[y][x] * jConv->kernel[y][x]; 228 double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x]; 229 double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; 230 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 215 231 #ifdef USE_VARIANCE 216 cc /= variance->kernel[y][x]; 217 #endif 218 sumCC += cc; 232 aa /= variance->kernel[y][x]; 233 bb /= variance->kernel[y][x]; 234 ab /= variance->kernel[y][x]; 235 #endif 236 sumAA += aa; 237 sumBB += bb; 238 sumAB += ab; 219 239 } 220 240 } … … 223 243 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 224 244 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 225 double value = sumCC * poly2[iTerm][jTerm]; 226 matrix->data.F64[iIndex][jIndex] = value; 227 matrix->data.F64[jIndex][iIndex] = value; 228 } 229 } 230 } 231 232 double sumRC = 0.0; // Sum of the reference-convolution products 233 double sumIC = 0.0; // Sum of the input-convolution products 234 double sumC = 0.0; // Sum of the convolution 245 double aa = sumAA * poly2[iTerm][jTerm]; 246 double bb = sumBB * poly2[iTerm][jTerm]; 247 double ab = sumAB * poly2[iTerm][jTerm]; 248 matrix1->data.F64[iIndex][jIndex] = aa; 249 matrix1->data.F64[jIndex][iIndex] = aa; 250 matrix2->data.F64[iIndex][jIndex] = bb; 251 matrix2->data.F64[jIndex][iIndex] = bb; 252 matrixX->data.F64[iIndex][jIndex] = ab; 253 } 254 } 255 } 256 for (int j = 0; j < i; j++) { 257 psKernel *jConv2 = convolutions2->data[j]; // Convolution 2 for index j 258 double sumAB = 0.0; // Sum of convolution products for matrix C 259 for (int y = - footprint; y <= footprint; y++) { 260 for (int x = - footprint; x <= footprint; x++) { 261 double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x]; 262 #ifdef USE_VARIANCE 263 ab /= variance->kernel[y][x]; 264 #endif 265 sumAB += ab; 266 } 267 } 268 269 // Spatial variation 270 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 271 for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) { 272 double ab = sumAB * poly2[iTerm][jTerm]; 273 matrixX->data.F64[iIndex][jIndex] = ab; 274 } 275 } 276 } 277 278 double sumAI2 = 0.0; // Sum of A.I_2 products (for vector 1) 279 double sumBI2 = 0.0; // Sum of B.I_2 products (for vector 2) 280 double sumAI1 = 0.0; // Sum of A.I_1 products (for matrix 1, normalisation) 281 double sumA = 0.0; // Sum of A (for matrix 1, background) 282 double sumBI1 = 0.0; // Sum of B.I_1 products (for matrix X, normalisation) 283 double sumB = 0.0; // Sum of B products (for matrix X, background) 284 double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) 285 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) 235 286 for (int y = - footprint; y <= footprint; y++) { 236 287 for (int x = - footprint; x <= footprint; x++) { 237 float conv = iConv->kernel[y][x]; 238 float in = input->kernel[y][x]; 239 float ref = reference->kernel[y][x]; 240 double ic = in * conv; 241 double rc = ref * conv; 242 double c = conv; 288 float a = iConv1->kernel[y][x]; 289 float b = iConv2->kernel[y][x]; 290 float i1 = image1->kernel[y][x]; 291 float i2 = image2->kernel[y][x]; 292 293 double ai2 = a * i2; 294 double bi2 = b * i2; 295 double ai1 = a * i1; 296 double bi1 = b * i1; 297 double i1i2 = i1 * i2; 298 243 299 #ifdef USE_VARIANCE 244 300 float varVal = 1.0 / variance->kernel[y][x]; 245 ic *= varVal; 246 rc *= varVal; 247 c *= varVal; 248 #endif 249 sumIC += ic; 250 sumRC += rc; 251 sumC += c; 301 ai2 *= varVal; 302 bi2 *= varVal; 303 ai1 *= varVal; 304 bi1 *= varVal; 305 i1i2 *= varVal; 306 a *= varVal; 307 b *= varVal; 308 i2 *= varVal; 309 #endif 310 311 sumAI2 += ai2; 312 sumBI2 += bi2; 313 sumAI1 += ai1; 314 sumA += a; 315 sumBI1 += bi1; 316 sumB += b; 317 sumI2 += i2; 318 sumI1I2 += i1i2; 252 319 } 253 320 } 254 321 // Spatial variation 255 322 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) { 256 double normTerm = sumRC * poly[iTerm]; 257 double bgTerm = sumC * poly[iTerm]; 258 matrix->data.F64[iIndex][normIndex] = normTerm; 259 matrix->data.F64[normIndex][iIndex] = normTerm; 260 matrix->data.F64[iIndex][bgIndex] = bgTerm; 261 matrix->data.F64[bgIndex][iIndex] = bgTerm; 262 vector->data.F64[iIndex] = sumIC * poly[iTerm]; 263 } 264 } 265 266 double sumRR = 0.0; // Sum of the reference product 267 double sumIR = 0.0; // Sum of the input-reference product 268 double sum1 = 0.0; // Sum of the background 269 double sumR = 0.0; // Sum of the reference 270 double sumI = 0.0; // Sum of the input 323 double ai2 = sumAI2 * poly[iTerm]; 324 double bi2 = sumBI2 * poly[iTerm]; 325 double ai1 = sumAI1 * poly[iTerm]; 326 double a = sumA * poly[iTerm]; 327 double bi1 = sumBI1 * poly[iTerm]; 328 double b = sumB * poly[iTerm]; 329 330 matrix1->data.F64[iIndex][normIndex] = ai1; 331 matrix1->data.F64[normIndex][iIndex] = ai1; 332 matrix1->data.F64[iIndex][bgIndex] = a; 333 matrix1->data.F64[bgIndex][iIndex] = a; 334 vector1->data.F64[iIndex] = ai2; 335 vector2->data.F64[iIndex] = bi2; 336 matrixX->data.F64[iIndex][normIndex] = bi1; 337 matrixX->data.F64[iIndex][bgIndex] = b; 338 } 339 } 340 341 double sumI1 = 0.0; // Sum of I_1 (for matrix 1, background-normalisation) 342 double sumI1I1 = 0.0; // Sum of I_1^2 (for matrix 1, normalisation-normalisation) 343 double sum1 = 0.0; // Sum of 1 (for matrix 1, background-background) 344 double sumI2 = 0.0; // Sum of I_2 (for vector 1, background) 345 double sumI1I2 = 0.0; // Sum of I_1.I_2 (for vector 1, normalisation) 271 346 for (int y = - footprint; y <= footprint; y++) { 272 347 for (int x = - footprint; x <= footprint; x++) { 273 double in = input->kernel[y][x]; 274 double ref = reference->kernel[y][x]; 275 double ir = in * ref; 276 double rr = PS_SQR(ref); 348 float i1 = image1->kernel[y][x]; 349 float i2 = image2->kernel[y][x]; 350 351 double i1i1 = i1 * i1; 352 double one = 1.0; 353 double i1i2 = i1 * i2; 354 277 355 #ifdef USE_VARIANCE 278 356 float varVal = 1.0 / variance->kernel[y][x]; 279 rr *= varVal; 280 ir *= varVal; 281 in *= varVal; 282 ref *= varVal; 283 #endif 284 sumRR += rr; 285 sumIR += ir; 286 sumR += ref; 287 sumI += in; 288 sum1 += 1.0; 289 } 290 } 291 matrix->data.F64[normIndex][normIndex] = sumRR; 292 matrix->data.F64[bgIndex][bgIndex] = sum1; 293 matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR; 294 vector->data.F64[normIndex] = sumIR; 295 vector->data.F64[bgIndex] = sumI; 357 i1 *= varVal; 358 i1i1 *= varVal; 359 one *= varVal; 360 i2 *= varVal; 361 i1i2 *= varVal; 362 #endif 363 364 sumI1 += i1; 365 sumI1I1 += i1i1; 366 sum1 += one; 367 sumI2 += i2; 368 sumI1I2 += i1i2; 369 } 370 } 371 matrix1->data.F64[bgIndex][normIndex] = sumI1; 372 matrix1->data.F64[normIndex][bgIndex] = sumI1; 373 matrix1->data.F64[normIndex][normIndex] = sumI1I1; 374 matrix1->data.F64[bgIndex][bgIndex] = sum1; 375 vector1->data.F64[bgIndex] = sumI2; 376 vector1->data.F64[normIndex] = sumI1I2; 296 377 297 378 return true; 298 379 } 299 #endif300 380 #endif 301 381 … … 321 401 } 322 402 403 #ifdef OLD_FUNCTIONS 323 404 // Calculate a single element of the least-squares matrix, with the polynomial expansions in one direction 324 405 static inline bool calculateMatrixElement1(psImage *matrix, // Matrix to calculate … … 684 765 return true; 685 766 } 686 767 #endif 687 768 688 769 // Add in penalty term to least-squares vector … … 860 941 switch (kernels->mode) { 861 942 case PM_SUBTRACTION_MODE_1: 862 #if 0943 #ifdef OLD_FUNCTIONS 863 944 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 864 945 stamp->variance, polyValues, footprint, true); … … 872 953 break; 873 954 case PM_SUBTRACTION_MODE_2: 874 #if 0955 #ifdef OLD_FUNCTIONS 875 956 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, 876 957 stamp->variance, polyValues, footprint, true); … … 894 975 psVectorInit(stamp->vector2, NAN); 895 976 #endif 977 #ifdef OLD_FUNCTIONS 896 978 status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, 897 979 stamp->variance, polyValues, footprint, true); … … 905 987 status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, 906 988 stamp->image2, stamp->variance, polyValues, footprint, false); 989 #else 990 status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2, 991 stamp->matrixX, stamp->image1, stamp->image2, stamp->variance, 992 stamp->convolutions1, stamp->convolutions2, kernels, polyValues, 993 footprint); 994 #endif 907 995 break; 908 996 default:
Note:
See TracChangeset
for help on using the changeset viewer.
