/* Copyright restrictions apply - see stsdas$copyright.stsdas */ # include # include # include "estreak.h" # define ITER 2 # define NSIGMA 2. /* G_MEDIAN: Scale the individual flats and compute their median. * * Adapted from the WFPC version written by J.C. Hsu. * * Image data resides in either temporary or raw files. Output array * 'median' with the computed median image must be allocated by caller. * * * * Revision history: * --------------- * 25 Apr 96 - Implementation (IB) * 07 Oct 96 - Revised after code review (IB) * 28 Oct 96 - Fixed g_select call (IB) * */ int g_median (IOControl *ioc, Bool tmp_file, int ngood, floatArray *median) { /* Parameters: * * IOControl ioc io: I/O control structure * Bool tmp_file i: read from _tmp images instead of _raw ? * int ngood i: min. good points needed * floatArray median o: 2-d median array (allocated by caller) */ /* Local variables */ int j, i, n; unsigned long k, l; int iter, image, remain; int x_w1, x_w2, y_w1, y_w2; /* Window */ float high, low; /* Sigma-clip limits */ double *sum; /* Accumulators */ int *npix; int blocksize; floatArray pix; /* Current image */ floatArray ratio; /* Ratio image */ floatArray block[MAX_FILES]; /* blkSize-line block buffer */ float temp[MAX_FILES+1]; /* Sorting vector */ float avg, stddev, med; float meanRatio[MAX_FILES]; /* Mean ratio for each image */ char *malloc_err = "Median computation: cannot allocate memory."; /* Function declarations */ void g_stat (floatArray *, int, int, int, int, Bool, float, float, float *, float *, int *, Bool); float g_select (unsigned long, unsigned long, float *); /* Window where to compute average. */ x_w1 = (int)(ioc->x_size * WIND1) - 1; x_w2 = (int)(ioc->x_size * WIND2) - 1; y_w1 = (int)(ioc->y_size * WIND1) - 1; y_w2 = (int)(ioc->y_size * WIND2) - 1; /* Step 1: Compute average of input arrays. Result is stored in median array.*/ /* Calloc accumulators. */ sum = (double *) calloc ((size_t)((long)ioc->y_size*(long)ioc->x_size), (size_t)sizeof(double)); if (sum == NULL) { g_error (malloc_err); return (1); } npix = (int *) calloc ((size_t)((long)ioc->y_size*(long)ioc->x_size), (size_t)sizeof(int)); if (npix == NULL) { g_error (malloc_err); return (1); } /* Loop over image list. */ for (image = 0; image < ioc->nimage; image++) { /* Read full image. */ if (g_openImage (ioc, image, tmp_file)) return (1); g_getBlock (ioc, 1, ioc->y_size, &pix); /* Accumulate image. */ for (l = 0; l < pix.bufsize; l++) { if (pix.data[l] != BADVAL) { sum[l] += (double)pix.data[l]; npix[l] += 1; } } /* Close image. */ g_closeImage (ioc); } /* Compute average image. */ for (l = 0; l < pix.bufsize; l++) { if (npix[l] > 0) median->data[l] = (float)(sum[l] / (double)npix[l]); else median->data[l] = BADVAL; } /* Free accumulator memory */ free (npix); free (sum); /* Step 2: Compute mean ratio for each image, with sigma-clipping. */ /* Skip if only one image. */ if (ioc->nimage > 1) { /* Alloc 2-d temporary area for storing pixel ratios. */ if (allocArray (&ratio, ioc->x_size, ioc->y_size)) { g_error (malloc_err); return (1); } /* Loop over image list. */ for (image = 0; image < ioc->nimage; image++) { /* Read full image. */ if (g_openImage (ioc, image, tmp_file)) return (1); g_getBlock (ioc, 1, ioc->y_size, &pix); /* Loop over lines and columns of central region */ for (j = y_w1; j < y_w2; j++) { for (i = x_w1; i < x_w2; i++) { /* Compute ratio and store in temporary area. */ if ((GPPix(median,i,j) != BADVAL) && (GPix(pix,i,j) != BADVAL)) GPix(ratio,i,j) = GPix(pix,i,j) / GPPix(median,i,j); else GPix(ratio,i,j) = BADVAL; } } /* Close image. */ g_closeImage (ioc); /* Compute average and stddev in central window. */ g_stat (&ratio, x_w1, x_w2, y_w1, y_w2, False, 0.0F,0.0F, &avg, &stddev, &n, False); /* Iterate for sigma-clipping. */ for (iter = 1; iter <= ITER; iter++) { /* Compute high and low limits. */ high = avg + NSIGMA * stddev; low = avg - NSIGMA * stddev; /* Compute stats in central region with sigma-clipping. */ g_stat (&ratio, x_w1, x_w2, y_w1, y_w2, True, high, low, &avg, &stddev, &n, False); } /* Store average ratio for current image. */ meanRatio[image] = avg; } /* Free temporary area. */ freeArray (&ratio); } else meanRatio[0] = 1.0F; /* Step 3: Determine the median value for each pixel. This is the most * memory-demanding step in the streakflat algorithm, since data from * _all_ images must reside in maim memory simultaneously. The approach * adopted here tries to balance memory and I/O requirements by keeping * in memory not the full images, but a chunk of blkSize image lines * instead. */ /* Alloc buffers holding blkSize-line chunks of each image. * Buffers are physically allocated only for non-memory-residing * images. */ for (image = 0; image < ioc->nimage; image++) { if (ioc->image[image].accessMode == Ddisk) { if (allocArray (&(block[image]), ioc->x_size, ioc->blkSize)) { g_error (malloc_err); return (1); } } else { block[image].nx = ioc->x_size; block[image].ny = ioc->blkSize; block[image].bufsize = (long)ioc->x_size * (long)ioc->blkSize; } } /* Loop over blocks. */ k = 0L; for (j = 1; j < ioc->y_size; j += ioc->blkSize) { /* This takes care of (eventually smaller) last block. */ remain = ioc->y_size - j + 1; blocksize = (remain < ioc->blkSize) ? remain : ioc->blkSize; /* Loop over image list and read current block from each image. */ for (image = 0; image < ioc->nimage; image++) { if (g_openImage (ioc, image, tmp_file)) return (1); g_getBlock (ioc, j, blocksize, &pix); /* If disk-resident, physically copy IMIO buffer into block. */ if (ioc->image[image].accessMode == Ddisk) memcpy (block[image].data , pix.data, (unsigned int)pix.bufsize*sizeof(float)); /* If memory-resident, just set the pointer. */ else block[image].data = pix.data; g_closeImage (ioc); /* This takes care of last block. */ block[image].ny = pix.ny; block[image].bufsize = pix.bufsize; } /* Loop over block pixels. */ for (l = 0; l < block[0].bufsize; l++) { /* n counts how many meaningful pixels across the image set.*/ n = 1; /* Loop over image set, computing for the current pixel the * ratio of image pixel by image mean ratio. Store result in * temporary vector for subsequent sorting. * * Note that we uset the temp array from its second element * on. The reason is that g_select is coded following the * 1-indexing method of Numerical Recipes. One could think * that by just passing temp-1 to the function would suffice. * We found in practice that this works with the Sun compiler * but fails with Compaq's. */ for (image = 0; image < ioc->nimage; image++) { if (block[image].data[l] != BADVAL) { temp[n] = block[image].data[l] / meanRatio[image]; n++; } } /* Now compute median for current pixel. */ if ((n == 0) || (n < ngood)) med = BADVAL; else if (n == 1) med = temp[1]; else { /* THIS NOT-SO-RIGOROUS MEDIAN COMPUTATION IS LEFT IN PLACE FOR COMPARISON WITH THE OLD STREAKFLAT TASK. THE CORRECT METHOD FOR SMALL NUMBER OF POINTS IS TO COMPUTE THE AVERAGE OF THE TWO MIDPOINTS WHEN N IS EVEN. */ /* if (((int)fmod((double)n,2.0)) || (n > 15)) */ med = g_select ((long)(n+1)/2, (long)n, temp); /* else { med = (g_select ((long)(n+1)/2, (long)n, temp) + g_select ((long)((n+1)/2-1), (long)n, temp)) * 0.5; } */ } /* Store median at corresponding pixel in output array. */ median->data[k++] = med; } } /* De-alloc block buffers. */ for (image = 0; image < ioc->nimage; image++) { if (ioc->image[image].accessMode == Ddisk) freeArray (&(block[image])); else { block[image].nx = 0; block[image].ny = 0; block[image].bufsize = 0L; } } /* Successfull return. */ return (0); }