- Timestamp:
- Jan 8, 2010, 2:24:01 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionVisual.c
r26505 r26547 242 242 243 243 for (int i = 0; i < kernels->num; i++) { 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 } 262 244 pmSubtractionKernelPreCalc *preCalc = kernels->preCalc->data[i]; 245 psKernel *kernel = preCalc->kernel; 246 247 int xSub = i % NXsub; 248 int ySub = i / NXsub; 249 250 int xPix = xSub * (2*footprint + 1 + 3) + footprint; 251 int yPix = ySub * (2*footprint + 1 + 3) + footprint; 252 253 double sum = 0.0; 254 for (int y = -footprint; y <= footprint; y++) { 255 for (int x = -footprint; x <= footprint; x++) { 256 output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; 257 sum += kernel->kernel[y][x]; 258 } 259 } 260 fprintf (stderr, "kernel %d, sum %f\n", i, sum); 261 } 262 263 263 pmVisualScaleImage(kapa1, output, "Image", 0, true); 264 264 pmVisualAskUser(&plotImage); … … 280 280 float maxFlux = NAN; 281 281 for (int i = 0; i < stamps->num; i++) { 282 283 284 285 286 287 288 289 290 291 292 293 282 pmSubtractionStamp *stamp = stamps->stamps->data[i]; 283 if (!isfinite(stamp->flux)) continue; 284 if (!stamp->convolutions1 && !stamp->convolutions2) continue; 285 if (!maxStamp) { 286 maxFlux = stamp->flux; 287 maxStamp = stamp; 288 continue; 289 } 290 if (stamp->flux > maxFlux) { 291 maxFlux = stamp->flux; 292 maxStamp = stamp; 293 } 294 294 } 295 295 296 296 if (!isfinite(maxStamp->flux)) { 297 297 fprintf (stderr, "no valid stamps?\n"); 298 298 } 299 299 … … 301 301 302 302 if (maxStamp->convolutions1) { 303 304 305 306 307 308 309 310 311 312 313 314 315 303 // output image is a grid of NXsub by NYsub sub-images 304 nKernels = maxStamp->convolutions1->n; 305 int NXsub = sqrt(nKernels); 306 int NYsub = nKernels / NXsub; 307 if (nKernels % NXsub) NYsub++; 308 309 int NXpix = NXsub * (2*footprint + 1 + 3); 310 int NYpix = NYsub * (2*footprint + 1 + 3); 311 312 psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); 313 psImageInit (output, 0.0); 314 315 for (int i = 0; i < nKernels; i++) { 316 316 psKernel *kernel = maxStamp->convolutions1->data[i]; 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 } 333 334 } 335 317 318 int xSub = i % NXsub; 319 int ySub = i / NXsub; 320 321 int xPix = xSub * (2*footprint + 1 + 3) + footprint; 322 int yPix = ySub * (2*footprint + 1 + 3) + footprint; 323 324 double sum = 0.0; 325 for (int y = -footprint; y <= footprint; y++) { 326 for (int x = -footprint; x <= footprint; x++) { 327 output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; 328 sum += kernel->kernel[y][x]; 329 } 330 } 331 fprintf (stderr, "kernel %d, sum %f\n", i, sum); 332 } 333 pmVisualScaleImage(kapa2, output, "Image", 0, true); 334 } 335 336 336 if (maxStamp->convolutions2) { 337 338 339 340 341 342 343 344 345 346 347 348 349 337 // output image is a grid of NXsub by NYsub sub-images 338 nKernels = maxStamp->convolutions2->n; 339 int NXsub = sqrt(nKernels); 340 int NYsub = nKernels / NXsub; 341 if (nKernels % NXsub) NYsub++; 342 343 int NXpix = NXsub * (2*footprint + 1 + 3); 344 int NYpix = NYsub * (2*footprint + 1 + 3); 345 346 psImage *output = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); 347 psImageInit (output, 0.0); 348 349 for (int i = 0; i < nKernels; i++) { 350 350 psKernel *kernel = maxStamp->convolutions2->data[i]; 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 } 367 368 } 369 351 352 int xSub = i % NXsub; 353 int ySub = i / NXsub; 354 355 int xPix = xSub * (2*footprint + 1 + 3) + footprint; 356 int yPix = ySub * (2*footprint + 1 + 3) + footprint; 357 358 double sum = 0.0; 359 for (int y = -footprint; y <= footprint; y++) { 360 for (int x = -footprint; x <= footprint; x++) { 361 output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; 362 sum += kernel->kernel[y][x]; 363 } 364 } 365 fprintf (stderr, "kernel %d, sum %f\n", i, sum); 366 } 367 pmVisualScaleImage(kapa2, output, "Image", 1, true); 368 } 369 370 370 pmVisualAskUser(&plotImage); 371 371 return true; … … 394 394 395 395 overlay[Noverlay].type = KII_OVERLAY_BOX; 396 397 398 399 400 401 402 403 396 if ((stamp->x < 1.0) && (stamp->y < 1.0)) { 397 // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y); 398 continue; 399 } 400 if (!isfinite(stamp->x) && !isfinite(stamp->y)) { 401 // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y); 402 continue; 403 } 404 404 overlay[Noverlay].x = stamp->x; 405 405 overlay[Noverlay].y = stamp->y; … … 437 437 float NXf = sqrt(stamps->num); 438 438 NX = (int) NXf == NXf ? NXf : NXf + 1.0; 439 439 440 440 float NYf = stamps->num / NX; 441 441 NY = (int) NYf == NY ? NYf : NYf + 1.0; … … 453 453 differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); 454 454 convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); 455 455 456 456 psImageInit (sourceImage, 0.0); 457 457 psImageInit (targetImage, 0.0); … … 479 479 sum = 0.0; 480 480 for (int y = -footprint; y <= footprint; y++) { 481 482 483 484 481 for (int x = -footprint; x <= footprint; x++) { 482 targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x]; 483 sum += targetImage->data.F32[y + NYpix][x + NXpix]; 484 } 485 485 } 486 486 targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; … … 489 489 sum = 0.0; 490 490 for (int y = -footprint; y <= footprint; y++) { 491 492 493 494 491 for (int x = -footprint; x <= footprint; x++) { 492 sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x]; 493 sum += sourceImage->data.F32[y + NYpix][x + NXpix]; 494 } 495 495 } 496 496 sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; … … 499 499 sum = 0.0; 500 500 for (int y = -footprint; y <= footprint; y++) { 501 502 503 504 501 for (int x = -footprint; x <= footprint; x++) { 502 convolutionImage->data.F32[y + NYpix][x + NXpix] = -convolution->kernel[y][x]; 503 sum += convolutionImage->data.F32[y + NYpix][x + NXpix]; 504 } 505 505 } 506 506 convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; 507 507 508 508 // insert the (difference) kernel into the (difference) image: 509 509 sum = 0.0; 510 510 for (int y = -footprint; y <= footprint; y++) { 511 512 513 514 511 for (int x = -footprint; x <= footprint; x++) { 512 differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm; 513 sum += differenceImage->data.F32[y + NYpix][x + NXpix]; 514 } 515 515 } 516 516 differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; … … 519 519 sum = 0.0; 520 520 for (int y = -footprint; y <= footprint; y++) { 521 522 523 524 521 for (int x = -footprint; x <= footprint; x++) { 522 residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm + convolution->kernel[y][x]; 523 sum += residualImage->data.F32[y + NYpix][x + NXpix]; 524 } 525 525 } 526 526 residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; … … 528 528 // insert the (fresidual) kernel into the (fresidual) image: 529 529 for (int y = -footprint; y <= footprint; y++) { 530 531 532 530 for (int x = -footprint; x <= footprint; x++) { 531 fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); 532 } 533 533 } 534 534 return true; … … 536 536 537 537 bool pmSubtractionVisualShowFit(double norm) { 538 539 if (!pmVisualIsVisual()) return true; 538 540 539 541 // XXX a dumb test : dump the residual image and exit 540 542 { 541 543 psMetadata *header = psMetadataAlloc(); 542 544 psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm); 543 psFits *fits = psFitsOpen("resid.fits", "w"); 544 psFitsWriteImage(fits, header, residualImage, 0, NULL); 545 psFitsClose(fits); 546 exit (0); 547 } 548 549 // if (!pmVisualIsVisual()) return true; 545 psFits *fits = psFitsOpen("resid.fits", "w"); 546 psFitsWriteImage(fits, header, residualImage, 0, NULL); 547 psFitsClose(fits); 548 exit (0); 549 } 550 550 551 if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; 551 552 if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false; … … 606 607 for (int i = 0; i < kernels->num; i++) { 607 608 x->data.F32[i] = i; 608 609 y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false); 609 610 graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]); 610 611 graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
Note:
See TracChangeset
for help on using the changeset viewer.
