- Timestamp:
- Jan 19, 2006, 4:38:28 PM (21 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.c
r5867 r6062 11 11 * @author GLG, MHPCC 12 12 * 13 * @version $Revision: 1.6.8.1.2. 1$ $Name: not supported by cvs2svn $14 * @date $Date: 200 5-12-31 04:35:58 $13 * @version $Revision: 1.6.8.1.2.2 $ $Name: not supported by cvs2svn $ 14 * @date $Date: 2006-01-20 02:38:28 $ 15 15 * 16 16 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 22 22 #endif 23 23 24 #include <assert.h> 24 25 #include "pmSubtractBias.h" 25 26 26 27 #define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 27 28 #define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 29 30 31 #define MAX(a,b) ((a) > (b) ? (a) : (b)) 32 #define MIN(a,b) ((a) < (b) ? (a) : (b)) 33 28 34 29 35 // XXX: put these in psConstants.h … … 56 62 }\ 57 63 64 65 void overscanOptionsFree(pmOverscanOptions *options) 66 { 67 psFree(options->stat); 68 psFree(options->poly); 69 psFree(options->spline); 70 } 71 72 pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat) 73 { 74 pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); 75 psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); 76 77 // Inputs 78 opts->single = single; 79 opts->fitType = fitType; 80 opts->order = order; 81 opts->stat = psMemIncrRefCounter(stat); 82 83 // Outputs 84 opts->poly = NULL; 85 opts->spline = NULL; 86 87 return opts; 88 } 89 90 58 91 /****************************************************************************** 59 92 psSubtractFrame(): this routine will take as input a readout for the input … … 61 94 place from the input image. 62 95 *****************************************************************************/ 63 static pmReadout *SubtractFrame(pmReadout *in, 64 const pmReadout *bias) 65 { 66 psS32 i; 67 psS32 j; 68 69 if (bias == NULL) { 70 psLogMsg(__func__, PS_LOG_NOTE, 71 "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL. Returning original image.\n"); 72 return(in); 73 } 74 75 76 if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { 77 psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); 78 return(in); 79 } 80 if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { 81 psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); 82 return(in); 83 } 84 85 for (i=0;i<in->image->numRows;i++) { 86 for (j=0;j<in->image->numCols;j++) { 87 in->image->data.F32[i][j]-= 88 bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0]; 89 if ((in->mask != NULL) && (bias->mask != NULL)) { 90 (in->mask->data.U8[i][j])|= 91 bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0]; 96 static bool SubtractFrame(pmReadout *in,// Input readout 97 const pmReadout *sub, // Readout to be subtracted from input 98 float scale // Scale to apply before subtracting 99 ) 100 { 101 assert(in); 102 assert(sub); 103 104 // Get the trim sections 105 psRegion *inTrimsec = psMetadataLookupPtr(NULL, in->parent->concepts, "CELL.TRIMSEC"); 106 psRegion *subTrimsec = psMetadataLookupPtr(NULL, sub->parent->concepts, "CELL.TRIMSEC"); 107 psImage *inImage = psImageSubset(in->image, *inTrimsec); // The input image 108 psImage *subImage = psImageSubset(sub->image, *subTrimsec); // The image to be subtracted 109 psImage *inMask = in->mask ? psImageSubset(in->mask, *inTrimsec) : NULL; // The input mask 110 psImage *subMask = sub->mask ? psImageSubset(sub->mask, *subTrimsec) : NULL; // The input mask 111 112 // Offsets of the cells 113 int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); 114 int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); 115 int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0"); 116 int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0"); 117 118 if ((inImage->numCols + x0in - x0sub) > subImage->numCols) { 119 psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n"); 120 return false; 121 } 122 if ((inImage->numRows + y0in - y0sub) > subImage->numRows) { 123 psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n"); 124 return false; 125 } 126 127 if (scale == 1.0) { 128 for (int i = 0; i < inImage->numRows; i++) { 129 for (int j = 0; j < inImage->numCols; j++) { 130 inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub]; 131 if (inMask && subMask) { 132 inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; 133 } 92 134 } 93 135 } 94 } 95 96 return(in); 97 } 98 136 } else { 137 for (int i = 0; i < inImage->numRows; i++) { 138 for (int j = 0; j < inImage->numCols; j++) { 139 inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale; 140 if (inMask && subMask) { 141 inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; 142 } 143 } 144 } 145 } 146 147 return true; 148 } 149 150 151 #if 0 99 152 /****************************************************************************** 100 153 ImageSubtractScalar(): subtract a scalar from the input image. … … 114 167 return(image); 115 168 } 169 #endif 116 170 117 171 /****************************************************************************** … … 180 234 181 235 236 #if 0 182 237 /****************************************************************************** 183 238 ScaleOverscanVector(): this routine takes as input an arbitrary vector, … … 283 338 } 284 339 340 #endif 341 342 // Produce an overscan vector from an array of pixels 343 static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options 344 const psArray *pixels, // Array of vectors containing the pixel values 345 psStats *myStats // Statistic to use in reducing the overscan 346 ) 347 { 348 // Reduce the overscans 349 psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row 350 psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate 351 psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting 352 for (int i = 0; i < pixels->n; i++) { 353 psVector *values = pixels->data[i]; // Vector with overscan values 354 if (values->n > 0) { 355 mask->data.U8[i] = 0; 356 ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] 357 psVectorStats(myStats, values, NULL, NULL, 0); 358 double reducedVal = NAN; // Result of statistics 359 if (! p_psGetStatValue(myStats, &reducedVal)) { 360 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result " 361 "of statistics on row %d.\n", i); 362 return NULL; 363 } 364 reduced->data.F32[i] = reducedVal; 365 // printf("%f ", reducedVal); 366 } else if (overscanOpts->fitType == PM_FIT_NONE) { 367 psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " 368 "image, and no fit is requested.\n"); 369 return NULL; 370 } else { 371 // We'll fit this one out 372 mask->data.U8[i] = 1; 373 } 374 375 // printf("\n"); 376 } 377 378 // Fit the overscan, if required 379 switch (overscanOpts->fitType) { 380 case PM_FIT_NONE: 381 // No fitting --- that's easy. 382 break; 383 case PM_FIT_POLY_ORD: 384 overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_ORD); 385 overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, 386 ordinate); 387 psFree(reduced); 388 reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); 389 break; 390 case PM_FIT_POLY_CHEBY: 391 overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_CHEB); 392 overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, 393 ordinate); 394 psFree(reduced); 395 reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); 396 break; 397 case PM_FIT_SPLINE: 398 // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an 399 // input spline 400 overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); 401 psFree(reduced); 402 reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate); 403 break; 404 default: 405 psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); 406 return NULL; 407 break; 408 } 409 410 psFree(ordinate); 411 psFree(mask); 412 413 return reduced; 414 } 415 416 417 285 418 /****************************************************************************** 286 419 XXX: The SDRS does not specify type support. F32 is implemented here. 287 420 *****************************************************************************/ 288 pmReadout *pmSubtractBias(pmReadout *in, 289 void *fitSpec, 290 const psList *overscans, 291 pmOverscanAxis overScanAxis, 292 psStats *stat, 293 psS32 nBinOrig, 294 pmFit fit, 295 const pmReadout *bias) 421 pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, 422 const pmReadout *bias, const pmReadout *dark) 296 423 { 297 424 psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, … … 301 428 PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); 302 429 303 // 304 // If the overscans != NULL, then check the type of each image. 305 // 306 if (overscans != NULL) { 307 psListElem *tmpOverscan = (psListElem *) overscans->head; 308 while (NULL != tmpOverscan) { 309 psImage *myOverscanImage = (psImage *) tmpOverscan->data; 310 PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); 311 tmpOverscan = tmpOverscan->next; 312 } 313 } 314 315 if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) { 316 psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE). Returning in image\n"); 317 return(in); 318 } 319 320 // Check for an unallowable pmFit. 321 if ((fit != PM_OVERSCAN_NONE) && 322 (fit != PM_OVERSCAN_ROWS) && 323 (fit != PM_OVERSCAN_COLUMNS) && 324 (fit != PM_OVERSCAN_ALL)) { 325 psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d). Returning in image.\n", fit); 326 return(in); 327 } 328 // Check for an unallowable pmOverscanAxis. 329 if ((overScanAxis != PM_OVERSCAN_NONE) && 330 (overScanAxis != PM_OVERSCAN_ROWS) && 331 (overScanAxis != PM_OVERSCAN_COLUMNS) && 332 (overScanAxis != PM_OVERSCAN_ALL)) { 333 psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d). Returning in image.\n", overScanAxis); 334 return(in); 335 } 336 psS32 i; 337 psS32 j; 338 psS32 numBins = 0; 339 static psVector *overscanVector = NULL; 340 psVector *tmpRow = NULL; 341 psVector *tmpCol = NULL; 342 psVector *myBin = NULL; 343 psVector *binVec = NULL; 344 psListElem *tmpOverscan = NULL; 345 double statValue; 346 psImage *myOverscanImage = NULL; 347 psPolynomial1D *myPoly = NULL; 348 psSpline1D *mySpline = NULL; 349 psS32 nBin; 350 351 // 352 // Create a static stats data structure and determine the highest 353 // priority stats option. 354 // 355 static psStats *myStats = NULL; 356 if (myStats == NULL) { 357 myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); 358 p_psMemSetPersistent(myStats, true); 359 } 360 if (stat != NULL) { 361 myStats->options = GenNewStatOptions(stat); 362 } 363 364 365 if (overScanAxis == PM_OVERSCAN_NONE) { 366 if (fit != PM_FIT_NONE) { 367 psLogMsg(__func__, PS_LOG_WARN, 368 "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE. Proceeding to full fram subtraction.\n"); 369 } 370 371 if (overscans != NULL) { 372 psLogMsg(__func__, PS_LOG_WARN, 373 "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL. Proceeding to full fram subtraction.\n"); 374 } 375 return(SubtractFrame(in, bias)); 376 } 377 378 if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) { 379 psLogMsg(__func__, PS_LOG_WARN, 380 "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE. Proceeding with the rest of the module.\n"); 381 } 382 383 384 // 385 // We subtract each overscan region from the image data. 386 // If we get here we know that overscans != NULL. 387 // 388 389 if (overScanAxis == PM_OVERSCAN_ALL) { 390 tmpOverscan = (psListElem *) overscans->head; 391 while (NULL != tmpOverscan) { 392 myOverscanImage = (psImage *) tmpOverscan->data; 393 394 PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); 395 psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); 396 if (rc == NULL) { 397 psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); 430 pmCell *cell = in->parent; // The parent cell 431 psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // The trim region 432 psImage *image = psImageSubset(in->image, *trimsec); // The image corresponding to the trim region 433 434 // Overscan processing 435 if (overscanOpts) { 436 // Check for an unallowable pmFit. 437 if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && 438 overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { 439 psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", overscanOpts->fitType); 440 return(in); 441 } 442 443 // Get the list of overscans 444 psList *overscanRegions = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); 445 psList *overscans = psListAlloc(NULL); // List of the overscan images 446 psListIterator *iter = psListIteratorAlloc(overscanRegions, PS_LIST_HEAD, false); // Iterator 447 psRegion *biassec = NULL; // A BIASSEC region from the list 448 while ((biassec = psListGetAndIncrement(iter))) { 449 psImage *overscan = psImageSubset(in->image, *biassec); 450 psListAdd(overscans, PS_LIST_TAIL, overscan); 451 psFree(overscan); 452 } 453 psFree(iter); 454 455 psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original 456 myStats->options = GenNewStatOptions(overscanOpts->stat); 457 458 // Reduce all overscan pixels to a single value 459 if (overscanOpts->single) { 460 psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); 461 pixels->n = 0; 462 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 463 psImage *overscan = NULL; // Overscan image from iterator 464 while ((overscan = psListGetAndIncrement(iter))) { 465 int index = pixels->n; // Index 466 pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); 467 // XXX Reimplement with memcpy 468 for (int i = 0; i < overscan->numRows; i++) { 469 for (int j = 0; j < overscan->numCols; j++) { 470 pixels->data.F32[index++] = overscan->data.F32[i][j]; 471 } 472 } 473 474 } 475 psFree(iter); 476 477 (void)psVectorStats(myStats, pixels, NULL, NULL, 0); 478 double reduced = NAN; // Result of statistics 479 if (! p_psGetStatValue(myStats, &reduced)) { 480 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning input image.\n"); 398 481 return(in); 399 482 } 400 if (false == p_psGetStatValue(myStats, &statValue)) { 401 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 402 return(in); 483 (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); 484 } else { 485 486 // We do the regular overscan subtraction 487 bool readRows = psMetadataLookupBool(NULL, cell->concepts, "CELL.READDIR"); // Read direction 488 489 if (readRows) { 490 // The read direction is rows 491 psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels 492 for (int i = 0; i < pixels->n; i++) { 493 psVector *values = psVectorAlloc(0, PS_TYPE_F32); 494 values->n = 0; 495 pixels->data[i] = values; 496 } 497 498 // Pull the pixels out into the vectors 499 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 500 psImage *overscan = NULL; // Overscan image from iterator 501 while ((overscan = psListGetAndIncrement(iter))) { 502 int diff = image->row0 - overscan->row0; // Offset between the two regions 503 for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) { 504 // i is row on overscan 505 // XXX Reimplement with memcpy 506 psVector *values = pixels->data[i]; 507 int index = values->n; // Index in the vector 508 values = psVectorRealloc(values, values->n + overscan->numCols); 509 for (int j = 0; j < overscan->numCols; j++) { 510 values->data.F32[index++] = overscan->data.F32[i][j]; 511 } 512 values->n += overscan->numCols; 513 pixels->data[i] = values; // Update the pointer in case it's moved 514 } 515 } 516 517 // Reduce the overscans 518 psVector *reduced = overscanVector(overscanOpts, pixels, myStats); 519 if (! reduced) { 520 return in; 521 } 522 523 // Subtract row by row 524 for (int i = 0; i < image->numRows; i++) { 525 for (int j = 0; j < image->numCols; j++) { 526 image->data.F32[i][j] -= reduced->data.F32[i]; 527 } 528 } 529 psFree(reduced); 530 531 } else { 532 // The read direction is columns 533 psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels 534 for (int i = 0; i < pixels->n; i++) { 535 psVector *values = psVectorAlloc(0, PS_TYPE_F32); 536 values->n = 0; 537 pixels->data[i] = values; 538 } 539 540 // Pull the pixels out into the vectors 541 psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator 542 psImage *overscan = NULL; // Overscan image from iterator 543 while ((overscan = psListGetAndIncrement(iter))) { 544 int diff = image->col0 - overscan->col0; // Offset between the two regions 545 for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) { 546 // i is column on overscan 547 // XXX Reimplement with memcpy 548 psVector *values = pixels->data[i]; 549 int index = values->n; // Index in the vector 550 values = psVectorRealloc(values, values->n + overscan->numRows); 551 for (int j = 0; j < overscan->numRows; j++) { 552 values->data.F32[index++] = overscan->data.F32[i][j]; 553 } 554 values->n += overscan->numRows; 555 pixels->data[i] = values; // Update the pointer in case it's moved 556 } 557 } 558 559 // Reduce the overscans 560 psVector *reduced = overscanVector(overscanOpts, pixels, myStats); 561 if (! reduced) { 562 return in; 563 } 564 565 // Subtract column by column 566 for (int i = 0; i < image->numCols; i++) { 567 for (int j = 0; j < image->numRows; j++) { 568 image->data.F32[j][i] -= reduced->data.F32[i]; 569 } 570 } 571 psFree(reduced); 403 572 } 404 ImageSubtractScalar(in->image, statValue); 405 406 tmpOverscan = tmpOverscan->next; 407 } 408 return(in); 409 } 410 411 // This check is redundant with above code. 412 if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) { 413 psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis); 414 return(in); 415 } 416 417 tmpOverscan = (psListElem *) overscans->head; 418 while (NULL != tmpOverscan) { 419 // PS_IMAGE_PRINT_F32_HIDEF(in->image); 420 myOverscanImage = (psImage *) tmpOverscan->data; 421 422 if (overScanAxis == PM_OVERSCAN_ROWS) { 423 if (myOverscanImage->numCols != (in->image)->numCols) { 424 psLogMsg(__func__, PS_LOG_WARN, 425 "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n", 426 myOverscanImage->numCols, in->image->numCols); 427 } 428 429 // We create a row vector and subtract this vector from image. 430 // XXX: Is there a better way to extract a psVector from a psImage without 431 // having to copy every element in that vector? 432 overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); 433 for (i=0;i<overscanVector->n;i++) { 434 overscanVector->data.F32[i] = 0.0; 435 } 436 tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); 437 438 // For each column of the input image, loop through every row, 439 // collect the pixel in that row, then performed the specified 440 // statistical op on those pixels. Store this in overscanVector. 441 for (i=0;i<myOverscanImage->numCols;i++) { 442 for (j=0;j<myOverscanImage->numRows;j++) { 443 tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i]; 444 } 445 psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); 446 if (rc == NULL) { 447 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); 448 return(in); 449 } 450 if (false == p_psGetStatValue(rc, &statValue)) { 451 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 452 return(in); 453 } 454 overscanVector->data.F32[i] = statValue; 455 } 456 psFree(tmpRow); 457 458 // Scale the overscan vector to the size of the input image. 459 if (overscanVector->n != in->image->numCols) { 460 if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { 461 psVector *newVec = ScaleOverscanVector(overscanVector, 462 in->image->numCols, 463 fitSpec, fit); 464 if (newVec == NULL) { 465 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); 466 return(in); 467 } 468 psFree(overscanVector); 469 overscanVector = newVec; 470 } else { 471 psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); 472 psFree(overscanVector); 473 return(in); 474 } 475 } 476 } 477 478 if (overScanAxis == PM_OVERSCAN_COLUMNS) { 479 if (myOverscanImage->numRows != (in->image)->numRows) { 480 psLogMsg(__func__, PS_LOG_WARN, 481 "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n", 482 myOverscanImage->numRows, in->image->numRows); 483 } 484 485 // We create a column vector and subtract this vector from image. 486 overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); 487 for (i=0;i<overscanVector->n;i++) { 488 overscanVector->data.F32[i] = 0.0; 489 } 490 tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); 491 492 // For each row of the input image, loop through every column, 493 // collect the pixel in that row, then performed the specified 494 // statistical op on those pixels. Store this in overscanVector. 495 for (i=0;i<myOverscanImage->numRows;i++) { 496 for (j=0;j<myOverscanImage->numCols;j++) { 497 tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j]; 498 } 499 psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); 500 if (rc == NULL) { 501 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); 502 return(in); 503 } 504 if (false == p_psGetStatValue(rc, &statValue)) { 505 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 506 return(in); 507 } 508 overscanVector->data.F32[i] = statValue; 509 } 510 psFree(tmpCol); 511 512 // Scale the overscan vector to the size of the input image. 513 if (overscanVector->n != in->image->numRows) { 514 if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { 515 psVector *newVec = ScaleOverscanVector(overscanVector, 516 in->image->numRows, 517 fitSpec, fit); 518 if (newVec == NULL) { 519 psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector. Returning in image.\n"); 520 return(in); 521 } 522 psFree(overscanVector); 523 overscanVector = newVec; 524 } else { 525 psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL. Returning in image.\n"); 526 psFree(overscanVector); 527 return(in); 528 } 529 } 530 } 531 532 // 533 // Re-bin the overscan vector (change its length). 534 // 535 // Only if nBinOrig > 1. 536 if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { 537 numBins = 1+((overscanVector->n)/nBinOrig); 538 myBin = psVectorAlloc(numBins, PS_TYPE_F32); 539 binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); 540 541 for (i=0;i<numBins;i++) { 542 for(j=0;j<nBinOrig;j++) { 543 if (overscanVector->n > ((i*nBinOrig)+j)) { 544 binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; 545 } else { 546 // XXX: we get here if nBinOrig does not evenly divide 547 // the overscanVector vector. This is the last bin. Should 548 // we change the binVec->n to acknowledge that? 549 binVec->n = j; 550 } 551 } 552 psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); 553 if (rc == NULL) { 554 psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); 555 return(in); 556 } 557 if (false == p_psGetStatValue(rc, &statValue)) { 558 psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); 559 return(in); 560 } 561 myBin->data.F32[i] = statValue; 562 } 563 564 // Change the effective size of overscanVector. 565 overscanVector->n = numBins; 566 for (i=0;i<numBins;i++) { 567 overscanVector->data.F32[i] = myBin->data.F32[i]; 568 } 569 psFree(binVec); 570 psFree(myBin); 571 nBin = nBinOrig; 572 } else { 573 nBin = 1; 574 } 575 576 // At this point the number of data points in overscanVector should be 577 // equal to the number of rows/columns (whatever is appropriate) in the 578 // image divided by numBins. 579 // 580 581 582 // 583 // This doesn't seem right. The only way to do a spline fit is if, 584 // by SDRS requirements, fitSpec is not-NULL> But in order for it 585 // to be non-NULL, someone must have called psSpline1DAlloc() with 586 // the min, max, and number of splines. 587 // 588 if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) { 589 // 590 // Fit a polynomial or spline to the overscan vector. 591 // 592 if (fit == PM_FIT_POLYNOMIAL) { 593 myPoly = (psPolynomial1D *) fitSpec; 594 myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); 595 if (myPoly == NULL) { 596 psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector. Returning in image.\n"); 597 psFree(overscanVector); 598 return(in); 599 } 600 } else if (fit == PM_FIT_SPLINE) { 601 // XXX: This makes no sense 602 // XXX: must free mySpline? 603 mySpline = (psSpline1D *) fitSpec; 604 mySpline = psVectorFitSpline1D(NULL, overscanVector); 605 if (mySpline == NULL) { 606 psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning in image.\n"); 607 psFree(overscanVector); 608 return(in); 609 } 610 } 611 612 // 613 // Subtract fitted overscan vector row-wise from the image. 614 // 615 if (overScanAxis == PM_OVERSCAN_ROWS) { 616 for (i=0;i<(in->image)->numCols;i++) { 617 psF32 tmpF32 = 0.0; 618 if (fit == PM_FIT_POLYNOMIAL) { 619 tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); 620 } else if (fit == PM_FIT_SPLINE) { 621 tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); 622 } 623 for (j=0;j<(in->image)->numRows;j++) { 624 (in->image)->data.F32[j][i]-= tmpF32; 625 } 626 } 627 } 628 629 // 630 // Subtract fitted overscan vector column-wise from the image. 631 // 632 if (overScanAxis == PM_OVERSCAN_COLUMNS) { 633 for (i=0;i<(in->image)->numRows;i++) { 634 psF32 tmpF32 = 0.0; 635 if (fit == PM_FIT_POLYNOMIAL) { 636 tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); 637 } else if (fit == PM_FIT_SPLINE) { 638 tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); 639 } 640 641 for (j=0;j<(in->image)->numCols;j++) { 642 (in->image)->data.F32[i][j]-= tmpF32; 643 } 644 } 645 } 646 } else { 647 // 648 // If we get here, then no polynomials were fit to the overscan 649 // vector. We simply subtract it, taking into account binning, 650 // from the image. 651 // 652 653 // 654 // Subtract overscan vector row-wise from the image. 655 // 656 if (overScanAxis == PM_OVERSCAN_ROWS) { 657 for (i=0;i<(in->image)->numCols;i++) { 658 for (j=0;j<(in->image)->numRows;j++) { 659 (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin]; 660 } 661 } 662 } 663 664 // 665 // Subtract overscan vector column-wise from the image. 666 // 667 if (overScanAxis == PM_OVERSCAN_COLUMNS) { 668 for (i=0;i<(in->image)->numRows;i++) { 669 for (j=0;j<(in->image)->numCols;j++) { 670 (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin]; 671 } 672 } 673 } 674 } 675 676 psFree(overscanVector); 677 678 tmpOverscan = tmpOverscan->next; 679 } 680 681 psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, 682 "---- pmSubtractBias() exit ----\n"); 683 684 if (bias != NULL) { 685 return(SubtractFrame(in, bias)); 686 } 687 return(in); 688 } 689 690 573 } 574 psFree(myStats); 575 psFree(overscans); 576 } // End of overscan subtraction 577 578 // Bias frame subtraction 579 if (bias) { 580 SubtractFrame(in, bias, 1.0); 581 } 582 583 if (dark) { 584 // Get the scaling 585 float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); 586 float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); 587 SubtractFrame(in, dark, inTime/darkTime); 588 } 589 590 return in; 591 } 592 593
Note:
See TracChangeset
for help on using the changeset viewer.
