1 /*
2 The following code is based on algorithms written by Richard White at STScI and made
3 available for use in CFITSIO in July 1999 and updated in January 2008.
4 */
5
6 # include <stdio.h>
7 # include <stdlib.h>
8 # include <math.h>
9 # include <limits.h>
10 # include <float.h>
11
12 #include "fitsio2.h"
13
14 /* nearest integer function */
15 # define NINT(x) ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5))
16
17 #define NULL_VALUE -2147483647 /* value used to represent undefined pixels */
18 #define ZERO_VALUE -2147483646 /* value used to represent zero-valued pixels */
19 #define N_RESERVED_VALUES 10 /* number of reserved values, starting with */
20 /* and including NULL_VALUE. These values */
21 /* may not be used to represent the quantized */
22 /* and scaled floating point pixel values */
23 /* If lossy Hcompression is used, and the */
24 /* array contains null values, then it is also */
25 /* possible for the compressed values to slightly */
26 /* exceed the range of the actual (lossless) values */
27 /* so we must reserve a little more space */
28
29 /* more than this many standard deviations from the mean is an outlier */
30 # define SIGMA_CLIP 5.
31 # define NITER 3 /* number of sigma-clipping iterations */
32
33 static int FnMeanSigma_short(short *array, long npix, int nullcheck,
34 short nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);
35 static int FnMeanSigma_int(int *array, long npix, int nullcheck,
36 int nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);
37 static int FnMeanSigma_float(float *array, long npix, int nullcheck,
38 float nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);
39 static int FnMeanSigma_double(double *array, long npix, int nullcheck,
40 double nullvalue, long *ngoodpix, double *mean, double *sigma, int *status);
41
42 static int FnNoise5_short(short *array, long nx, long ny, int nullcheck,
43 short nullvalue, long *ngood, short *minval, short *maxval,
44 double *n2, double *n3, double *n5, int *status);
45 static int FnNoise5_int(int *array, long nx, long ny, int nullcheck,
46 int nullvalue, long *ngood, int *minval, int *maxval,
47 double *n2, double *n3, double *n5, int *status);
48 static int FnNoise5_float(float *array, long nx, long ny, int nullcheck,
49 float nullvalue, long *ngood, float *minval, float *maxval,
50 double *n2, double *n3, double *n5, int *status);
51 static int FnNoise5_double(double *array, long nx, long ny, int nullcheck,
52 double nullvalue, long *ngood, double *minval, double *maxval,
53 double *n2, double *n3, double *n5, int *status);
54
55 static int FnNoise3_short(short *array, long nx, long ny, int nullcheck,
56 short nullvalue, long *ngood, short *minval, short *maxval, double *noise, int *status);
57 static int FnNoise3_int(int *array, long nx, long ny, int nullcheck,
58 int nullvalue, long *ngood, int *minval, int *maxval, double *noise, int *status);
59 static int FnNoise3_float(float *array, long nx, long ny, int nullcheck,
60 float nullvalue, long *ngood, float *minval, float *maxval, double *noise, int *status);
61 static int FnNoise3_double(double *array, long nx, long ny, int nullcheck,
62 double nullvalue, long *ngood, double *minval, double *maxval, double *noise, int *status);
63
64 static int FnNoise1_short(short *array, long nx, long ny,
65 int nullcheck, short nullvalue, double *noise, int *status);
66 static int FnNoise1_int(int *array, long nx, long ny,
67 int nullcheck, int nullvalue, double *noise, int *status);
68 static int FnNoise1_float(float *array, long nx, long ny,
69 int nullcheck, float nullvalue, double *noise, int *status);
70 static int FnNoise1_double(double *array, long nx, long ny,
71 int nullcheck, double nullvalue, double *noise, int *status);
72
73 static int FnCompare_short (const void *, const void *);
74 static int FnCompare_int (const void *, const void *);
75 static int FnCompare_float (const void *, const void *);
76 static int FnCompare_double (const void *, const void *);
77 static float quick_select_float(float arr[], int n);
78 static short quick_select_short(short arr[], int n);
79 static int quick_select_int(int arr[], int n);
80 static LONGLONG quick_select_longlong(LONGLONG arr[], int n);
81 static double quick_select_double(double arr[], int n);
82
83 /*---------------------------------------------------------------------------*/
fits_quantize_float(long row,float fdata[],long nxpix,long nypix,int nullcheck,float in_null_value,float qlevel,int dither_method,int idata[],double * bscale,double * bzero,int * iminval,int * imaxval)84 int fits_quantize_float (long row, float fdata[], long nxpix, long nypix, int nullcheck,
85 float in_null_value, float qlevel, int dither_method, int idata[], double *bscale,
86 double *bzero, int *iminval, int *imaxval) {
87
88 /* arguments:
89 long row i: if positive, used to calculate random dithering seed value
90 (this is only used when dithering the quantized values)
91 float fdata[] i: array of image pixels to be compressed
92 long nxpix i: number of pixels in each row of fdata
93 long nypix i: number of rows in fdata
94 nullcheck i: check for nullvalues in fdata?
95 float in_null_value i: value used to represent undefined pixels in fdata
96 float qlevel i: quantization level
97 int dither_method i; which dithering method to use
98 int idata[] o: values of fdata after applying bzero and bscale
99 double bscale o: scale factor
100 double bzero o: zero offset
101 int iminval o: minimum quantized value that is returned
102 int imaxval o: maximum quantized value that is returned
103
104 The function value will be one if the input fdata were copied to idata;
105 in this case the parameters bscale and bzero can be used to convert back to
106 nearly the original floating point values: fdata ~= idata * bscale + bzero.
107 If the function value is zero, the data were not copied to idata.
108 */
109
110 int status, iseed = 0;
111 long i, nx, ngood = 0;
112 double stdev, noise2, noise3, noise5; /* MAD 2nd, 3rd, and 5th order noise values */
113 float minval = 0., maxval = 0.; /* min & max of fdata */
114 double delta; /* bscale, 1 in idata = delta in fdata */
115 double zeropt; /* bzero */
116 double temp;
117 int nextrand = 0;
118 extern float *fits_rand_value; /* this is defined in imcompress.c */
119 LONGLONG iqfactor;
120
121 nx = nxpix * nypix;
122 if (nx <= 1) {
123 *bscale = 1.;
124 *bzero = 0.;
125 return (0);
126 }
127
128 if (qlevel >= 0.) {
129
130 /* estimate background noise using MAD pixel differences */
131 FnNoise5_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
132 &minval, &maxval, &noise2, &noise3, &noise5, &status);
133
134 if (nullcheck && ngood == 0) { /* special case of an image filled with Nulls */
135 /* set parameters to dummy values, which are not used */
136 minval = 0.;
137 maxval = 1.;
138 stdev = 1;
139 } else {
140
141 /* use the minimum of noise2, noise3, and noise5 as the best noise value */
142 stdev = noise3;
143 if (noise2 != 0. && noise2 < stdev) stdev = noise2;
144 if (noise5 != 0. && noise5 < stdev) stdev = noise5;
145 }
146
147 if (qlevel == 0.)
148 delta = stdev / 4.; /* default quantization */
149 else
150 delta = stdev / qlevel;
151
152 if (delta == 0.)
153 return (0); /* don't quantize */
154
155 } else {
156 /* negative value represents the absolute quantization level */
157 delta = -qlevel;
158
159 /* only nned to calculate the min and max values */
160 FnNoise3_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
161 &minval, &maxval, 0, &status);
162 }
163
164 /* check that the range of quantized levels is not > range of int */
165 if ((maxval - minval) / delta > 2. * 2147483647. - N_RESERVED_VALUES )
166 return (0); /* don't quantize */
167
168 if (row > 0) { /* we need to dither the quantized values */
169 if (!fits_rand_value)
170 if (fits_init_randoms()) return(MEMORY_ALLOCATION);
171
172 /* initialize the index to the next random number in the list */
173 iseed = (int) ((row - 1) % N_RANDOM);
174 nextrand = (int) (fits_rand_value[iseed] * 500.);
175 }
176
177 if (ngood == nx) { /* don't have to check for nulls */
178 /* return all positive values, if possible since some */
179 /* compression algorithms either only work for positive integers, */
180 /* or are more efficient. */
181
182 if (dither_method == SUBTRACTIVE_DITHER_2)
183 {
184 /* shift the range to be close to the value used to represent zeros */
185 zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
186 }
187 else if ((maxval - minval) / delta < 2147483647. - N_RESERVED_VALUES )
188 {
189 zeropt = minval;
190 /* fudge the zero point so it is an integer multiple of delta */
191 /* This helps to ensure the same scaling will be performed if the */
192 /* file undergoes multiple fpack/funpack cycles */
193 iqfactor = (LONGLONG) (zeropt/delta + 0.5);
194 zeropt = iqfactor * delta;
195 }
196 else
197 {
198 /* center the quantized levels around zero */
199 zeropt = (minval + maxval) / 2.;
200 }
201
202 if (row > 0) { /* dither the values when quantizing */
203 for (i = 0; i < nx; i++) {
204
205 if (dither_method == SUBTRACTIVE_DITHER_2 && fdata[i] == 0.0) {
206 idata[i] = ZERO_VALUE;
207 } else {
208 idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
209 }
210
211 nextrand++;
212 if (nextrand == N_RANDOM) {
213 iseed++;
214 if (iseed == N_RANDOM) iseed = 0;
215 nextrand = (int) (fits_rand_value[iseed] * 500);
216 }
217 }
218 } else { /* do not dither the values */
219
220 for (i = 0; i < nx; i++) {
221 idata[i] = NINT ((fdata[i] - zeropt) / delta);
222 }
223 }
224 }
225 else {
226 /* data contains null values; shift the range to be */
227 /* close to the value used to represent null values */
228 zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
229
230 if (row > 0) { /* dither the values */
231 for (i = 0; i < nx; i++) {
232 if (fdata[i] != in_null_value) {
233 if (dither_method == SUBTRACTIVE_DITHER_2 && fdata[i] == 0.0) {
234 idata[i] = ZERO_VALUE;
235 } else {
236 idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
237 }
238 } else {
239 idata[i] = NULL_VALUE;
240 }
241
242 /* increment the random number index, regardless */
243 nextrand++;
244 if (nextrand == N_RANDOM) {
245 iseed++;
246 if (iseed == N_RANDOM) iseed = 0;
247 nextrand = (int) (fits_rand_value[iseed] * 500);
248 }
249 }
250 } else { /* do not dither the values */
251 for (i = 0; i < nx; i++) {
252
253 if (fdata[i] != in_null_value) {
254 idata[i] = NINT((fdata[i] - zeropt) / delta);
255 } else {
256 idata[i] = NULL_VALUE;
257 }
258 }
259 }
260 }
261
262 /* calc min and max values */
263 temp = (minval - zeropt) / delta;
264 *iminval = NINT (temp);
265 temp = (maxval - zeropt) / delta;
266 *imaxval = NINT (temp);
267
268 *bscale = delta;
269 *bzero = zeropt;
270 return (1); /* yes, data have been quantized */
271 }
272 /*---------------------------------------------------------------------------*/
fits_quantize_double(long row,double fdata[],long nxpix,long nypix,int nullcheck,double in_null_value,float qlevel,int dither_method,int idata[],double * bscale,double * bzero,int * iminval,int * imaxval)273 int fits_quantize_double (long row, double fdata[], long nxpix, long nypix, int nullcheck,
274 double in_null_value, float qlevel, int dither_method, int idata[], double *bscale,
275 double *bzero, int *iminval, int *imaxval) {
276
277 /* arguments:
278 long row i: tile number = row number in the binary table
279 (this is only used when dithering the quantized values)
280 double fdata[] i: array of image pixels to be compressed
281 long nxpix i: number of pixels in each row of fdata
282 long nypix i: number of rows in fdata
283 nullcheck i: check for nullvalues in fdata?
284 double in_null_value i: value used to represent undefined pixels in fdata
285 float qlevel i: quantization level
286 int dither_method i; which dithering method to use
287 int idata[] o: values of fdata after applying bzero and bscale
288 double bscale o: scale factor
289 double bzero o: zero offset
290 int iminval o: minimum quantized value that is returned
291 int imaxval o: maximum quantized value that is returned
292
293 The function value will be one if the input fdata were copied to idata;
294 in this case the parameters bscale and bzero can be used to convert back to
295 nearly the original floating point values: fdata ~= idata * bscale + bzero.
296 If the function value is zero, the data were not copied to idata.
297 */
298
299 int status, iseed = 0;
300 long i, nx, ngood = 0;
301 double stdev, noise2 = 0., noise3 = 0., noise5 = 0.; /* MAD 2nd, 3rd, and 5th order noise values */
302 double minval = 0., maxval = 0.; /* min & max of fdata */
303 double delta; /* bscale, 1 in idata = delta in fdata */
304 double zeropt; /* bzero */
305 double temp;
306 int nextrand = 0;
307 extern float *fits_rand_value;
308 LONGLONG iqfactor;
309
310 nx = nxpix * nypix;
311 if (nx <= 1) {
312 *bscale = 1.;
313 *bzero = 0.;
314 return (0);
315 }
316
317 if (qlevel >= 0.) {
318
319 /* estimate background noise using MAD pixel differences */
320 FnNoise5_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
321 &minval, &maxval, &noise2, &noise3, &noise5, &status);
322
323 if (nullcheck && ngood == 0) { /* special case of an image filled with Nulls */
324 /* set parameters to dummy values, which are not used */
325 minval = 0.;
326 maxval = 1.;
327 stdev = 1;
328 } else {
329
330 /* use the minimum of noise2, noise3, and noise5 as the best noise value */
331 stdev = noise3;
332 if (noise2 != 0. && noise2 < stdev) stdev = noise2;
333 if (noise5 != 0. && noise5 < stdev) stdev = noise5;
334 }
335
336 if (qlevel == 0.)
337 delta = stdev / 4.; /* default quantization */
338 else
339 delta = stdev / qlevel;
340
341 if (delta == 0.)
342 return (0); /* don't quantize */
343
344 } else {
345 /* negative value represents the absolute quantization level */
346 delta = -qlevel;
347
348 /* only nned to calculate the min and max values */
349 FnNoise3_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood,
350 &minval, &maxval, 0, &status);
351 }
352
353 /* check that the range of quantized levels is not > range of int */
354 if ((maxval - minval) / delta > 2. * 2147483647. - N_RESERVED_VALUES )
355 return (0); /* don't quantize */
356
357 if (row > 0) { /* we need to dither the quantized values */
358 if (!fits_rand_value)
359 if (fits_init_randoms()) return(MEMORY_ALLOCATION);
360
361 /* initialize the index to the next random number in the list */
362 iseed = (int) ((row - 1) % N_RANDOM);
363 nextrand = (int) (fits_rand_value[iseed] * 500);
364 }
365
366 if (ngood == nx) { /* don't have to check for nulls */
367 /* return all positive values, if possible since some */
368 /* compression algorithms either only work for positive integers, */
369 /* or are more efficient. */
370
371 if (dither_method == SUBTRACTIVE_DITHER_2)
372 {
373 /* shift the range to be close to the value used to represent zeros */
374 zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
375 }
376 else if ((maxval - minval) / delta < 2147483647. - N_RESERVED_VALUES )
377 {
378 zeropt = minval;
379 /* fudge the zero point so it is an integer multiple of delta */
380 /* This helps to ensure the same scaling will be performed if the */
381 /* file undergoes multiple fpack/funpack cycles */
382 iqfactor = (LONGLONG) (zeropt/delta + 0.5);
383 zeropt = iqfactor * delta;
384 }
385 else
386 {
387 /* center the quantized levels around zero */
388 zeropt = (minval + maxval) / 2.;
389 }
390
391 if (row > 0) { /* dither the values when quantizing */
392 for (i = 0; i < nx; i++) {
393
394 if (dither_method == SUBTRACTIVE_DITHER_2 && fdata[i] == 0.0) {
395 idata[i] = ZERO_VALUE;
396 } else {
397 idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
398 }
399
400 nextrand++;
401 if (nextrand == N_RANDOM) {
402 iseed++;
403 if (iseed == N_RANDOM) iseed = 0;
404 nextrand = (int) (fits_rand_value[iseed] * 500);
405 }
406 }
407 } else { /* do not dither the values */
408
409 for (i = 0; i < nx; i++) {
410 idata[i] = NINT ((fdata[i] - zeropt) / delta);
411 }
412 }
413 }
414 else {
415 /* data contains null values; shift the range to be */
416 /* close to the value used to represent null values */
417 zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
418
419 if (row > 0) { /* dither the values */
420 for (i = 0; i < nx; i++) {
421 if (fdata[i] != in_null_value) {
422 if (dither_method == SUBTRACTIVE_DITHER_2 && fdata[i] == 0.0) {
423 idata[i] = ZERO_VALUE;
424 } else {
425 idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5);
426 }
427 } else {
428 idata[i] = NULL_VALUE;
429 }
430
431 /* increment the random number index, regardless */
432 nextrand++;
433 if (nextrand == N_RANDOM) {
434 iseed++;
435 if (iseed == N_RANDOM) iseed = 0;
436 nextrand = (int) (fits_rand_value[iseed] * 500);
437 }
438 }
439 } else { /* do not dither the values */
440 for (i = 0; i < nx; i++) {
441 if (fdata[i] != in_null_value)
442 idata[i] = NINT((fdata[i] - zeropt) / delta);
443 else
444 idata[i] = NULL_VALUE;
445 }
446 }
447 }
448
449 /* calc min and max values */
450 temp = (minval - zeropt) / delta;
451 *iminval = NINT (temp);
452 temp = (maxval - zeropt) / delta;
453 *imaxval = NINT (temp);
454
455 *bscale = delta;
456 *bzero = zeropt;
457
458 return (1); /* yes, data have been quantized */
459 }
460 /*--------------------------------------------------------------------------*/
fits_img_stats_short(short * array,long nx,long ny,int nullcheck,short nullvalue,long * ngoodpix,short * minvalue,short * maxvalue,double * mean,double * sigma,double * noise1,double * noise2,double * noise3,double * noise5,int * status)461 int fits_img_stats_short(short *array, /* 2 dimensional array of image pixels */
462 long nx, /* number of pixels in each row of the image */
463 long ny, /* number of rows in the image */
464 /* (if this is a 3D image, then ny should be the */
465 /* product of the no. of rows times the no. of planes) */
466 int nullcheck, /* check for null values, if true */
467 short nullvalue, /* value of null pixels, if nullcheck is true */
468
469 /* returned parameters (if the pointer is not null) */
470 long *ngoodpix, /* number of non-null pixels in the image */
471 short *minvalue, /* returned minimum non-null value in the array */
472 short *maxvalue, /* returned maximum non-null value in the array */
473 double *mean, /* returned mean value of all non-null pixels */
474 double *sigma, /* returned R.M.S. value of all non-null pixels */
475 double *noise1, /* 1st order estimate of noise in image background level */
476 double *noise2, /* 2nd order estimate of noise in image background level */
477 double *noise3, /* 3rd order estimate of noise in image background level */
478 double *noise5, /* 5th order estimate of noise in image background level */
479 int *status) /* error status */
480
481 /*
482 Compute statistics of the input short integer image.
483 */
484 {
485 long ngood;
486 short minval = 0, maxval = 0;
487 double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
488
489 /* need to calculate mean and/or sigma and/or limits? */
490 if (mean || sigma ) {
491 FnMeanSigma_short(array, nx * ny, nullcheck, nullvalue,
492 &ngood, &xmean, &xsigma, status);
493
494 if (ngoodpix) *ngoodpix = ngood;
495 if (mean) *mean = xmean;
496 if (sigma) *sigma = xsigma;
497 }
498
499 if (noise1) {
500 FnNoise1_short(array, nx, ny, nullcheck, nullvalue,
501 &xnoise, status);
502
503 *noise1 = xnoise;
504 }
505
506 if (minvalue || maxvalue || noise3) {
507 FnNoise5_short(array, nx, ny, nullcheck, nullvalue,
508 &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
509
510 if (ngoodpix) *ngoodpix = ngood;
511 if (minvalue) *minvalue= minval;
512 if (maxvalue) *maxvalue = maxval;
513 if (noise2) *noise2 = xnoise2;
514 if (noise3) *noise3 = xnoise3;
515 if (noise5) *noise5 = xnoise5;
516 }
517 return(*status);
518 }
519 /*--------------------------------------------------------------------------*/
fits_img_stats_int(int * array,long nx,long ny,int nullcheck,int nullvalue,long * ngoodpix,int * minvalue,int * maxvalue,double * mean,double * sigma,double * noise1,double * noise2,double * noise3,double * noise5,int * status)520 int fits_img_stats_int(int *array, /* 2 dimensional array of image pixels */
521 long nx, /* number of pixels in each row of the image */
522 long ny, /* number of rows in the image */
523 /* (if this is a 3D image, then ny should be the */
524 /* product of the no. of rows times the no. of planes) */
525 int nullcheck, /* check for null values, if true */
526 int nullvalue, /* value of null pixels, if nullcheck is true */
527
528 /* returned parameters (if the pointer is not null) */
529 long *ngoodpix, /* number of non-null pixels in the image */
530 int *minvalue, /* returned minimum non-null value in the array */
531 int *maxvalue, /* returned maximum non-null value in the array */
532 double *mean, /* returned mean value of all non-null pixels */
533 double *sigma, /* returned R.M.S. value of all non-null pixels */
534 double *noise1, /* 1st order estimate of noise in image background level */
535 double *noise2, /* 2nd order estimate of noise in image background level */
536 double *noise3, /* 3rd order estimate of noise in image background level */
537 double *noise5, /* 5th order estimate of noise in image background level */
538 int *status) /* error status */
539
540 /*
541 Compute statistics of the input integer image.
542 */
543 {
544 long ngood;
545 int minval = 0, maxval = 0;
546 double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
547
548 /* need to calculate mean and/or sigma and/or limits? */
549 if (mean || sigma ) {
550 FnMeanSigma_int(array, nx * ny, nullcheck, nullvalue,
551 &ngood, &xmean, &xsigma, status);
552
553 if (ngoodpix) *ngoodpix = ngood;
554 if (mean) *mean = xmean;
555 if (sigma) *sigma = xsigma;
556 }
557
558 if (noise1) {
559 FnNoise1_int(array, nx, ny, nullcheck, nullvalue,
560 &xnoise, status);
561
562 *noise1 = xnoise;
563 }
564
565 if (minvalue || maxvalue || noise3) {
566 FnNoise5_int(array, nx, ny, nullcheck, nullvalue,
567 &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
568
569 if (ngoodpix) *ngoodpix = ngood;
570 if (minvalue) *minvalue= minval;
571 if (maxvalue) *maxvalue = maxval;
572 if (noise2) *noise2 = xnoise2;
573 if (noise3) *noise3 = xnoise3;
574 if (noise5) *noise5 = xnoise5;
575 }
576 return(*status);
577 }
578 /*--------------------------------------------------------------------------*/
fits_img_stats_float(float * array,long nx,long ny,int nullcheck,float nullvalue,long * ngoodpix,float * minvalue,float * maxvalue,double * mean,double * sigma,double * noise1,double * noise2,double * noise3,double * noise5,int * status)579 int fits_img_stats_float(float *array, /* 2 dimensional array of image pixels */
580 long nx, /* number of pixels in each row of the image */
581 long ny, /* number of rows in the image */
582 /* (if this is a 3D image, then ny should be the */
583 /* product of the no. of rows times the no. of planes) */
584 int nullcheck, /* check for null values, if true */
585 float nullvalue, /* value of null pixels, if nullcheck is true */
586
587 /* returned parameters (if the pointer is not null) */
588 long *ngoodpix, /* number of non-null pixels in the image */
589 float *minvalue, /* returned minimum non-null value in the array */
590 float *maxvalue, /* returned maximum non-null value in the array */
591 double *mean, /* returned mean value of all non-null pixels */
592 double *sigma, /* returned R.M.S. value of all non-null pixels */
593 double *noise1, /* 1st order estimate of noise in image background level */
594 double *noise2, /* 2nd order estimate of noise in image background level */
595 double *noise3, /* 3rd order estimate of noise in image background level */
596 double *noise5, /* 5th order estimate of noise in image background level */
597 int *status) /* error status */
598
599 /*
600 Compute statistics of the input float image.
601 */
602 {
603 long ngood;
604 float minval, maxval;
605 double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.;
606
607 /* need to calculate mean and/or sigma and/or limits? */
608 if (mean || sigma ) {
609 FnMeanSigma_float(array, nx * ny, nullcheck, nullvalue,
610 &ngood, &xmean, &xsigma, status);
611
612 if (ngoodpix) *ngoodpix = ngood;
613 if (mean) *mean = xmean;
614 if (sigma) *sigma = xsigma;
615 }
616
617 if (noise1) {
618 FnNoise1_float(array, nx, ny, nullcheck, nullvalue,
619 &xnoise, status);
620
621 *noise1 = xnoise;
622 }
623
624 if (minvalue || maxvalue || noise3) {
625 FnNoise5_float(array, nx, ny, nullcheck, nullvalue,
626 &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status);
627
628 if (ngoodpix) *ngoodpix = ngood;
629 if (minvalue) *minvalue= minval;
630 if (maxvalue) *maxvalue = maxval;
631 if (noise2) *noise2 = xnoise2;
632 if (noise3) *noise3 = xnoise3;
633 if (noise5) *noise5 = xnoise5;
634 }
635 return(*status);
636 }
637 /*--------------------------------------------------------------------------*/
FnMeanSigma_short(short * array,long npix,int nullcheck,short nullvalue,long * ngoodpix,double * mean,double * sigma,int * status)638 static int FnMeanSigma_short
639 (short *array, /* 2 dimensional array of image pixels */
640 long npix, /* number of pixels in the image */
641 int nullcheck, /* check for null values, if true */
642 short nullvalue, /* value of null pixels, if nullcheck is true */
643
644 /* returned parameters */
645
646 long *ngoodpix, /* number of non-null pixels in the image */
647 double *mean, /* returned mean value of all non-null pixels */
648 double *sigma, /* returned R.M.S. value of all non-null pixels */
649 int *status) /* error status */
650
651 /*
652 Compute mean and RMS sigma of the non-null pixels in the input array.
653 */
654 {
655 long ii, ngood = 0;
656 short *value;
657 double sum = 0., sum2 = 0., xtemp;
658
659 value = array;
660
661 if (nullcheck) {
662 for (ii = 0; ii < npix; ii++, value++) {
663 if (*value != nullvalue) {
664 ngood++;
665 xtemp = (double) *value;
666 sum += xtemp;
667 sum2 += (xtemp * xtemp);
668 }
669 }
670 } else {
671 ngood = npix;
672 for (ii = 0; ii < npix; ii++, value++) {
673 xtemp = (double) *value;
674 sum += xtemp;
675 sum2 += (xtemp * xtemp);
676 }
677 }
678
679 if (ngood > 1) {
680 if (ngoodpix) *ngoodpix = ngood;
681 xtemp = sum / ngood;
682 if (mean) *mean = xtemp;
683 if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
684 } else if (ngood == 1){
685 if (ngoodpix) *ngoodpix = 1;
686 if (mean) *mean = sum;
687 if (sigma) *sigma = 0.0;
688 } else {
689 if (ngoodpix) *ngoodpix = 0;
690 if (mean) *mean = 0.;
691 if (sigma) *sigma = 0.;
692 }
693 return(*status);
694 }
695 /*--------------------------------------------------------------------------*/
FnMeanSigma_int(int * array,long npix,int nullcheck,int nullvalue,long * ngoodpix,double * mean,double * sigma,int * status)696 static int FnMeanSigma_int
697 (int *array, /* 2 dimensional array of image pixels */
698 long npix, /* number of pixels in the image */
699 int nullcheck, /* check for null values, if true */
700 int nullvalue, /* value of null pixels, if nullcheck is true */
701
702 /* returned parameters */
703
704 long *ngoodpix, /* number of non-null pixels in the image */
705 double *mean, /* returned mean value of all non-null pixels */
706 double *sigma, /* returned R.M.S. value of all non-null pixels */
707 int *status) /* error status */
708
709 /*
710 Compute mean and RMS sigma of the non-null pixels in the input array.
711 */
712 {
713 long ii, ngood = 0;
714 int *value;
715 double sum = 0., sum2 = 0., xtemp;
716
717 value = array;
718
719 if (nullcheck) {
720 for (ii = 0; ii < npix; ii++, value++) {
721 if (*value != nullvalue) {
722 ngood++;
723 xtemp = (double) *value;
724 sum += xtemp;
725 sum2 += (xtemp * xtemp);
726 }
727 }
728 } else {
729 ngood = npix;
730 for (ii = 0; ii < npix; ii++, value++) {
731 xtemp = (double) *value;
732 sum += xtemp;
733 sum2 += (xtemp * xtemp);
734 }
735 }
736
737 if (ngood > 1) {
738 if (ngoodpix) *ngoodpix = ngood;
739 xtemp = sum / ngood;
740 if (mean) *mean = xtemp;
741 if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
742 } else if (ngood == 1){
743 if (ngoodpix) *ngoodpix = 1;
744 if (mean) *mean = sum;
745 if (sigma) *sigma = 0.0;
746 } else {
747 if (ngoodpix) *ngoodpix = 0;
748 if (mean) *mean = 0.;
749 if (sigma) *sigma = 0.;
750 }
751 return(*status);
752 }
753 /*--------------------------------------------------------------------------*/
FnMeanSigma_float(float * array,long npix,int nullcheck,float nullvalue,long * ngoodpix,double * mean,double * sigma,int * status)754 static int FnMeanSigma_float
755 (float *array, /* 2 dimensional array of image pixels */
756 long npix, /* number of pixels in the image */
757 int nullcheck, /* check for null values, if true */
758 float nullvalue, /* value of null pixels, if nullcheck is true */
759
760 /* returned parameters */
761
762 long *ngoodpix, /* number of non-null pixels in the image */
763 double *mean, /* returned mean value of all non-null pixels */
764 double *sigma, /* returned R.M.S. value of all non-null pixels */
765 int *status) /* error status */
766
767 /*
768 Compute mean and RMS sigma of the non-null pixels in the input array.
769 */
770 {
771 long ii, ngood = 0;
772 float *value;
773 double sum = 0., sum2 = 0., xtemp;
774
775 value = array;
776
777 if (nullcheck) {
778 for (ii = 0; ii < npix; ii++, value++) {
779 if (*value != nullvalue) {
780 ngood++;
781 xtemp = (double) *value;
782 sum += xtemp;
783 sum2 += (xtemp * xtemp);
784 }
785 }
786 } else {
787 ngood = npix;
788 for (ii = 0; ii < npix; ii++, value++) {
789 xtemp = (double) *value;
790 sum += xtemp;
791 sum2 += (xtemp * xtemp);
792 }
793 }
794
795 if (ngood > 1) {
796 if (ngoodpix) *ngoodpix = ngood;
797 xtemp = sum / ngood;
798 if (mean) *mean = xtemp;
799 if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
800 } else if (ngood == 1){
801 if (ngoodpix) *ngoodpix = 1;
802 if (mean) *mean = sum;
803 if (sigma) *sigma = 0.0;
804 } else {
805 if (ngoodpix) *ngoodpix = 0;
806 if (mean) *mean = 0.;
807 if (sigma) *sigma = 0.;
808 }
809 return(*status);
810 }
811 /*--------------------------------------------------------------------------*/
FnMeanSigma_double(double * array,long npix,int nullcheck,double nullvalue,long * ngoodpix,double * mean,double * sigma,int * status)812 static int FnMeanSigma_double
813 (double *array, /* 2 dimensional array of image pixels */
814 long npix, /* number of pixels in the image */
815 int nullcheck, /* check for null values, if true */
816 double nullvalue, /* value of null pixels, if nullcheck is true */
817
818 /* returned parameters */
819
820 long *ngoodpix, /* number of non-null pixels in the image */
821 double *mean, /* returned mean value of all non-null pixels */
822 double *sigma, /* returned R.M.S. value of all non-null pixels */
823 int *status) /* error status */
824
825 /*
826 Compute mean and RMS sigma of the non-null pixels in the input array.
827 */
828 {
829 long ii, ngood = 0;
830 double *value;
831 double sum = 0., sum2 = 0., xtemp;
832
833 value = array;
834
835 if (nullcheck) {
836 for (ii = 0; ii < npix; ii++, value++) {
837 if (*value != nullvalue) {
838 ngood++;
839 xtemp = *value;
840 sum += xtemp;
841 sum2 += (xtemp * xtemp);
842 }
843 }
844 } else {
845 ngood = npix;
846 for (ii = 0; ii < npix; ii++, value++) {
847 xtemp = *value;
848 sum += xtemp;
849 sum2 += (xtemp * xtemp);
850 }
851 }
852
853 if (ngood > 1) {
854 if (ngoodpix) *ngoodpix = ngood;
855 xtemp = sum / ngood;
856 if (mean) *mean = xtemp;
857 if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp));
858 } else if (ngood == 1){
859 if (ngoodpix) *ngoodpix = 1;
860 if (mean) *mean = sum;
861 if (sigma) *sigma = 0.0;
862 } else {
863 if (ngoodpix) *ngoodpix = 0;
864 if (mean) *mean = 0.;
865 if (sigma) *sigma = 0.;
866 }
867 return(*status);
868 }
869 /*--------------------------------------------------------------------------*/
FnNoise5_short(short * array,long nx,long ny,int nullcheck,short nullvalue,long * ngood,short * minval,short * maxval,double * noise2,double * noise3,double * noise5,int * status)870 static int FnNoise5_short
871 (short *array, /* 2 dimensional array of image pixels */
872 long nx, /* number of pixels in each row of the image */
873 long ny, /* number of rows in the image */
874 int nullcheck, /* check for null values, if true */
875 short nullvalue, /* value of null pixels, if nullcheck is true */
876 /* returned parameters */
877 long *ngood, /* number of good, non-null pixels? */
878 short *minval, /* minimum non-null value */
879 short *maxval, /* maximum non-null value */
880 double *noise2, /* returned 2nd order MAD of all non-null pixels */
881 double *noise3, /* returned 3rd order MAD of all non-null pixels */
882 double *noise5, /* returned 5th order MAD of all non-null pixels */
883 int *status) /* error status */
884
885 /*
886 Estimate the median and background noise in the input image using 2nd, 3rd and 5th
887 order Median Absolute Differences.
888
889 The noise in the background of the image is calculated using the MAD algorithms
890 developed for deriving the signal to noise ratio in spectra
891 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
892
893 3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
894
895 The returned estimates are the median of the values that are computed for each
896 row of the image.
897 */
898 {
899 long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
900 int *differences2, *differences3, *differences5;
901 short *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
902 short xminval = SHRT_MAX, xmaxval = SHRT_MIN;
903 int do_range = 0;
904 double *diffs2, *diffs3, *diffs5;
905 double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
906
907 if (nx < 9) {
908 /* treat entire array as an image with a single row */
909 nx = nx * ny;
910 ny = 1;
911 }
912
913 /* rows must have at least 9 pixels */
914 if (nx < 9) {
915
916 for (ii = 0; ii < nx; ii++) {
917 if (nullcheck && array[ii] == nullvalue)
918 continue;
919 else {
920 if (array[ii] < xminval) xminval = array[ii];
921 if (array[ii] > xmaxval) xmaxval = array[ii];
922 ngoodpix++;
923 }
924 }
925 if (minval) *minval = xminval;
926 if (maxval) *maxval = xmaxval;
927 if (ngood) *ngood = ngoodpix;
928 if (noise2) *noise2 = 0.;
929 if (noise3) *noise3 = 0.;
930 if (noise5) *noise5 = 0.;
931 return(*status);
932 }
933
934 /* do we need to compute the min and max value? */
935 if (minval || maxval) do_range = 1;
936
937 /* allocate arrays used to compute the median and noise estimates */
938 differences2 = calloc(nx, sizeof(int));
939 if (!differences2) {
940 *status = MEMORY_ALLOCATION;
941 return(*status);
942 }
943 differences3 = calloc(nx, sizeof(int));
944 if (!differences3) {
945 free(differences2);
946 *status = MEMORY_ALLOCATION;
947 return(*status);
948 }
949 differences5 = calloc(nx, sizeof(int));
950 if (!differences5) {
951 free(differences2);
952 free(differences3);
953 *status = MEMORY_ALLOCATION;
954 return(*status);
955 }
956
957 diffs2 = calloc(ny, sizeof(double));
958 if (!diffs2) {
959 free(differences2);
960 free(differences3);
961 free(differences5);
962 *status = MEMORY_ALLOCATION;
963 return(*status);
964 }
965
966 diffs3 = calloc(ny, sizeof(double));
967 if (!diffs3) {
968 free(differences2);
969 free(differences3);
970 free(differences5);
971 free(diffs2);
972 *status = MEMORY_ALLOCATION;
973 return(*status);
974 }
975
976 diffs5 = calloc(ny, sizeof(double));
977 if (!diffs5) {
978 free(differences2);
979 free(differences3);
980 free(differences5);
981 free(diffs2);
982 free(diffs3);
983 *status = MEMORY_ALLOCATION;
984 return(*status);
985 }
986
987 /* loop over each row of the image */
988 for (jj=0; jj < ny; jj++) {
989
990 rowpix = array + (jj * nx); /* point to first pixel in the row */
991
992 /***** find the first valid pixel in row */
993 ii = 0;
994 if (nullcheck)
995 while (ii < nx && rowpix[ii] == nullvalue) ii++;
996
997 if (ii == nx) continue; /* hit end of row */
998 v1 = rowpix[ii]; /* store the good pixel value */
999 ngoodpix++;
1000
1001 if (do_range) {
1002 if (v1 < xminval) xminval = v1;
1003 if (v1 > xmaxval) xmaxval = v1;
1004 }
1005
1006 /***** find the 2nd valid pixel in row (which we will skip over) */
1007 ii++;
1008 if (nullcheck)
1009 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1010
1011 if (ii == nx) continue; /* hit end of row */
1012 v2 = rowpix[ii]; /* store the good pixel value */
1013 ngoodpix++;
1014
1015 if (do_range) {
1016 if (v2 < xminval) xminval = v2;
1017 if (v2 > xmaxval) xmaxval = v2;
1018 }
1019
1020 /***** find the 3rd valid pixel in row */
1021 ii++;
1022 if (nullcheck)
1023 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1024
1025 if (ii == nx) continue; /* hit end of row */
1026 v3 = rowpix[ii]; /* store the good pixel value */
1027 ngoodpix++;
1028
1029 if (do_range) {
1030 if (v3 < xminval) xminval = v3;
1031 if (v3 > xmaxval) xmaxval = v3;
1032 }
1033
1034 /* find the 4nd valid pixel in row (to be skipped) */
1035 ii++;
1036 if (nullcheck)
1037 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1038
1039 if (ii == nx) continue; /* hit end of row */
1040 v4 = rowpix[ii]; /* store the good pixel value */
1041 ngoodpix++;
1042
1043 if (do_range) {
1044 if (v4 < xminval) xminval = v4;
1045 if (v4 > xmaxval) xmaxval = v4;
1046 }
1047
1048 /* find the 5th valid pixel in row (to be skipped) */
1049 ii++;
1050 if (nullcheck)
1051 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1052
1053 if (ii == nx) continue; /* hit end of row */
1054 v5 = rowpix[ii]; /* store the good pixel value */
1055 ngoodpix++;
1056
1057 if (do_range) {
1058 if (v5 < xminval) xminval = v5;
1059 if (v5 > xmaxval) xmaxval = v5;
1060 }
1061
1062 /* find the 6th valid pixel in row (to be skipped) */
1063 ii++;
1064 if (nullcheck)
1065 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1066
1067 if (ii == nx) continue; /* hit end of row */
1068 v6 = rowpix[ii]; /* store the good pixel value */
1069 ngoodpix++;
1070
1071 if (do_range) {
1072 if (v6 < xminval) xminval = v6;
1073 if (v6 > xmaxval) xmaxval = v6;
1074 }
1075
1076 /* find the 7th valid pixel in row (to be skipped) */
1077 ii++;
1078 if (nullcheck)
1079 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1080
1081 if (ii == nx) continue; /* hit end of row */
1082 v7 = rowpix[ii]; /* store the good pixel value */
1083 ngoodpix++;
1084
1085 if (do_range) {
1086 if (v7 < xminval) xminval = v7;
1087 if (v7 > xmaxval) xmaxval = v7;
1088 }
1089
1090 /* find the 8th valid pixel in row (to be skipped) */
1091 ii++;
1092 if (nullcheck)
1093 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1094
1095 if (ii == nx) continue; /* hit end of row */
1096 v8 = rowpix[ii]; /* store the good pixel value */
1097 ngoodpix++;
1098
1099 if (do_range) {
1100 if (v8 < xminval) xminval = v8;
1101 if (v8 > xmaxval) xmaxval = v8;
1102 }
1103 /* now populate the differences arrays */
1104 /* for the remaining pixels in the row */
1105 nvals = 0;
1106 nvals2 = 0;
1107 for (ii++; ii < nx; ii++) {
1108
1109 /* find the next valid pixel in row */
1110 if (nullcheck)
1111 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1112
1113 if (ii == nx) break; /* hit end of row */
1114 v9 = rowpix[ii]; /* store the good pixel value */
1115
1116 if (do_range) {
1117 if (v9 < xminval) xminval = v9;
1118 if (v9 > xmaxval) xmaxval = v9;
1119 }
1120
1121 /* construct array of absolute differences */
1122
1123 if (!(v5 == v6 && v6 == v7) ) {
1124 differences2[nvals2] = abs((int) v5 - (int) v7);
1125 nvals2++;
1126 }
1127
1128 if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
1129 differences3[nvals] = abs((2 * (int) v5) - (int) v3 - (int) v7);
1130 differences5[nvals] = abs((6 * (int) v5) - (4 * (int) v3) - (4 * (int) v7) + (int) v1 + (int) v9);
1131 nvals++;
1132 } else {
1133 /* ignore constant background regions */
1134 ngoodpix++;
1135 }
1136
1137 /* shift over 1 pixel */
1138 v1 = v2;
1139 v2 = v3;
1140 v3 = v4;
1141 v4 = v5;
1142 v5 = v6;
1143 v6 = v7;
1144 v7 = v8;
1145 v8 = v9;
1146 } /* end of loop over pixels in the row */
1147
1148 /* compute the median diffs */
1149 /* Note that there are 8 more pixel values than there are diffs values. */
1150 ngoodpix += nvals;
1151
1152 if (nvals == 0) {
1153 continue; /* cannot compute medians on this row */
1154 } else if (nvals == 1) {
1155 if (nvals2 == 1) {
1156 diffs2[nrows2] = differences2[0];
1157 nrows2++;
1158 }
1159
1160 diffs3[nrows] = differences3[0];
1161 diffs5[nrows] = differences5[0];
1162 } else {
1163 /* quick_select returns the median MUCH faster than using qsort */
1164 if (nvals2 > 1) {
1165 diffs2[nrows2] = quick_select_int(differences2, nvals);
1166 nrows2++;
1167 }
1168
1169 diffs3[nrows] = quick_select_int(differences3, nvals);
1170 diffs5[nrows] = quick_select_int(differences5, nvals);
1171 }
1172
1173 nrows++;
1174 } /* end of loop over rows */
1175
1176 /* compute median of the values for each row */
1177 if (nrows == 0) {
1178 xnoise3 = 0;
1179 xnoise5 = 0;
1180 } else if (nrows == 1) {
1181 xnoise3 = diffs3[0];
1182 xnoise5 = diffs5[0];
1183 } else {
1184 qsort(diffs3, nrows, sizeof(double), FnCompare_double);
1185 qsort(diffs5, nrows, sizeof(double), FnCompare_double);
1186 xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
1187 xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
1188 }
1189
1190 if (nrows2 == 0) {
1191 xnoise2 = 0;
1192 } else if (nrows2 == 1) {
1193 xnoise2 = diffs2[0];
1194 } else {
1195 qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
1196 xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
1197 }
1198
1199 if (ngood) *ngood = ngoodpix;
1200 if (minval) *minval = xminval;
1201 if (maxval) *maxval = xmaxval;
1202 if (noise2) *noise2 = 1.0483579 * xnoise2;
1203 if (noise3) *noise3 = 0.6052697 * xnoise3;
1204 if (noise5) *noise5 = 0.1772048 * xnoise5;
1205
1206 free(diffs5);
1207 free(diffs3);
1208 free(diffs2);
1209 free(differences5);
1210 free(differences3);
1211 free(differences2);
1212
1213 return(*status);
1214 }
1215 /*--------------------------------------------------------------------------*/
FnNoise5_int(int * array,long nx,long ny,int nullcheck,int nullvalue,long * ngood,int * minval,int * maxval,double * noise2,double * noise3,double * noise5,int * status)1216 static int FnNoise5_int
1217 (int *array, /* 2 dimensional array of image pixels */
1218 long nx, /* number of pixels in each row of the image */
1219 long ny, /* number of rows in the image */
1220 int nullcheck, /* check for null values, if true */
1221 int nullvalue, /* value of null pixels, if nullcheck is true */
1222 /* returned parameters */
1223 long *ngood, /* number of good, non-null pixels? */
1224 int *minval, /* minimum non-null value */
1225 int *maxval, /* maximum non-null value */
1226 double *noise2, /* returned 2nd order MAD of all non-null pixels */
1227 double *noise3, /* returned 3rd order MAD of all non-null pixels */
1228 double *noise5, /* returned 5th order MAD of all non-null pixels */
1229 int *status) /* error status */
1230
1231 /*
1232 Estimate the median and background noise in the input image using 2nd, 3rd and 5th
1233 order Median Absolute Differences.
1234
1235 The noise in the background of the image is calculated using the MAD algorithms
1236 developed for deriving the signal to noise ratio in spectra
1237 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
1238
1239 3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
1240
1241 The returned estimates are the median of the values that are computed for each
1242 row of the image.
1243 */
1244 {
1245 long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
1246 LONGLONG *differences2, *differences3, *differences5, tdiff;
1247 int *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
1248 int xminval = INT_MAX, xmaxval = INT_MIN;
1249 int do_range = 0;
1250 double *diffs2, *diffs3, *diffs5;
1251 double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
1252
1253 if (nx < 9) {
1254 /* treat entire array as an image with a single row */
1255 nx = nx * ny;
1256 ny = 1;
1257 }
1258
1259 /* rows must have at least 9 pixels */
1260 if (nx < 9) {
1261
1262 for (ii = 0; ii < nx; ii++) {
1263 if (nullcheck && array[ii] == nullvalue)
1264 continue;
1265 else {
1266 if (array[ii] < xminval) xminval = array[ii];
1267 if (array[ii] > xmaxval) xmaxval = array[ii];
1268 ngoodpix++;
1269 }
1270 }
1271 if (minval) *minval = xminval;
1272 if (maxval) *maxval = xmaxval;
1273 if (ngood) *ngood = ngoodpix;
1274 if (noise2) *noise2 = 0.;
1275 if (noise3) *noise3 = 0.;
1276 if (noise5) *noise5 = 0.;
1277 return(*status);
1278 }
1279
1280 /* do we need to compute the min and max value? */
1281 if (minval || maxval) do_range = 1;
1282
1283 /* allocate arrays used to compute the median and noise estimates */
1284 differences2 = calloc(nx, sizeof(LONGLONG));
1285 if (!differences2) {
1286 *status = MEMORY_ALLOCATION;
1287 return(*status);
1288 }
1289 differences3 = calloc(nx, sizeof(LONGLONG));
1290 if (!differences3) {
1291 free(differences2);
1292 *status = MEMORY_ALLOCATION;
1293 return(*status);
1294 }
1295 differences5 = calloc(nx, sizeof(LONGLONG));
1296 if (!differences5) {
1297 free(differences2);
1298 free(differences3);
1299 *status = MEMORY_ALLOCATION;
1300 return(*status);
1301 }
1302
1303 diffs2 = calloc(ny, sizeof(double));
1304 if (!diffs2) {
1305 free(differences2);
1306 free(differences3);
1307 free(differences5);
1308 *status = MEMORY_ALLOCATION;
1309 return(*status);
1310 }
1311
1312 diffs3 = calloc(ny, sizeof(double));
1313 if (!diffs3) {
1314 free(differences2);
1315 free(differences3);
1316 free(differences5);
1317 free(diffs2);
1318 *status = MEMORY_ALLOCATION;
1319 return(*status);
1320 }
1321
1322 diffs5 = calloc(ny, sizeof(double));
1323 if (!diffs5) {
1324 free(differences2);
1325 free(differences3);
1326 free(differences5);
1327 free(diffs2);
1328 free(diffs3);
1329 *status = MEMORY_ALLOCATION;
1330 return(*status);
1331 }
1332
1333 /* loop over each row of the image */
1334 for (jj=0; jj < ny; jj++) {
1335
1336 rowpix = array + (jj * nx); /* point to first pixel in the row */
1337
1338 /***** find the first valid pixel in row */
1339 ii = 0;
1340 if (nullcheck)
1341 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1342
1343 if (ii == nx) continue; /* hit end of row */
1344 v1 = rowpix[ii]; /* store the good pixel value */
1345 ngoodpix++;
1346
1347 if (do_range) {
1348 if (v1 < xminval) xminval = v1;
1349 if (v1 > xmaxval) xmaxval = v1;
1350 }
1351
1352 /***** find the 2nd valid pixel in row (which we will skip over) */
1353 ii++;
1354 if (nullcheck)
1355 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1356
1357 if (ii == nx) continue; /* hit end of row */
1358 v2 = rowpix[ii]; /* store the good pixel value */
1359 ngoodpix++;
1360
1361 if (do_range) {
1362 if (v2 < xminval) xminval = v2;
1363 if (v2 > xmaxval) xmaxval = v2;
1364 }
1365
1366 /***** find the 3rd valid pixel in row */
1367 ii++;
1368 if (nullcheck)
1369 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1370
1371 if (ii == nx) continue; /* hit end of row */
1372 v3 = rowpix[ii]; /* store the good pixel value */
1373 ngoodpix++;
1374
1375 if (do_range) {
1376 if (v3 < xminval) xminval = v3;
1377 if (v3 > xmaxval) xmaxval = v3;
1378 }
1379
1380 /* find the 4nd valid pixel in row (to be skipped) */
1381 ii++;
1382 if (nullcheck)
1383 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1384
1385 if (ii == nx) continue; /* hit end of row */
1386 v4 = rowpix[ii]; /* store the good pixel value */
1387 ngoodpix++;
1388
1389 if (do_range) {
1390 if (v4 < xminval) xminval = v4;
1391 if (v4 > xmaxval) xmaxval = v4;
1392 }
1393
1394 /* find the 5th valid pixel in row (to be skipped) */
1395 ii++;
1396 if (nullcheck)
1397 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1398
1399 if (ii == nx) continue; /* hit end of row */
1400 v5 = rowpix[ii]; /* store the good pixel value */
1401 ngoodpix++;
1402
1403 if (do_range) {
1404 if (v5 < xminval) xminval = v5;
1405 if (v5 > xmaxval) xmaxval = v5;
1406 }
1407
1408 /* find the 6th valid pixel in row (to be skipped) */
1409 ii++;
1410 if (nullcheck)
1411 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1412
1413 if (ii == nx) continue; /* hit end of row */
1414 v6 = rowpix[ii]; /* store the good pixel value */
1415 ngoodpix++;
1416
1417 if (do_range) {
1418 if (v6 < xminval) xminval = v6;
1419 if (v6 > xmaxval) xmaxval = v6;
1420 }
1421
1422 /* find the 7th valid pixel in row (to be skipped) */
1423 ii++;
1424 if (nullcheck)
1425 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1426
1427 if (ii == nx) continue; /* hit end of row */
1428 v7 = rowpix[ii]; /* store the good pixel value */
1429 ngoodpix++;
1430
1431 if (do_range) {
1432 if (v7 < xminval) xminval = v7;
1433 if (v7 > xmaxval) xmaxval = v7;
1434 }
1435
1436 /* find the 8th valid pixel in row (to be skipped) */
1437 ii++;
1438 if (nullcheck)
1439 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1440
1441 if (ii == nx) continue; /* hit end of row */
1442 v8 = rowpix[ii]; /* store the good pixel value */
1443 ngoodpix++;
1444
1445 if (do_range) {
1446 if (v8 < xminval) xminval = v8;
1447 if (v8 > xmaxval) xmaxval = v8;
1448 }
1449 /* now populate the differences arrays */
1450 /* for the remaining pixels in the row */
1451 nvals = 0;
1452 nvals2 = 0;
1453 for (ii++; ii < nx; ii++) {
1454
1455 /* find the next valid pixel in row */
1456 if (nullcheck)
1457 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1458
1459 if (ii == nx) break; /* hit end of row */
1460 v9 = rowpix[ii]; /* store the good pixel value */
1461
1462 if (do_range) {
1463 if (v9 < xminval) xminval = v9;
1464 if (v9 > xmaxval) xmaxval = v9;
1465 }
1466
1467 /* construct array of absolute differences */
1468
1469 if (!(v5 == v6 && v6 == v7) ) {
1470 tdiff = (LONGLONG) v5 - (LONGLONG) v7;
1471 if (tdiff < 0)
1472 differences2[nvals2] = -1 * tdiff;
1473 else
1474 differences2[nvals2] = tdiff;
1475
1476 nvals2++;
1477 }
1478
1479 if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
1480 tdiff = (2 * (LONGLONG) v5) - (LONGLONG) v3 - (LONGLONG) v7;
1481 if (tdiff < 0)
1482 differences3[nvals] = -1 * tdiff;
1483 else
1484 differences3[nvals] = tdiff;
1485
1486 tdiff = (6 * (LONGLONG) v5) - (4 * (LONGLONG) v3) - (4 * (LONGLONG) v7) + (LONGLONG) v1 + (LONGLONG) v9;
1487 if (tdiff < 0)
1488 differences5[nvals] = -1 * tdiff;
1489 else
1490 differences5[nvals] = tdiff;
1491
1492 nvals++;
1493 } else {
1494 /* ignore constant background regions */
1495 ngoodpix++;
1496 }
1497
1498 /* shift over 1 pixel */
1499 v1 = v2;
1500 v2 = v3;
1501 v3 = v4;
1502 v4 = v5;
1503 v5 = v6;
1504 v6 = v7;
1505 v7 = v8;
1506 v8 = v9;
1507 } /* end of loop over pixels in the row */
1508
1509 /* compute the median diffs */
1510 /* Note that there are 8 more pixel values than there are diffs values. */
1511 ngoodpix += nvals;
1512
1513 if (nvals == 0) {
1514 continue; /* cannot compute medians on this row */
1515 } else if (nvals == 1) {
1516 if (nvals2 == 1) {
1517 diffs2[nrows2] = (double) differences2[0];
1518 nrows2++;
1519 }
1520
1521 diffs3[nrows] = (double) differences3[0];
1522 diffs5[nrows] = (double) differences5[0];
1523 } else {
1524 /* quick_select returns the median MUCH faster than using qsort */
1525 if (nvals2 > 1) {
1526 diffs2[nrows2] = (double) quick_select_longlong(differences2, nvals);
1527 nrows2++;
1528 }
1529
1530 diffs3[nrows] = (double) quick_select_longlong(differences3, nvals);
1531 diffs5[nrows] = (double) quick_select_longlong(differences5, nvals);
1532 }
1533
1534 nrows++;
1535 } /* end of loop over rows */
1536
1537 /* compute median of the values for each row */
1538 if (nrows == 0) {
1539 xnoise3 = 0;
1540 xnoise5 = 0;
1541 } else if (nrows == 1) {
1542 xnoise3 = diffs3[0];
1543 xnoise5 = diffs5[0];
1544 } else {
1545 qsort(diffs3, nrows, sizeof(double), FnCompare_double);
1546 qsort(diffs5, nrows, sizeof(double), FnCompare_double);
1547 xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
1548 xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
1549 }
1550
1551 if (nrows2 == 0) {
1552 xnoise2 = 0;
1553 } else if (nrows2 == 1) {
1554 xnoise2 = diffs2[0];
1555 } else {
1556 qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
1557 xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
1558 }
1559
1560 if (ngood) *ngood = ngoodpix;
1561 if (minval) *minval = xminval;
1562 if (maxval) *maxval = xmaxval;
1563 if (noise2) *noise2 = 1.0483579 * xnoise2;
1564 if (noise3) *noise3 = 0.6052697 * xnoise3;
1565 if (noise5) *noise5 = 0.1772048 * xnoise5;
1566
1567 free(diffs5);
1568 free(diffs3);
1569 free(diffs2);
1570 free(differences5);
1571 free(differences3);
1572 free(differences2);
1573
1574 return(*status);
1575 }
1576 /*--------------------------------------------------------------------------*/
FnNoise5_float(float * array,long nx,long ny,int nullcheck,float nullvalue,long * ngood,float * minval,float * maxval,double * noise2,double * noise3,double * noise5,int * status)1577 static int FnNoise5_float
1578 (float *array, /* 2 dimensional array of image pixels */
1579 long nx, /* number of pixels in each row of the image */
1580 long ny, /* number of rows in the image */
1581 int nullcheck, /* check for null values, if true */
1582 float nullvalue, /* value of null pixels, if nullcheck is true */
1583 /* returned parameters */
1584 long *ngood, /* number of good, non-null pixels? */
1585 float *minval, /* minimum non-null value */
1586 float *maxval, /* maximum non-null value */
1587 double *noise2, /* returned 2nd order MAD of all non-null pixels */
1588 double *noise3, /* returned 3rd order MAD of all non-null pixels */
1589 double *noise5, /* returned 5th order MAD of all non-null pixels */
1590 int *status) /* error status */
1591
1592 /*
1593 Estimate the median and background noise in the input image using 2nd, 3rd and 5th
1594 order Median Absolute Differences.
1595
1596 The noise in the background of the image is calculated using the MAD algorithms
1597 developed for deriving the signal to noise ratio in spectra
1598 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
1599
1600 3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
1601
1602 The returned estimates are the median of the values that are computed for each
1603 row of the image.
1604 */
1605 {
1606 long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
1607 float *differences2, *differences3, *differences5;
1608 float *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
1609 float xminval = FLT_MAX, xmaxval = -FLT_MAX;
1610 int do_range = 0;
1611 double *diffs2, *diffs3, *diffs5;
1612 double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
1613
1614 if (nx < 9) {
1615 /* treat entire array as an image with a single row */
1616 nx = nx * ny;
1617 ny = 1;
1618 }
1619
1620 /* rows must have at least 9 pixels */
1621 if (nx < 9) {
1622
1623 for (ii = 0; ii < nx; ii++) {
1624 if (nullcheck && array[ii] == nullvalue)
1625 continue;
1626 else {
1627 if (array[ii] < xminval) xminval = array[ii];
1628 if (array[ii] > xmaxval) xmaxval = array[ii];
1629 ngoodpix++;
1630 }
1631 }
1632 if (minval) *minval = xminval;
1633 if (maxval) *maxval = xmaxval;
1634 if (ngood) *ngood = ngoodpix;
1635 if (noise2) *noise2 = 0.;
1636 if (noise3) *noise3 = 0.;
1637 if (noise5) *noise5 = 0.;
1638 return(*status);
1639 }
1640
1641 /* do we need to compute the min and max value? */
1642 if (minval || maxval) do_range = 1;
1643
1644 /* allocate arrays used to compute the median and noise estimates */
1645 differences2 = calloc(nx, sizeof(float));
1646 if (!differences2) {
1647 *status = MEMORY_ALLOCATION;
1648 return(*status);
1649 }
1650 differences3 = calloc(nx, sizeof(float));
1651 if (!differences3) {
1652 free(differences2);
1653 *status = MEMORY_ALLOCATION;
1654 return(*status);
1655 }
1656 differences5 = calloc(nx, sizeof(float));
1657 if (!differences5) {
1658 free(differences2);
1659 free(differences3);
1660 *status = MEMORY_ALLOCATION;
1661 return(*status);
1662 }
1663
1664 diffs2 = calloc(ny, sizeof(double));
1665 if (!diffs2) {
1666 free(differences2);
1667 free(differences3);
1668 free(differences5);
1669 *status = MEMORY_ALLOCATION;
1670 return(*status);
1671 }
1672
1673 diffs3 = calloc(ny, sizeof(double));
1674 if (!diffs3) {
1675 free(differences2);
1676 free(differences3);
1677 free(differences5);
1678 free(diffs2);
1679 *status = MEMORY_ALLOCATION;
1680 return(*status);
1681 }
1682
1683 diffs5 = calloc(ny, sizeof(double));
1684 if (!diffs5) {
1685 free(differences2);
1686 free(differences3);
1687 free(differences5);
1688 free(diffs2);
1689 free(diffs3);
1690 *status = MEMORY_ALLOCATION;
1691 return(*status);
1692 }
1693
1694 /* loop over each row of the image */
1695 for (jj=0; jj < ny; jj++) {
1696
1697 rowpix = array + (jj * nx); /* point to first pixel in the row */
1698
1699 /***** find the first valid pixel in row */
1700 ii = 0;
1701 if (nullcheck)
1702 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1703
1704 if (ii == nx) continue; /* hit end of row */
1705 v1 = rowpix[ii]; /* store the good pixel value */
1706 ngoodpix++;
1707
1708 if (do_range) {
1709 if (v1 < xminval) xminval = v1;
1710 if (v1 > xmaxval) xmaxval = v1;
1711 }
1712
1713 /***** find the 2nd valid pixel in row (which we will skip over) */
1714 ii++;
1715 if (nullcheck)
1716 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1717
1718 if (ii == nx) continue; /* hit end of row */
1719 v2 = rowpix[ii]; /* store the good pixel value */
1720 ngoodpix++;
1721
1722 if (do_range) {
1723 if (v2 < xminval) xminval = v2;
1724 if (v2 > xmaxval) xmaxval = v2;
1725 }
1726
1727 /***** find the 3rd valid pixel in row */
1728 ii++;
1729 if (nullcheck)
1730 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1731
1732 if (ii == nx) continue; /* hit end of row */
1733 v3 = rowpix[ii]; /* store the good pixel value */
1734 ngoodpix++;
1735
1736 if (do_range) {
1737 if (v3 < xminval) xminval = v3;
1738 if (v3 > xmaxval) xmaxval = v3;
1739 }
1740
1741 /* find the 4nd valid pixel in row (to be skipped) */
1742 ii++;
1743 if (nullcheck)
1744 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1745
1746 if (ii == nx) continue; /* hit end of row */
1747 v4 = rowpix[ii]; /* store the good pixel value */
1748 ngoodpix++;
1749
1750 if (do_range) {
1751 if (v4 < xminval) xminval = v4;
1752 if (v4 > xmaxval) xmaxval = v4;
1753 }
1754
1755 /* find the 5th valid pixel in row (to be skipped) */
1756 ii++;
1757 if (nullcheck)
1758 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1759
1760 if (ii == nx) continue; /* hit end of row */
1761 v5 = rowpix[ii]; /* store the good pixel value */
1762 ngoodpix++;
1763
1764 if (do_range) {
1765 if (v5 < xminval) xminval = v5;
1766 if (v5 > xmaxval) xmaxval = v5;
1767 }
1768
1769 /* find the 6th valid pixel in row (to be skipped) */
1770 ii++;
1771 if (nullcheck)
1772 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1773
1774 if (ii == nx) continue; /* hit end of row */
1775 v6 = rowpix[ii]; /* store the good pixel value */
1776 ngoodpix++;
1777
1778 if (do_range) {
1779 if (v6 < xminval) xminval = v6;
1780 if (v6 > xmaxval) xmaxval = v6;
1781 }
1782
1783 /* find the 7th valid pixel in row (to be skipped) */
1784 ii++;
1785 if (nullcheck)
1786 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1787
1788 if (ii == nx) continue; /* hit end of row */
1789 v7 = rowpix[ii]; /* store the good pixel value */
1790 ngoodpix++;
1791
1792 if (do_range) {
1793 if (v7 < xminval) xminval = v7;
1794 if (v7 > xmaxval) xmaxval = v7;
1795 }
1796
1797 /* find the 8th valid pixel in row (to be skipped) */
1798 ii++;
1799 if (nullcheck)
1800 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1801
1802 if (ii == nx) continue; /* hit end of row */
1803 v8 = rowpix[ii]; /* store the good pixel value */
1804 ngoodpix++;
1805
1806 if (do_range) {
1807 if (v8 < xminval) xminval = v8;
1808 if (v8 > xmaxval) xmaxval = v8;
1809 }
1810 /* now populate the differences arrays */
1811 /* for the remaining pixels in the row */
1812 nvals = 0;
1813 nvals2 = 0;
1814 for (ii++; ii < nx; ii++) {
1815
1816 /* find the next valid pixel in row */
1817 if (nullcheck)
1818 while (ii < nx && rowpix[ii] == nullvalue) ii++;
1819
1820 if (ii == nx) break; /* hit end of row */
1821 v9 = rowpix[ii]; /* store the good pixel value */
1822
1823 if (do_range) {
1824 if (v9 < xminval) xminval = v9;
1825 if (v9 > xmaxval) xmaxval = v9;
1826 }
1827
1828 /* construct array of absolute differences */
1829
1830 if (!(v5 == v6 && v6 == v7) ) {
1831 differences2[nvals2] = (float) fabs(v5 - v7);
1832 nvals2++;
1833 }
1834
1835 if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
1836 differences3[nvals] = (float) fabs((2 * v5) - v3 - v7);
1837 differences5[nvals] = (float) fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9);
1838 nvals++;
1839 } else {
1840 /* ignore constant background regions */
1841 ngoodpix++;
1842 }
1843
1844 /* shift over 1 pixel */
1845 v1 = v2;
1846 v2 = v3;
1847 v3 = v4;
1848 v4 = v5;
1849 v5 = v6;
1850 v6 = v7;
1851 v7 = v8;
1852 v8 = v9;
1853 } /* end of loop over pixels in the row */
1854
1855 /* compute the median diffs */
1856 /* Note that there are 8 more pixel values than there are diffs values. */
1857 ngoodpix += nvals;
1858
1859 if (nvals == 0) {
1860 continue; /* cannot compute medians on this row */
1861 } else if (nvals == 1) {
1862 if (nvals2 == 1) {
1863 diffs2[nrows2] = differences2[0];
1864 nrows2++;
1865 }
1866
1867 diffs3[nrows] = differences3[0];
1868 diffs5[nrows] = differences5[0];
1869 } else {
1870 /* quick_select returns the median MUCH faster than using qsort */
1871 if (nvals2 > 1) {
1872 diffs2[nrows2] = quick_select_float(differences2, nvals);
1873 nrows2++;
1874 }
1875
1876 diffs3[nrows] = quick_select_float(differences3, nvals);
1877 diffs5[nrows] = quick_select_float(differences5, nvals);
1878 }
1879
1880 nrows++;
1881 } /* end of loop over rows */
1882
1883 /* compute median of the values for each row */
1884 if (nrows == 0) {
1885 xnoise3 = 0;
1886 xnoise5 = 0;
1887 } else if (nrows == 1) {
1888 xnoise3 = diffs3[0];
1889 xnoise5 = diffs5[0];
1890 } else {
1891 qsort(diffs3, nrows, sizeof(double), FnCompare_double);
1892 qsort(diffs5, nrows, sizeof(double), FnCompare_double);
1893 xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
1894 xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
1895 }
1896
1897 if (nrows2 == 0) {
1898 xnoise2 = 0;
1899 } else if (nrows2 == 1) {
1900 xnoise2 = diffs2[0];
1901 } else {
1902 qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
1903 xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
1904 }
1905
1906 if (ngood) *ngood = ngoodpix;
1907 if (minval) *minval = xminval;
1908 if (maxval) *maxval = xmaxval;
1909 if (noise2) *noise2 = 1.0483579 * xnoise2;
1910 if (noise3) *noise3 = 0.6052697 * xnoise3;
1911 if (noise5) *noise5 = 0.1772048 * xnoise5;
1912
1913 free(diffs5);
1914 free(diffs3);
1915 free(diffs2);
1916 free(differences5);
1917 free(differences3);
1918 free(differences2);
1919
1920 return(*status);
1921 }
1922 /*--------------------------------------------------------------------------*/
FnNoise5_double(double * array,long nx,long ny,int nullcheck,double nullvalue,long * ngood,double * minval,double * maxval,double * noise2,double * noise3,double * noise5,int * status)1923 static int FnNoise5_double
1924 (double *array, /* 2 dimensional array of image pixels */
1925 long nx, /* number of pixels in each row of the image */
1926 long ny, /* number of rows in the image */
1927 int nullcheck, /* check for null values, if true */
1928 double nullvalue, /* value of null pixels, if nullcheck is true */
1929 /* returned parameters */
1930 long *ngood, /* number of good, non-null pixels? */
1931 double *minval, /* minimum non-null value */
1932 double *maxval, /* maximum non-null value */
1933 double *noise2, /* returned 2nd order MAD of all non-null pixels */
1934 double *noise3, /* returned 3rd order MAD of all non-null pixels */
1935 double *noise5, /* returned 5th order MAD of all non-null pixels */
1936 int *status) /* error status */
1937
1938 /*
1939 Estimate the median and background noise in the input image using 2nd, 3rd and 5th
1940 order Median Absolute Differences.
1941
1942 The noise in the background of the image is calculated using the MAD algorithms
1943 developed for deriving the signal to noise ratio in spectra
1944 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
1945
1946 3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
1947
1948 The returned estimates are the median of the values that are computed for each
1949 row of the image.
1950 */
1951 {
1952 long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0;
1953 double *differences2, *differences3, *differences5;
1954 double *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9;
1955 double xminval = DBL_MAX, xmaxval = -DBL_MAX;
1956 int do_range = 0;
1957 double *diffs2, *diffs3, *diffs5;
1958 double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0;
1959
1960 if (nx < 9) {
1961 /* treat entire array as an image with a single row */
1962 nx = nx * ny;
1963 ny = 1;
1964 }
1965
1966 /* rows must have at least 9 pixels */
1967 if (nx < 9) {
1968
1969 for (ii = 0; ii < nx; ii++) {
1970 if (nullcheck && array[ii] == nullvalue)
1971 continue;
1972 else {
1973 if (array[ii] < xminval) xminval = array[ii];
1974 if (array[ii] > xmaxval) xmaxval = array[ii];
1975 ngoodpix++;
1976 }
1977 }
1978 if (minval) *minval = xminval;
1979 if (maxval) *maxval = xmaxval;
1980 if (ngood) *ngood = ngoodpix;
1981 if (noise2) *noise2 = 0.;
1982 if (noise3) *noise3 = 0.;
1983 if (noise5) *noise5 = 0.;
1984 return(*status);
1985 }
1986
1987 /* do we need to compute the min and max value? */
1988 if (minval || maxval) do_range = 1;
1989
1990 /* allocate arrays used to compute the median and noise estimates */
1991 differences2 = calloc(nx, sizeof(double));
1992 if (!differences2) {
1993 *status = MEMORY_ALLOCATION;
1994 return(*status);
1995 }
1996 differences3 = calloc(nx, sizeof(double));
1997 if (!differences3) {
1998 free(differences2);
1999 *status = MEMORY_ALLOCATION;
2000 return(*status);
2001 }
2002 differences5 = calloc(nx, sizeof(double));
2003 if (!differences5) {
2004 free(differences2);
2005 free(differences3);
2006 *status = MEMORY_ALLOCATION;
2007 return(*status);
2008 }
2009
2010 diffs2 = calloc(ny, sizeof(double));
2011 if (!diffs2) {
2012 free(differences2);
2013 free(differences3);
2014 free(differences5);
2015 *status = MEMORY_ALLOCATION;
2016 return(*status);
2017 }
2018
2019 diffs3 = calloc(ny, sizeof(double));
2020 if (!diffs3) {
2021 free(differences2);
2022 free(differences3);
2023 free(differences5);
2024 free(diffs2);
2025 *status = MEMORY_ALLOCATION;
2026 return(*status);
2027 }
2028
2029 diffs5 = calloc(ny, sizeof(double));
2030 if (!diffs5) {
2031 free(differences2);
2032 free(differences3);
2033 free(differences5);
2034 free(diffs2);
2035 free(diffs3);
2036 *status = MEMORY_ALLOCATION;
2037 return(*status);
2038 }
2039
2040 /* loop over each row of the image */
2041 for (jj=0; jj < ny; jj++) {
2042
2043 rowpix = array + (jj * nx); /* point to first pixel in the row */
2044
2045 /***** find the first valid pixel in row */
2046 ii = 0;
2047 if (nullcheck)
2048 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2049
2050 if (ii == nx) continue; /* hit end of row */
2051 v1 = rowpix[ii]; /* store the good pixel value */
2052 ngoodpix++;
2053
2054 if (do_range) {
2055 if (v1 < xminval) xminval = v1;
2056 if (v1 > xmaxval) xmaxval = v1;
2057 }
2058
2059 /***** find the 2nd valid pixel in row (which we will skip over) */
2060 ii++;
2061 if (nullcheck)
2062 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2063
2064 if (ii == nx) continue; /* hit end of row */
2065 v2 = rowpix[ii]; /* store the good pixel value */
2066 ngoodpix++;
2067
2068 if (do_range) {
2069 if (v2 < xminval) xminval = v2;
2070 if (v2 > xmaxval) xmaxval = v2;
2071 }
2072
2073 /***** find the 3rd valid pixel in row */
2074 ii++;
2075 if (nullcheck)
2076 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2077
2078 if (ii == nx) continue; /* hit end of row */
2079 v3 = rowpix[ii]; /* store the good pixel value */
2080 ngoodpix++;
2081
2082 if (do_range) {
2083 if (v3 < xminval) xminval = v3;
2084 if (v3 > xmaxval) xmaxval = v3;
2085 }
2086
2087 /* find the 4nd valid pixel in row (to be skipped) */
2088 ii++;
2089 if (nullcheck)
2090 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2091
2092 if (ii == nx) continue; /* hit end of row */
2093 v4 = rowpix[ii]; /* store the good pixel value */
2094 ngoodpix++;
2095
2096 if (do_range) {
2097 if (v4 < xminval) xminval = v4;
2098 if (v4 > xmaxval) xmaxval = v4;
2099 }
2100
2101 /* find the 5th valid pixel in row (to be skipped) */
2102 ii++;
2103 if (nullcheck)
2104 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2105
2106 if (ii == nx) continue; /* hit end of row */
2107 v5 = rowpix[ii]; /* store the good pixel value */
2108 ngoodpix++;
2109
2110 if (do_range) {
2111 if (v5 < xminval) xminval = v5;
2112 if (v5 > xmaxval) xmaxval = v5;
2113 }
2114
2115 /* find the 6th valid pixel in row (to be skipped) */
2116 ii++;
2117 if (nullcheck)
2118 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2119
2120 if (ii == nx) continue; /* hit end of row */
2121 v6 = rowpix[ii]; /* store the good pixel value */
2122 ngoodpix++;
2123
2124 if (do_range) {
2125 if (v6 < xminval) xminval = v6;
2126 if (v6 > xmaxval) xmaxval = v6;
2127 }
2128
2129 /* find the 7th valid pixel in row (to be skipped) */
2130 ii++;
2131 if (nullcheck)
2132 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2133
2134 if (ii == nx) continue; /* hit end of row */
2135 v7 = rowpix[ii]; /* store the good pixel value */
2136 ngoodpix++;
2137
2138 if (do_range) {
2139 if (v7 < xminval) xminval = v7;
2140 if (v7 > xmaxval) xmaxval = v7;
2141 }
2142
2143 /* find the 8th valid pixel in row (to be skipped) */
2144 ii++;
2145 if (nullcheck)
2146 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2147
2148 if (ii == nx) continue; /* hit end of row */
2149 v8 = rowpix[ii]; /* store the good pixel value */
2150 ngoodpix++;
2151
2152 if (do_range) {
2153 if (v8 < xminval) xminval = v8;
2154 if (v8 > xmaxval) xmaxval = v8;
2155 }
2156 /* now populate the differences arrays */
2157 /* for the remaining pixels in the row */
2158 nvals = 0;
2159 nvals2 = 0;
2160 for (ii++; ii < nx; ii++) {
2161
2162 /* find the next valid pixel in row */
2163 if (nullcheck)
2164 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2165
2166 if (ii == nx) break; /* hit end of row */
2167 v9 = rowpix[ii]; /* store the good pixel value */
2168
2169 if (do_range) {
2170 if (v9 < xminval) xminval = v9;
2171 if (v9 > xmaxval) xmaxval = v9;
2172 }
2173
2174 /* construct array of absolute differences */
2175
2176 if (!(v5 == v6 && v6 == v7) ) {
2177 differences2[nvals2] = fabs(v5 - v7);
2178 nvals2++;
2179 }
2180
2181 if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) {
2182 differences3[nvals] = fabs((2 * v5) - v3 - v7);
2183 differences5[nvals] = fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9);
2184 nvals++;
2185 } else {
2186 /* ignore constant background regions */
2187 ngoodpix++;
2188 }
2189
2190 /* shift over 1 pixel */
2191 v1 = v2;
2192 v2 = v3;
2193 v3 = v4;
2194 v4 = v5;
2195 v5 = v6;
2196 v6 = v7;
2197 v7 = v8;
2198 v8 = v9;
2199 } /* end of loop over pixels in the row */
2200
2201 /* compute the median diffs */
2202 /* Note that there are 8 more pixel values than there are diffs values. */
2203 ngoodpix += nvals;
2204
2205 if (nvals == 0) {
2206 continue; /* cannot compute medians on this row */
2207 } else if (nvals == 1) {
2208 if (nvals2 == 1) {
2209 diffs2[nrows2] = differences2[0];
2210 nrows2++;
2211 }
2212
2213 diffs3[nrows] = differences3[0];
2214 diffs5[nrows] = differences5[0];
2215 } else {
2216 /* quick_select returns the median MUCH faster than using qsort */
2217 if (nvals2 > 1) {
2218 diffs2[nrows2] = quick_select_double(differences2, nvals);
2219 nrows2++;
2220 }
2221
2222 diffs3[nrows] = quick_select_double(differences3, nvals);
2223 diffs5[nrows] = quick_select_double(differences5, nvals);
2224 }
2225
2226 nrows++;
2227 } /* end of loop over rows */
2228
2229 /* compute median of the values for each row */
2230 if (nrows == 0) {
2231 xnoise3 = 0;
2232 xnoise5 = 0;
2233 } else if (nrows == 1) {
2234 xnoise3 = diffs3[0];
2235 xnoise5 = diffs5[0];
2236 } else {
2237 qsort(diffs3, nrows, sizeof(double), FnCompare_double);
2238 qsort(diffs5, nrows, sizeof(double), FnCompare_double);
2239 xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.;
2240 xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.;
2241 }
2242
2243 if (nrows2 == 0) {
2244 xnoise2 = 0;
2245 } else if (nrows2 == 1) {
2246 xnoise2 = diffs2[0];
2247 } else {
2248 qsort(diffs2, nrows2, sizeof(double), FnCompare_double);
2249 xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.;
2250 }
2251
2252 if (ngood) *ngood = ngoodpix;
2253 if (minval) *minval = xminval;
2254 if (maxval) *maxval = xmaxval;
2255 if (noise2) *noise2 = 1.0483579 * xnoise2;
2256 if (noise3) *noise3 = 0.6052697 * xnoise3;
2257 if (noise5) *noise5 = 0.1772048 * xnoise5;
2258
2259 free(diffs5);
2260 free(diffs3);
2261 free(diffs2);
2262 free(differences5);
2263 free(differences3);
2264 free(differences2);
2265
2266 return(*status);
2267 }
2268 /*--------------------------------------------------------------------------*/
FnNoise3_short(short * array,long nx,long ny,int nullcheck,short nullvalue,long * ngood,short * minval,short * maxval,double * noise,int * status)2269 static int FnNoise3_short
2270 (short *array, /* 2 dimensional array of image pixels */
2271 long nx, /* number of pixels in each row of the image */
2272 long ny, /* number of rows in the image */
2273 int nullcheck, /* check for null values, if true */
2274 short nullvalue, /* value of null pixels, if nullcheck is true */
2275 /* returned parameters */
2276 long *ngood, /* number of good, non-null pixels? */
2277 short *minval, /* minimum non-null value */
2278 short *maxval, /* maximum non-null value */
2279 double *noise, /* returned R.M.S. value of all non-null pixels */
2280 int *status) /* error status */
2281
2282 /*
2283 Estimate the median and background noise in the input image using 3rd order differences.
2284
2285 The noise in the background of the image is calculated using the 3rd order algorithm
2286 developed for deriving the signal to noise ratio in spectra
2287 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
2288
2289 noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
2290
2291 The returned estimates are the median of the values that are computed for each
2292 row of the image.
2293 */
2294 {
2295 long ii, jj, nrows = 0, nvals, ngoodpix = 0;
2296 short *differences, *rowpix, v1, v2, v3, v4, v5;
2297 short xminval = SHRT_MAX, xmaxval = SHRT_MIN, do_range = 0;
2298 double *diffs, xnoise = 0, sigma;
2299
2300 if (nx < 5) {
2301 /* treat entire array as an image with a single row */
2302 nx = nx * ny;
2303 ny = 1;
2304 }
2305
2306 /* rows must have at least 5 pixels */
2307 if (nx < 5) {
2308
2309 for (ii = 0; ii < nx; ii++) {
2310 if (nullcheck && array[ii] == nullvalue)
2311 continue;
2312 else {
2313 if (array[ii] < xminval) xminval = array[ii];
2314 if (array[ii] > xmaxval) xmaxval = array[ii];
2315 ngoodpix++;
2316 }
2317 }
2318 if (minval) *minval = xminval;
2319 if (maxval) *maxval = xmaxval;
2320 if (ngood) *ngood = ngoodpix;
2321 if (noise) *noise = 0.;
2322 return(*status);
2323 }
2324
2325 /* do we need to compute the min and max value? */
2326 if (minval || maxval) do_range = 1;
2327
2328 /* allocate arrays used to compute the median and noise estimates */
2329 differences = calloc(nx, sizeof(short));
2330 if (!differences) {
2331 *status = MEMORY_ALLOCATION;
2332 return(*status);
2333 }
2334
2335 diffs = calloc(ny, sizeof(double));
2336 if (!diffs) {
2337 free(differences);
2338 *status = MEMORY_ALLOCATION;
2339 return(*status);
2340 }
2341
2342 /* loop over each row of the image */
2343 for (jj=0; jj < ny; jj++) {
2344
2345 rowpix = array + (jj * nx); /* point to first pixel in the row */
2346
2347 /***** find the first valid pixel in row */
2348 ii = 0;
2349 if (nullcheck)
2350 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2351
2352 if (ii == nx) continue; /* hit end of row */
2353 v1 = rowpix[ii]; /* store the good pixel value */
2354
2355 if (do_range) {
2356 if (v1 < xminval) xminval = v1;
2357 if (v1 > xmaxval) xmaxval = v1;
2358 }
2359
2360 /***** find the 2nd valid pixel in row (which we will skip over) */
2361 ii++;
2362 if (nullcheck)
2363 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2364
2365 if (ii == nx) continue; /* hit end of row */
2366 v2 = rowpix[ii]; /* store the good pixel value */
2367
2368 if (do_range) {
2369 if (v2 < xminval) xminval = v2;
2370 if (v2 > xmaxval) xmaxval = v2;
2371 }
2372
2373 /***** find the 3rd valid pixel in row */
2374 ii++;
2375 if (nullcheck)
2376 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2377
2378 if (ii == nx) continue; /* hit end of row */
2379 v3 = rowpix[ii]; /* store the good pixel value */
2380
2381 if (do_range) {
2382 if (v3 < xminval) xminval = v3;
2383 if (v3 > xmaxval) xmaxval = v3;
2384 }
2385
2386 /* find the 4nd valid pixel in row (to be skipped) */
2387 ii++;
2388 if (nullcheck)
2389 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2390
2391 if (ii == nx) continue; /* hit end of row */
2392 v4 = rowpix[ii]; /* store the good pixel value */
2393
2394 if (do_range) {
2395 if (v4 < xminval) xminval = v4;
2396 if (v4 > xmaxval) xmaxval = v4;
2397 }
2398
2399 /* now populate the differences arrays */
2400 /* for the remaining pixels in the row */
2401 nvals = 0;
2402 for (ii++; ii < nx; ii++) {
2403
2404 /* find the next valid pixel in row */
2405 if (nullcheck)
2406 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2407
2408 if (ii == nx) break; /* hit end of row */
2409 v5 = rowpix[ii]; /* store the good pixel value */
2410
2411 if (do_range) {
2412 if (v5 < xminval) xminval = v5;
2413 if (v5 > xmaxval) xmaxval = v5;
2414 }
2415
2416 /* construct array of 3rd order absolute differences */
2417 if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
2418 differences[nvals] = abs((2 * v3) - v1 - v5);
2419 nvals++;
2420 } else {
2421 /* ignore constant background regions */
2422 ngoodpix++;
2423 }
2424
2425
2426 /* shift over 1 pixel */
2427 v1 = v2;
2428 v2 = v3;
2429 v3 = v4;
2430 v4 = v5;
2431 } /* end of loop over pixels in the row */
2432
2433 /* compute the 3rd order diffs */
2434 /* Note that there are 4 more pixel values than there are diffs values. */
2435 ngoodpix += (nvals + 4);
2436
2437 if (nvals == 0) {
2438 continue; /* cannot compute medians on this row */
2439 } else if (nvals == 1) {
2440 diffs[nrows] = differences[0];
2441 } else {
2442 /* quick_select returns the median MUCH faster than using qsort */
2443 diffs[nrows] = quick_select_short(differences, nvals);
2444 }
2445
2446 nrows++;
2447 } /* end of loop over rows */
2448
2449 /* compute median of the values for each row */
2450 if (nrows == 0) {
2451 xnoise = 0;
2452 } else if (nrows == 1) {
2453 xnoise = diffs[0];
2454 } else {
2455
2456
2457 qsort(diffs, nrows, sizeof(double), FnCompare_double);
2458 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
2459
2460 FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status);
2461
2462 /* do a 4.5 sigma rejection of outliers */
2463 jj = 0;
2464 sigma = 4.5 * sigma;
2465 for (ii = 0; ii < nrows; ii++) {
2466 if ( fabs(diffs[ii] - xnoise) <= sigma) {
2467 if (jj != ii)
2468 diffs[jj] = diffs[ii];
2469 jj++;
2470 }
2471 }
2472 if (ii != jj)
2473 FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status);
2474 }
2475
2476 if (ngood) *ngood = ngoodpix;
2477 if (minval) *minval = xminval;
2478 if (maxval) *maxval = xmaxval;
2479 if (noise) *noise = 0.6052697 * xnoise;
2480
2481 free(diffs);
2482 free(differences);
2483
2484 return(*status);
2485 }
2486 /*--------------------------------------------------------------------------*/
FnNoise3_int(int * array,long nx,long ny,int nullcheck,int nullvalue,long * ngood,int * minval,int * maxval,double * noise,int * status)2487 static int FnNoise3_int
2488 (int *array, /* 2 dimensional array of image pixels */
2489 long nx, /* number of pixels in each row of the image */
2490 long ny, /* number of rows in the image */
2491 int nullcheck, /* check for null values, if true */
2492 int nullvalue, /* value of null pixels, if nullcheck is true */
2493 /* returned parameters */
2494 long *ngood, /* number of good, non-null pixels? */
2495 int *minval, /* minimum non-null value */
2496 int *maxval, /* maximum non-null value */
2497 double *noise, /* returned R.M.S. value of all non-null pixels */
2498 int *status) /* error status */
2499
2500 /*
2501 Estimate the background noise in the input image using 3rd order differences.
2502
2503 The noise in the background of the image is calculated using the 3rd order algorithm
2504 developed for deriving the signal to noise ratio in spectra
2505 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
2506
2507 noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
2508
2509 The returned estimates are the median of the values that are computed for each
2510 row of the image.
2511 */
2512 {
2513 long ii, jj, nrows = 0, nvals, ngoodpix = 0;
2514 int *differences, *rowpix, v1, v2, v3, v4, v5;
2515 int xminval = INT_MAX, xmaxval = INT_MIN, do_range = 0;
2516 double *diffs, xnoise = 0, sigma;
2517
2518 if (nx < 5) {
2519 /* treat entire array as an image with a single row */
2520 nx = nx * ny;
2521 ny = 1;
2522 }
2523
2524 /* rows must have at least 5 pixels */
2525 if (nx < 5) {
2526
2527 for (ii = 0; ii < nx; ii++) {
2528 if (nullcheck && array[ii] == nullvalue)
2529 continue;
2530 else {
2531 if (array[ii] < xminval) xminval = array[ii];
2532 if (array[ii] > xmaxval) xmaxval = array[ii];
2533 ngoodpix++;
2534 }
2535 }
2536 if (minval) *minval = xminval;
2537 if (maxval) *maxval = xmaxval;
2538 if (ngood) *ngood = ngoodpix;
2539 if (noise) *noise = 0.;
2540 return(*status);
2541 }
2542
2543 /* do we need to compute the min and max value? */
2544 if (minval || maxval) do_range = 1;
2545
2546 /* allocate arrays used to compute the median and noise estimates */
2547 differences = calloc(nx, sizeof(int));
2548 if (!differences) {
2549 *status = MEMORY_ALLOCATION;
2550 return(*status);
2551 }
2552
2553 diffs = calloc(ny, sizeof(double));
2554 if (!diffs) {
2555 free(differences);
2556 *status = MEMORY_ALLOCATION;
2557 return(*status);
2558 }
2559
2560 /* loop over each row of the image */
2561 for (jj=0; jj < ny; jj++) {
2562
2563 rowpix = array + (jj * nx); /* point to first pixel in the row */
2564
2565 /***** find the first valid pixel in row */
2566 ii = 0;
2567 if (nullcheck)
2568 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2569
2570 if (ii == nx) continue; /* hit end of row */
2571 v1 = rowpix[ii]; /* store the good pixel value */
2572
2573 if (do_range) {
2574 if (v1 < xminval) xminval = v1;
2575 if (v1 > xmaxval) xmaxval = v1;
2576 }
2577
2578 /***** find the 2nd valid pixel in row (which we will skip over) */
2579 ii++;
2580 if (nullcheck)
2581 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2582
2583 if (ii == nx) continue; /* hit end of row */
2584 v2 = rowpix[ii]; /* store the good pixel value */
2585
2586 if (do_range) {
2587 if (v2 < xminval) xminval = v2;
2588 if (v2 > xmaxval) xmaxval = v2;
2589 }
2590
2591 /***** find the 3rd valid pixel in row */
2592 ii++;
2593 if (nullcheck)
2594 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2595
2596 if (ii == nx) continue; /* hit end of row */
2597 v3 = rowpix[ii]; /* store the good pixel value */
2598
2599 if (do_range) {
2600 if (v3 < xminval) xminval = v3;
2601 if (v3 > xmaxval) xmaxval = v3;
2602 }
2603
2604 /* find the 4nd valid pixel in row (to be skipped) */
2605 ii++;
2606 if (nullcheck)
2607 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2608
2609 if (ii == nx) continue; /* hit end of row */
2610 v4 = rowpix[ii]; /* store the good pixel value */
2611
2612 if (do_range) {
2613 if (v4 < xminval) xminval = v4;
2614 if (v4 > xmaxval) xmaxval = v4;
2615 }
2616
2617 /* now populate the differences arrays */
2618 /* for the remaining pixels in the row */
2619 nvals = 0;
2620 for (ii++; ii < nx; ii++) {
2621
2622 /* find the next valid pixel in row */
2623 if (nullcheck)
2624 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2625
2626 if (ii == nx) break; /* hit end of row */
2627 v5 = rowpix[ii]; /* store the good pixel value */
2628
2629 if (do_range) {
2630 if (v5 < xminval) xminval = v5;
2631 if (v5 > xmaxval) xmaxval = v5;
2632 }
2633
2634 /* construct array of 3rd order absolute differences */
2635 if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
2636 differences[nvals] = abs((2 * v3) - v1 - v5);
2637 nvals++;
2638 } else {
2639 /* ignore constant background regions */
2640 ngoodpix++;
2641 }
2642
2643 /* shift over 1 pixel */
2644 v1 = v2;
2645 v2 = v3;
2646 v3 = v4;
2647 v4 = v5;
2648 } /* end of loop over pixels in the row */
2649
2650 /* compute the 3rd order diffs */
2651 /* Note that there are 4 more pixel values than there are diffs values. */
2652 ngoodpix += (nvals + 4);
2653
2654 if (nvals == 0) {
2655 continue; /* cannot compute medians on this row */
2656 } else if (nvals == 1) {
2657 diffs[nrows] = differences[0];
2658 } else {
2659 /* quick_select returns the median MUCH faster than using qsort */
2660 diffs[nrows] = quick_select_int(differences, nvals);
2661 }
2662
2663 nrows++;
2664 } /* end of loop over rows */
2665
2666 /* compute median of the values for each row */
2667 if (nrows == 0) {
2668 xnoise = 0;
2669 } else if (nrows == 1) {
2670 xnoise = diffs[0];
2671 } else {
2672
2673 qsort(diffs, nrows, sizeof(double), FnCompare_double);
2674 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
2675
2676 FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status);
2677
2678 /* do a 4.5 sigma rejection of outliers */
2679 jj = 0;
2680 sigma = 4.5 * sigma;
2681 for (ii = 0; ii < nrows; ii++) {
2682 if ( fabs(diffs[ii] - xnoise) <= sigma) {
2683 if (jj != ii)
2684 diffs[jj] = diffs[ii];
2685 jj++;
2686 }
2687 }
2688 if (ii != jj)
2689 FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status);
2690 }
2691
2692 if (ngood) *ngood = ngoodpix;
2693 if (minval) *minval = xminval;
2694 if (maxval) *maxval = xmaxval;
2695 if (noise) *noise = 0.6052697 * xnoise;
2696
2697 free(diffs);
2698 free(differences);
2699
2700 return(*status);
2701 }
2702 /*--------------------------------------------------------------------------*/
FnNoise3_float(float * array,long nx,long ny,int nullcheck,float nullvalue,long * ngood,float * minval,float * maxval,double * noise,int * status)2703 static int FnNoise3_float
2704 (float *array, /* 2 dimensional array of image pixels */
2705 long nx, /* number of pixels in each row of the image */
2706 long ny, /* number of rows in the image */
2707 int nullcheck, /* check for null values, if true */
2708 float nullvalue, /* value of null pixels, if nullcheck is true */
2709 /* returned parameters */
2710 long *ngood, /* number of good, non-null pixels? */
2711 float *minval, /* minimum non-null value */
2712 float *maxval, /* maximum non-null value */
2713 double *noise, /* returned R.M.S. value of all non-null pixels */
2714 int *status) /* error status */
2715
2716 /*
2717 Estimate the median and background noise in the input image using 3rd order differences.
2718
2719 The noise in the background of the image is calculated using the 3rd order algorithm
2720 developed for deriving the signal to noise ratio in spectra
2721 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
2722
2723 noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
2724
2725 The returned estimates are the median of the values that are computed for each
2726 row of the image.
2727 */
2728 {
2729 long ii, jj, nrows = 0, nvals, ngoodpix = 0;
2730 float *differences, *rowpix, v1, v2, v3, v4, v5;
2731 float xminval = FLT_MAX, xmaxval = -FLT_MAX;
2732 int do_range = 0;
2733 double *diffs, xnoise = 0;
2734
2735 if (nx < 5) {
2736 /* treat entire array as an image with a single row */
2737 nx = nx * ny;
2738 ny = 1;
2739 }
2740
2741 /* rows must have at least 5 pixels to calc noise, so just calc min, max, ngood */
2742 if (nx < 5) {
2743
2744 for (ii = 0; ii < nx; ii++) {
2745 if (nullcheck && array[ii] == nullvalue)
2746 continue;
2747 else {
2748 if (array[ii] < xminval) xminval = array[ii];
2749 if (array[ii] > xmaxval) xmaxval = array[ii];
2750 ngoodpix++;
2751 }
2752 }
2753 if (minval) *minval = xminval;
2754 if (maxval) *maxval = xmaxval;
2755 if (ngood) *ngood = ngoodpix;
2756 if (noise) *noise = 0.;
2757 return(*status);
2758 }
2759
2760 /* do we need to compute the min and max value? */
2761 if (minval || maxval) do_range = 1;
2762
2763 /* allocate arrays used to compute the median and noise estimates */
2764 if (noise) {
2765 differences = calloc(nx, sizeof(float));
2766 if (!differences) {
2767 *status = MEMORY_ALLOCATION;
2768 return(*status);
2769 }
2770
2771 diffs = calloc(ny, sizeof(double));
2772 if (!diffs) {
2773 free(differences);
2774 *status = MEMORY_ALLOCATION;
2775 return(*status);
2776 }
2777 }
2778
2779 /* loop over each row of the image */
2780 for (jj=0; jj < ny; jj++) {
2781
2782 rowpix = array + (jj * nx); /* point to first pixel in the row */
2783
2784 /***** find the first valid pixel in row */
2785 ii = 0;
2786 if (nullcheck)
2787 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2788
2789 if (ii == nx) continue; /* hit end of row */
2790 v1 = rowpix[ii]; /* store the good pixel value */
2791
2792 if (do_range) {
2793 if (v1 < xminval) xminval = v1;
2794 if (v1 > xmaxval) xmaxval = v1;
2795 }
2796
2797 /***** find the 2nd valid pixel in row (which we will skip over) */
2798 ii++;
2799 if (nullcheck)
2800 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2801
2802 if (ii == nx) continue; /* hit end of row */
2803 v2 = rowpix[ii]; /* store the good pixel value */
2804
2805 if (do_range) {
2806 if (v2 < xminval) xminval = v2;
2807 if (v2 > xmaxval) xmaxval = v2;
2808 }
2809
2810 /***** find the 3rd valid pixel in row */
2811 ii++;
2812 if (nullcheck)
2813 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2814
2815 if (ii == nx) continue; /* hit end of row */
2816 v3 = rowpix[ii]; /* store the good pixel value */
2817
2818 if (do_range) {
2819 if (v3 < xminval) xminval = v3;
2820 if (v3 > xmaxval) xmaxval = v3;
2821 }
2822
2823 /* find the 4nd valid pixel in row (to be skipped) */
2824 ii++;
2825 if (nullcheck)
2826 while (ii < nx && rowpix[ii] == nullvalue) ii++;
2827
2828 if (ii == nx) continue; /* hit end of row */
2829 v4 = rowpix[ii]; /* store the good pixel value */
2830
2831 if (do_range) {
2832 if (v4 < xminval) xminval = v4;
2833 if (v4 > xmaxval) xmaxval = v4;
2834 }
2835
2836 /* now populate the differences arrays */
2837 /* for the remaining pixels in the row */
2838 nvals = 0;
2839 for (ii++; ii < nx; ii++) {
2840
2841 /* find the next valid pixel in row */
2842 if (nullcheck)
2843 while (ii < nx && rowpix[ii] == nullvalue) {
2844 ii++;
2845 }
2846
2847 if (ii == nx) break; /* hit end of row */
2848 v5 = rowpix[ii]; /* store the good pixel value */
2849
2850 if (do_range) {
2851 if (v5 < xminval) xminval = v5;
2852 if (v5 > xmaxval) xmaxval = v5;
2853 }
2854
2855 /* construct array of 3rd order absolute differences */
2856 if (noise) {
2857 if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
2858
2859 differences[nvals] = (float) fabs((2. * v3) - v1 - v5);
2860 nvals++;
2861 } else {
2862 /* ignore constant background regions */
2863 ngoodpix++;
2864 }
2865 } else {
2866 /* just increment the number of non-null pixels */
2867 ngoodpix++;
2868 }
2869
2870 /* shift over 1 pixel */
2871 v1 = v2;
2872 v2 = v3;
2873 v3 = v4;
2874 v4 = v5;
2875 } /* end of loop over pixels in the row */
2876
2877 /* compute the 3rd order diffs */
2878 /* Note that there are 4 more pixel values than there are diffs values. */
2879 ngoodpix += (nvals + 4);
2880
2881 if (noise) {
2882 if (nvals == 0) {
2883 continue; /* cannot compute medians on this row */
2884 } else if (nvals == 1) {
2885 diffs[nrows] = differences[0];
2886 } else {
2887 /* quick_select returns the median MUCH faster than using qsort */
2888 diffs[nrows] = quick_select_float(differences, nvals);
2889 }
2890 }
2891 nrows++;
2892 } /* end of loop over rows */
2893
2894 /* compute median of the values for each row */
2895 if (noise) {
2896 if (nrows == 0) {
2897 xnoise = 0;
2898 } else if (nrows == 1) {
2899 xnoise = diffs[0];
2900 } else {
2901 qsort(diffs, nrows, sizeof(double), FnCompare_double);
2902 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
2903 }
2904 }
2905
2906 if (ngood) *ngood = ngoodpix;
2907 if (minval) *minval = xminval;
2908 if (maxval) *maxval = xmaxval;
2909 if (noise) {
2910 *noise = 0.6052697 * xnoise;
2911 free(diffs);
2912 free(differences);
2913 }
2914
2915 return(*status);
2916 }
2917 /*--------------------------------------------------------------------------*/
FnNoise3_double(double * array,long nx,long ny,int nullcheck,double nullvalue,long * ngood,double * minval,double * maxval,double * noise,int * status)2918 static int FnNoise3_double
2919 (double *array, /* 2 dimensional array of image pixels */
2920 long nx, /* number of pixels in each row of the image */
2921 long ny, /* number of rows in the image */
2922 int nullcheck, /* check for null values, if true */
2923 double nullvalue, /* value of null pixels, if nullcheck is true */
2924 /* returned parameters */
2925 long *ngood, /* number of good, non-null pixels? */
2926 double *minval, /* minimum non-null value */
2927 double *maxval, /* maximum non-null value */
2928 double *noise, /* returned R.M.S. value of all non-null pixels */
2929 int *status) /* error status */
2930
2931 /*
2932 Estimate the median and background noise in the input image using 3rd order differences.
2933
2934 The noise in the background of the image is calculated using the 3rd order algorithm
2935 developed for deriving the signal to noise ratio in spectra
2936 (see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/)
2937
2938 noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2)))
2939
2940 The returned estimates are the median of the values that are computed for each
2941 row of the image.
2942 */
2943 {
2944 long ii, jj, nrows = 0, nvals, ngoodpix = 0;
2945 double *differences, *rowpix, v1, v2, v3, v4, v5;
2946 double xminval = DBL_MAX, xmaxval = -DBL_MAX;
2947 int do_range = 0;
2948 double *diffs, xnoise = 0;
2949
2950 if (nx < 5) {
2951 /* treat entire array as an image with a single row */
2952 nx = nx * ny;
2953 ny = 1;
2954 }
2955
2956 /* rows must have at least 5 pixels */
2957 if (nx < 5) {
2958
2959 for (ii = 0; ii < nx; ii++) {
2960 if (nullcheck && array[ii] == nullvalue)
2961 continue;
2962 else {
2963 if (array[ii] < xminval) xminval = array[ii];
2964 if (array[ii] > xmaxval) xmaxval = array[ii];
2965 ngoodpix++;
2966 }
2967 }
2968 if (minval) *minval = xminval;
2969 if (maxval) *maxval = xmaxval;
2970 if (ngood) *ngood = ngoodpix;
2971 if (noise) *noise = 0.;
2972 return(*status);
2973 }
2974
2975 /* do we need to compute the min and max value? */
2976 if (minval || maxval) do_range = 1;
2977
2978 /* allocate arrays used to compute the median and noise estimates */
2979 if (noise) {
2980 differences = calloc(nx, sizeof(double));
2981 if (!differences) {
2982 *status = MEMORY_ALLOCATION;
2983 return(*status);
2984 }
2985
2986 diffs = calloc(ny, sizeof(double));
2987 if (!diffs) {
2988 free(differences);
2989 *status = MEMORY_ALLOCATION;
2990 return(*status);
2991 }
2992 }
2993
2994 /* loop over each row of the image */
2995 for (jj=0; jj < ny; jj++) {
2996
2997 rowpix = array + (jj * nx); /* point to first pixel in the row */
2998
2999 /***** find the first valid pixel in row */
3000 ii = 0;
3001 if (nullcheck)
3002 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3003
3004 if (ii == nx) continue; /* hit end of row */
3005 v1 = rowpix[ii]; /* store the good pixel value */
3006
3007 if (do_range) {
3008 if (v1 < xminval) xminval = v1;
3009 if (v1 > xmaxval) xmaxval = v1;
3010 }
3011
3012 /***** find the 2nd valid pixel in row (which we will skip over) */
3013 ii++;
3014 if (nullcheck)
3015 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3016
3017 if (ii == nx) continue; /* hit end of row */
3018 v2 = rowpix[ii]; /* store the good pixel value */
3019
3020 if (do_range) {
3021 if (v2 < xminval) xminval = v2;
3022 if (v2 > xmaxval) xmaxval = v2;
3023 }
3024
3025 /***** find the 3rd valid pixel in row */
3026 ii++;
3027 if (nullcheck)
3028 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3029
3030 if (ii == nx) continue; /* hit end of row */
3031 v3 = rowpix[ii]; /* store the good pixel value */
3032
3033 if (do_range) {
3034 if (v3 < xminval) xminval = v3;
3035 if (v3 > xmaxval) xmaxval = v3;
3036 }
3037
3038 /* find the 4nd valid pixel in row (to be skipped) */
3039 ii++;
3040 if (nullcheck)
3041 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3042
3043 if (ii == nx) continue; /* hit end of row */
3044 v4 = rowpix[ii]; /* store the good pixel value */
3045
3046 if (do_range) {
3047 if (v4 < xminval) xminval = v4;
3048 if (v4 > xmaxval) xmaxval = v4;
3049 }
3050
3051 /* now populate the differences arrays */
3052 /* for the remaining pixels in the row */
3053 nvals = 0;
3054 for (ii++; ii < nx; ii++) {
3055
3056 /* find the next valid pixel in row */
3057 if (nullcheck)
3058 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3059
3060 if (ii == nx) break; /* hit end of row */
3061 v5 = rowpix[ii]; /* store the good pixel value */
3062
3063 if (do_range) {
3064 if (v5 < xminval) xminval = v5;
3065 if (v5 > xmaxval) xmaxval = v5;
3066 }
3067
3068 /* construct array of 3rd order absolute differences */
3069 if (noise) {
3070 if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) {
3071
3072 differences[nvals] = fabs((2. * v3) - v1 - v5);
3073 nvals++;
3074 } else {
3075 /* ignore constant background regions */
3076 ngoodpix++;
3077 }
3078 } else {
3079 /* just increment the number of non-null pixels */
3080 ngoodpix++;
3081 }
3082
3083 /* shift over 1 pixel */
3084 v1 = v2;
3085 v2 = v3;
3086 v3 = v4;
3087 v4 = v5;
3088 } /* end of loop over pixels in the row */
3089
3090 /* compute the 3rd order diffs */
3091 /* Note that there are 4 more pixel values than there are diffs values. */
3092 ngoodpix += (nvals + 4);
3093
3094 if (noise) {
3095 if (nvals == 0) {
3096 continue; /* cannot compute medians on this row */
3097 } else if (nvals == 1) {
3098 diffs[nrows] = differences[0];
3099 } else {
3100 /* quick_select returns the median MUCH faster than using qsort */
3101 diffs[nrows] = quick_select_double(differences, nvals);
3102 }
3103 }
3104 nrows++;
3105 } /* end of loop over rows */
3106
3107 /* compute median of the values for each row */
3108 if (noise) {
3109 if (nrows == 0) {
3110 xnoise = 0;
3111 } else if (nrows == 1) {
3112 xnoise = diffs[0];
3113 } else {
3114 qsort(diffs, nrows, sizeof(double), FnCompare_double);
3115 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
3116 }
3117 }
3118
3119 if (ngood) *ngood = ngoodpix;
3120 if (minval) *minval = xminval;
3121 if (maxval) *maxval = xmaxval;
3122 if (noise) {
3123 *noise = 0.6052697 * xnoise;
3124 free(diffs);
3125 free(differences);
3126 }
3127
3128 return(*status);
3129 }
3130 /*--------------------------------------------------------------------------*/
FnNoise1_short(short * array,long nx,long ny,int nullcheck,short nullvalue,double * noise,int * status)3131 static int FnNoise1_short
3132 (short *array, /* 2 dimensional array of image pixels */
3133 long nx, /* number of pixels in each row of the image */
3134 long ny, /* number of rows in the image */
3135 int nullcheck, /* check for null values, if true */
3136 short nullvalue, /* value of null pixels, if nullcheck is true */
3137 /* returned parameters */
3138 double *noise, /* returned R.M.S. value of all non-null pixels */
3139 int *status) /* error status */
3140 /*
3141 Estimate the background noise in the input image using sigma of 1st order differences.
3142
3143 noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
3144
3145 The returned estimate is the median of the values that are computed for each
3146 row of the image.
3147 */
3148 {
3149 int iter;
3150 long ii, jj, kk, nrows = 0, nvals;
3151 short *differences, *rowpix, v1;
3152 double *diffs, xnoise, mean, stdev;
3153
3154 /* rows must have at least 3 pixels to estimate noise */
3155 if (nx < 3) {
3156 *noise = 0;
3157 return(*status);
3158 }
3159
3160 /* allocate arrays used to compute the median and noise estimates */
3161 differences = calloc(nx, sizeof(short));
3162 if (!differences) {
3163 *status = MEMORY_ALLOCATION;
3164 return(*status);
3165 }
3166
3167 diffs = calloc(ny, sizeof(double));
3168 if (!diffs) {
3169 free(differences);
3170 *status = MEMORY_ALLOCATION;
3171 return(*status);
3172 }
3173
3174 /* loop over each row of the image */
3175 for (jj=0; jj < ny; jj++) {
3176
3177 rowpix = array + (jj * nx); /* point to first pixel in the row */
3178
3179 /***** find the first valid pixel in row */
3180 ii = 0;
3181 if (nullcheck)
3182 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3183
3184 if (ii == nx) continue; /* hit end of row */
3185 v1 = rowpix[ii]; /* store the good pixel value */
3186
3187 /* now continue populating the differences arrays */
3188 /* for the remaining pixels in the row */
3189 nvals = 0;
3190 for (ii++; ii < nx; ii++) {
3191
3192 /* find the next valid pixel in row */
3193 if (nullcheck)
3194 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3195
3196 if (ii == nx) break; /* hit end of row */
3197
3198 /* construct array of 1st order differences */
3199 differences[nvals] = v1 - rowpix[ii];
3200
3201 nvals++;
3202 /* shift over 1 pixel */
3203 v1 = rowpix[ii];
3204 } /* end of loop over pixels in the row */
3205
3206 if (nvals < 2)
3207 continue;
3208 else {
3209
3210 FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3211
3212 if (stdev > 0.) {
3213 for (iter = 0; iter < NITER; iter++) {
3214 kk = 0;
3215 for (ii = 0; ii < nvals; ii++) {
3216 if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
3217 if (kk < ii)
3218 differences[kk] = differences[ii];
3219 kk++;
3220 }
3221 }
3222 if (kk == nvals) break;
3223
3224 nvals = kk;
3225 FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3226 }
3227 }
3228
3229 diffs[nrows] = stdev;
3230 nrows++;
3231 }
3232 } /* end of loop over rows */
3233
3234 /* compute median of the values for each row */
3235 if (nrows == 0) {
3236 xnoise = 0;
3237 } else if (nrows == 1) {
3238 xnoise = diffs[0];
3239 } else {
3240 qsort(diffs, nrows, sizeof(double), FnCompare_double);
3241 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
3242 }
3243
3244 *noise = .70710678 * xnoise;
3245
3246 free(diffs);
3247 free(differences);
3248
3249 return(*status);
3250 }
3251 /*--------------------------------------------------------------------------*/
FnNoise1_int(int * array,long nx,long ny,int nullcheck,int nullvalue,double * noise,int * status)3252 static int FnNoise1_int
3253 (int *array, /* 2 dimensional array of image pixels */
3254 long nx, /* number of pixels in each row of the image */
3255 long ny, /* number of rows in the image */
3256 int nullcheck, /* check for null values, if true */
3257 int nullvalue, /* value of null pixels, if nullcheck is true */
3258 /* returned parameters */
3259 double *noise, /* returned R.M.S. value of all non-null pixels */
3260 int *status) /* error status */
3261 /*
3262 Estimate the background noise in the input image using sigma of 1st order differences.
3263
3264 noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
3265
3266 The returned estimate is the median of the values that are computed for each
3267 row of the image.
3268 */
3269 {
3270 int iter;
3271 long ii, jj, kk, nrows = 0, nvals;
3272 int *differences, *rowpix, v1;
3273 double *diffs, xnoise, mean, stdev;
3274
3275 /* rows must have at least 3 pixels to estimate noise */
3276 if (nx < 3) {
3277 *noise = 0;
3278 return(*status);
3279 }
3280
3281 /* allocate arrays used to compute the median and noise estimates */
3282 differences = calloc(nx, sizeof(int));
3283 if (!differences) {
3284 *status = MEMORY_ALLOCATION;
3285 return(*status);
3286 }
3287
3288 diffs = calloc(ny, sizeof(double));
3289 if (!diffs) {
3290 free(differences);
3291 *status = MEMORY_ALLOCATION;
3292 return(*status);
3293 }
3294
3295 /* loop over each row of the image */
3296 for (jj=0; jj < ny; jj++) {
3297
3298 rowpix = array + (jj * nx); /* point to first pixel in the row */
3299
3300 /***** find the first valid pixel in row */
3301 ii = 0;
3302 if (nullcheck)
3303 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3304
3305 if (ii == nx) continue; /* hit end of row */
3306 v1 = rowpix[ii]; /* store the good pixel value */
3307
3308 /* now continue populating the differences arrays */
3309 /* for the remaining pixels in the row */
3310 nvals = 0;
3311 for (ii++; ii < nx; ii++) {
3312
3313 /* find the next valid pixel in row */
3314 if (nullcheck)
3315 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3316
3317 if (ii == nx) break; /* hit end of row */
3318
3319 /* construct array of 1st order differences */
3320 differences[nvals] = v1 - rowpix[ii];
3321
3322 nvals++;
3323 /* shift over 1 pixel */
3324 v1 = rowpix[ii];
3325 } /* end of loop over pixels in the row */
3326
3327 if (nvals < 2)
3328 continue;
3329 else {
3330
3331 FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3332
3333 if (stdev > 0.) {
3334 for (iter = 0; iter < NITER; iter++) {
3335 kk = 0;
3336 for (ii = 0; ii < nvals; ii++) {
3337 if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
3338 if (kk < ii)
3339 differences[kk] = differences[ii];
3340 kk++;
3341 }
3342 }
3343 if (kk == nvals) break;
3344
3345 nvals = kk;
3346 FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3347 }
3348 }
3349
3350 diffs[nrows] = stdev;
3351 nrows++;
3352 }
3353 } /* end of loop over rows */
3354
3355 /* compute median of the values for each row */
3356 if (nrows == 0) {
3357 xnoise = 0;
3358 } else if (nrows == 1) {
3359 xnoise = diffs[0];
3360 } else {
3361 qsort(diffs, nrows, sizeof(double), FnCompare_double);
3362 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
3363 }
3364
3365 *noise = .70710678 * xnoise;
3366
3367 free(diffs);
3368 free(differences);
3369
3370 return(*status);
3371 }
3372 /*--------------------------------------------------------------------------*/
FnNoise1_float(float * array,long nx,long ny,int nullcheck,float nullvalue,double * noise,int * status)3373 static int FnNoise1_float
3374 (float *array, /* 2 dimensional array of image pixels */
3375 long nx, /* number of pixels in each row of the image */
3376 long ny, /* number of rows in the image */
3377 int nullcheck, /* check for null values, if true */
3378 float nullvalue, /* value of null pixels, if nullcheck is true */
3379 /* returned parameters */
3380 double *noise, /* returned R.M.S. value of all non-null pixels */
3381 int *status) /* error status */
3382 /*
3383 Estimate the background noise in the input image using sigma of 1st order differences.
3384
3385 noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
3386
3387 The returned estimate is the median of the values that are computed for each
3388 row of the image.
3389 */
3390 {
3391 int iter;
3392 long ii, jj, kk, nrows = 0, nvals;
3393 float *differences, *rowpix, v1;
3394 double *diffs, xnoise, mean, stdev;
3395
3396 /* rows must have at least 3 pixels to estimate noise */
3397 if (nx < 3) {
3398 *noise = 0;
3399 return(*status);
3400 }
3401
3402 /* allocate arrays used to compute the median and noise estimates */
3403 differences = calloc(nx, sizeof(float));
3404 if (!differences) {
3405 *status = MEMORY_ALLOCATION;
3406 return(*status);
3407 }
3408
3409 diffs = calloc(ny, sizeof(double));
3410 if (!diffs) {
3411 free(differences);
3412 *status = MEMORY_ALLOCATION;
3413 return(*status);
3414 }
3415
3416 /* loop over each row of the image */
3417 for (jj=0; jj < ny; jj++) {
3418
3419 rowpix = array + (jj * nx); /* point to first pixel in the row */
3420
3421 /***** find the first valid pixel in row */
3422 ii = 0;
3423 if (nullcheck)
3424 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3425
3426 if (ii == nx) continue; /* hit end of row */
3427 v1 = rowpix[ii]; /* store the good pixel value */
3428
3429 /* now continue populating the differences arrays */
3430 /* for the remaining pixels in the row */
3431 nvals = 0;
3432 for (ii++; ii < nx; ii++) {
3433
3434 /* find the next valid pixel in row */
3435 if (nullcheck)
3436 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3437
3438 if (ii == nx) break; /* hit end of row */
3439
3440 /* construct array of 1st order differences */
3441 differences[nvals] = v1 - rowpix[ii];
3442
3443 nvals++;
3444 /* shift over 1 pixel */
3445 v1 = rowpix[ii];
3446 } /* end of loop over pixels in the row */
3447
3448 if (nvals < 2)
3449 continue;
3450 else {
3451
3452 FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3453
3454 if (stdev > 0.) {
3455 for (iter = 0; iter < NITER; iter++) {
3456 kk = 0;
3457 for (ii = 0; ii < nvals; ii++) {
3458 if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
3459 if (kk < ii)
3460 differences[kk] = differences[ii];
3461 kk++;
3462 }
3463 }
3464 if (kk == nvals) break;
3465
3466 nvals = kk;
3467 FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3468 }
3469 }
3470
3471 diffs[nrows] = stdev;
3472 nrows++;
3473 }
3474 } /* end of loop over rows */
3475
3476 /* compute median of the values for each row */
3477 if (nrows == 0) {
3478 xnoise = 0;
3479 } else if (nrows == 1) {
3480 xnoise = diffs[0];
3481 } else {
3482 qsort(diffs, nrows, sizeof(double), FnCompare_double);
3483 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
3484 }
3485
3486 *noise = .70710678 * xnoise;
3487
3488 free(diffs);
3489 free(differences);
3490
3491 return(*status);
3492 }
3493 /*--------------------------------------------------------------------------*/
FnNoise1_double(double * array,long nx,long ny,int nullcheck,double nullvalue,double * noise,int * status)3494 static int FnNoise1_double
3495 (double *array, /* 2 dimensional array of image pixels */
3496 long nx, /* number of pixels in each row of the image */
3497 long ny, /* number of rows in the image */
3498 int nullcheck, /* check for null values, if true */
3499 double nullvalue, /* value of null pixels, if nullcheck is true */
3500 /* returned parameters */
3501 double *noise, /* returned R.M.S. value of all non-null pixels */
3502 int *status) /* error status */
3503 /*
3504 Estimate the background noise in the input image using sigma of 1st order differences.
3505
3506 noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1])
3507
3508 The returned estimate is the median of the values that are computed for each
3509 row of the image.
3510 */
3511 {
3512 int iter;
3513 long ii, jj, kk, nrows = 0, nvals;
3514 double *differences, *rowpix, v1;
3515 double *diffs, xnoise, mean, stdev;
3516
3517 /* rows must have at least 3 pixels to estimate noise */
3518 if (nx < 3) {
3519 *noise = 0;
3520 return(*status);
3521 }
3522
3523 /* allocate arrays used to compute the median and noise estimates */
3524 differences = calloc(nx, sizeof(double));
3525 if (!differences) {
3526 *status = MEMORY_ALLOCATION;
3527 return(*status);
3528 }
3529
3530 diffs = calloc(ny, sizeof(double));
3531 if (!diffs) {
3532 free(differences);
3533 *status = MEMORY_ALLOCATION;
3534 return(*status);
3535 }
3536
3537 /* loop over each row of the image */
3538 for (jj=0; jj < ny; jj++) {
3539
3540 rowpix = array + (jj * nx); /* point to first pixel in the row */
3541
3542 /***** find the first valid pixel in row */
3543 ii = 0;
3544 if (nullcheck)
3545 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3546
3547 if (ii == nx) continue; /* hit end of row */
3548 v1 = rowpix[ii]; /* store the good pixel value */
3549
3550 /* now continue populating the differences arrays */
3551 /* for the remaining pixels in the row */
3552 nvals = 0;
3553 for (ii++; ii < nx; ii++) {
3554
3555 /* find the next valid pixel in row */
3556 if (nullcheck)
3557 while (ii < nx && rowpix[ii] == nullvalue) ii++;
3558
3559 if (ii == nx) break; /* hit end of row */
3560
3561 /* construct array of 1st order differences */
3562 differences[nvals] = v1 - rowpix[ii];
3563
3564 nvals++;
3565 /* shift over 1 pixel */
3566 v1 = rowpix[ii];
3567 } /* end of loop over pixels in the row */
3568
3569 if (nvals < 2)
3570 continue;
3571 else {
3572
3573 FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3574
3575 if (stdev > 0.) {
3576 for (iter = 0; iter < NITER; iter++) {
3577 kk = 0;
3578 for (ii = 0; ii < nvals; ii++) {
3579 if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) {
3580 if (kk < ii)
3581 differences[kk] = differences[ii];
3582 kk++;
3583 }
3584 }
3585 if (kk == nvals) break;
3586
3587 nvals = kk;
3588 FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status);
3589 }
3590 }
3591
3592 diffs[nrows] = stdev;
3593 nrows++;
3594 }
3595 } /* end of loop over rows */
3596
3597 /* compute median of the values for each row */
3598 if (nrows == 0) {
3599 xnoise = 0;
3600 } else if (nrows == 1) {
3601 xnoise = diffs[0];
3602 } else {
3603 qsort(diffs, nrows, sizeof(double), FnCompare_double);
3604 xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.;
3605 }
3606
3607 *noise = .70710678 * xnoise;
3608
3609 free(diffs);
3610 free(differences);
3611
3612 return(*status);
3613 }
3614 /*--------------------------------------------------------------------------*/
FnCompare_short(const void * v1,const void * v2)3615 static int FnCompare_short(const void *v1, const void *v2)
3616 {
3617 const short *i1 = v1;
3618 const short *i2 = v2;
3619
3620 if (*i1 < *i2)
3621 return(-1);
3622 else if (*i1 > *i2)
3623 return(1);
3624 else
3625 return(0);
3626 }
3627 /*--------------------------------------------------------------------------*/
FnCompare_int(const void * v1,const void * v2)3628 static int FnCompare_int(const void *v1, const void *v2)
3629 {
3630 const int *i1 = v1;
3631 const int *i2 = v2;
3632
3633 if (*i1 < *i2)
3634 return(-1);
3635 else if (*i1 > *i2)
3636 return(1);
3637 else
3638 return(0);
3639 }
3640 /*--------------------------------------------------------------------------*/
FnCompare_float(const void * v1,const void * v2)3641 static int FnCompare_float(const void *v1, const void *v2)
3642 {
3643 const float *i1 = v1;
3644 const float *i2 = v2;
3645
3646 if (*i1 < *i2)
3647 return(-1);
3648 else if (*i1 > *i2)
3649 return(1);
3650 else
3651 return(0);
3652 }
3653 /*--------------------------------------------------------------------------*/
FnCompare_double(const void * v1,const void * v2)3654 static int FnCompare_double(const void *v1, const void *v2)
3655 {
3656 const double *i1 = v1;
3657 const double *i2 = v2;
3658
3659 if (*i1 < *i2)
3660 return(-1);
3661 else if (*i1 > *i2)
3662 return(1);
3663 else
3664 return(0);
3665 }
3666 /*--------------------------------------------------------------------------*/
3667
3668 /*
3669 * These Quickselect routines are based on the algorithm described in
3670 * "Numerical recipes in C", Second Edition,
3671 * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
3672 * This code by Nicolas Devillard - 1998. Public domain.
3673 */
3674
3675 /*--------------------------------------------------------------------------*/
3676
3677 #define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
3678
quick_select_float(float arr[],int n)3679 static float quick_select_float(float arr[], int n)
3680 {
3681 int low, high ;
3682 int median;
3683 int middle, ll, hh;
3684
3685 low = 0 ; high = n-1 ; median = (low + high) / 2;
3686 for (;;) {
3687 if (high <= low) /* One element only */
3688 return arr[median] ;
3689
3690 if (high == low + 1) { /* Two elements only */
3691 if (arr[low] > arr[high])
3692 ELEM_SWAP(arr[low], arr[high]) ;
3693 return arr[median] ;
3694 }
3695
3696 /* Find median of low, middle and high items; swap into position low */
3697 middle = (low + high) / 2;
3698 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
3699 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
3700 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
3701
3702 /* Swap low item (now in position middle) into position (low+1) */
3703 ELEM_SWAP(arr[middle], arr[low+1]) ;
3704
3705 /* Nibble from each end towards middle, swapping items when stuck */
3706 ll = low + 1;
3707 hh = high;
3708 for (;;) {
3709 do ll++; while (arr[low] > arr[ll]) ;
3710 do hh--; while (arr[hh] > arr[low]) ;
3711
3712 if (hh < ll)
3713 break;
3714
3715 ELEM_SWAP(arr[ll], arr[hh]) ;
3716 }
3717
3718 /* Swap middle item (in position low) back into correct position */
3719 ELEM_SWAP(arr[low], arr[hh]) ;
3720
3721 /* Re-set active partition */
3722 if (hh <= median)
3723 low = ll;
3724 if (hh >= median)
3725 high = hh - 1;
3726 }
3727 }
3728
3729 #undef ELEM_SWAP
3730
3731 /*--------------------------------------------------------------------------*/
3732
3733 #define ELEM_SWAP(a,b) { register short t=(a);(a)=(b);(b)=t; }
3734
quick_select_short(short arr[],int n)3735 static short quick_select_short(short arr[], int n)
3736 {
3737 int low, high ;
3738 int median;
3739 int middle, ll, hh;
3740
3741 low = 0 ; high = n-1 ; median = (low + high) / 2;
3742 for (;;) {
3743 if (high <= low) /* One element only */
3744 return arr[median] ;
3745
3746 if (high == low + 1) { /* Two elements only */
3747 if (arr[low] > arr[high])
3748 ELEM_SWAP(arr[low], arr[high]) ;
3749 return arr[median] ;
3750 }
3751
3752 /* Find median of low, middle and high items; swap into position low */
3753 middle = (low + high) / 2;
3754 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
3755 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
3756 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
3757
3758 /* Swap low item (now in position middle) into position (low+1) */
3759 ELEM_SWAP(arr[middle], arr[low+1]) ;
3760
3761 /* Nibble from each end towards middle, swapping items when stuck */
3762 ll = low + 1;
3763 hh = high;
3764 for (;;) {
3765 do ll++; while (arr[low] > arr[ll]) ;
3766 do hh--; while (arr[hh] > arr[low]) ;
3767
3768 if (hh < ll)
3769 break;
3770
3771 ELEM_SWAP(arr[ll], arr[hh]) ;
3772 }
3773
3774 /* Swap middle item (in position low) back into correct position */
3775 ELEM_SWAP(arr[low], arr[hh]) ;
3776
3777 /* Re-set active partition */
3778 if (hh <= median)
3779 low = ll;
3780 if (hh >= median)
3781 high = hh - 1;
3782 }
3783 }
3784
3785 #undef ELEM_SWAP
3786
3787 /*--------------------------------------------------------------------------*/
3788
3789 #define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; }
3790
quick_select_int(int arr[],int n)3791 static int quick_select_int(int arr[], int n)
3792 {
3793 int low, high ;
3794 int median;
3795 int middle, ll, hh;
3796
3797 low = 0 ; high = n-1 ; median = (low + high) / 2;
3798 for (;;) {
3799 if (high <= low) /* One element only */
3800 return arr[median] ;
3801
3802 if (high == low + 1) { /* Two elements only */
3803 if (arr[low] > arr[high])
3804 ELEM_SWAP(arr[low], arr[high]) ;
3805 return arr[median] ;
3806 }
3807
3808 /* Find median of low, middle and high items; swap into position low */
3809 middle = (low + high) / 2;
3810 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
3811 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
3812 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
3813
3814 /* Swap low item (now in position middle) into position (low+1) */
3815 ELEM_SWAP(arr[middle], arr[low+1]) ;
3816
3817 /* Nibble from each end towards middle, swapping items when stuck */
3818 ll = low + 1;
3819 hh = high;
3820 for (;;) {
3821 do ll++; while (arr[low] > arr[ll]) ;
3822 do hh--; while (arr[hh] > arr[low]) ;
3823
3824 if (hh < ll)
3825 break;
3826
3827 ELEM_SWAP(arr[ll], arr[hh]) ;
3828 }
3829
3830 /* Swap middle item (in position low) back into correct position */
3831 ELEM_SWAP(arr[low], arr[hh]) ;
3832
3833 /* Re-set active partition */
3834 if (hh <= median)
3835 low = ll;
3836 if (hh >= median)
3837 high = hh - 1;
3838 }
3839 }
3840
3841 #undef ELEM_SWAP
3842
3843 /*--------------------------------------------------------------------------*/
3844
3845 #define ELEM_SWAP(a,b) { register LONGLONG t=(a);(a)=(b);(b)=t; }
3846
quick_select_longlong(LONGLONG arr[],int n)3847 static LONGLONG quick_select_longlong(LONGLONG arr[], int n)
3848 {
3849 int low, high ;
3850 int median;
3851 int middle, ll, hh;
3852
3853 low = 0 ; high = n-1 ; median = (low + high) / 2;
3854 for (;;) {
3855 if (high <= low) /* One element only */
3856 return arr[median] ;
3857
3858 if (high == low + 1) { /* Two elements only */
3859 if (arr[low] > arr[high])
3860 ELEM_SWAP(arr[low], arr[high]) ;
3861 return arr[median] ;
3862 }
3863
3864 /* Find median of low, middle and high items; swap into position low */
3865 middle = (low + high) / 2;
3866 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
3867 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
3868 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
3869
3870 /* Swap low item (now in position middle) into position (low+1) */
3871 ELEM_SWAP(arr[middle], arr[low+1]) ;
3872
3873 /* Nibble from each end towards middle, swapping items when stuck */
3874 ll = low + 1;
3875 hh = high;
3876 for (;;) {
3877 do ll++; while (arr[low] > arr[ll]) ;
3878 do hh--; while (arr[hh] > arr[low]) ;
3879
3880 if (hh < ll)
3881 break;
3882
3883 ELEM_SWAP(arr[ll], arr[hh]) ;
3884 }
3885
3886 /* Swap middle item (in position low) back into correct position */
3887 ELEM_SWAP(arr[low], arr[hh]) ;
3888
3889 /* Re-set active partition */
3890 if (hh <= median)
3891 low = ll;
3892 if (hh >= median)
3893 high = hh - 1;
3894 }
3895 }
3896
3897 #undef ELEM_SWAP
3898
3899 /*--------------------------------------------------------------------------*/
3900
3901 #define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; }
3902
quick_select_double(double arr[],int n)3903 static double quick_select_double(double arr[], int n)
3904 {
3905 int low, high ;
3906 int median;
3907 int middle, ll, hh;
3908
3909 low = 0 ; high = n-1 ; median = (low + high) / 2;
3910 for (;;) {
3911 if (high <= low) /* One element only */
3912 return arr[median] ;
3913
3914 if (high == low + 1) { /* Two elements only */
3915 if (arr[low] > arr[high])
3916 ELEM_SWAP(arr[low], arr[high]) ;
3917 return arr[median] ;
3918 }
3919
3920 /* Find median of low, middle and high items; swap into position low */
3921 middle = (low + high) / 2;
3922 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
3923 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
3924 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
3925
3926 /* Swap low item (now in position middle) into position (low+1) */
3927 ELEM_SWAP(arr[middle], arr[low+1]) ;
3928
3929 /* Nibble from each end towards middle, swapping items when stuck */
3930 ll = low + 1;
3931 hh = high;
3932 for (;;) {
3933 do ll++; while (arr[low] > arr[ll]) ;
3934 do hh--; while (arr[hh] > arr[low]) ;
3935
3936 if (hh < ll)
3937 break;
3938
3939 ELEM_SWAP(arr[ll], arr[hh]) ;
3940 }
3941
3942 /* Swap middle item (in position low) back into correct position */
3943 ELEM_SWAP(arr[low], arr[hh]) ;
3944
3945 /* Re-set active partition */
3946 if (hh <= median)
3947 low = ll;
3948 if (hh >= median)
3949 high = hh - 1;
3950 }
3951 }
3952
3953 #undef ELEM_SWAP
3954
3955
3956