File: tone.c

package info (click to toggle)
fitspng 1.4-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 1,804 kB
  • sloc: sh: 1,232; ansic: 860; makefile: 14
file content (122 lines) | stat: -rw-r--r-- 3,254 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
/*

  FITSPNG  --  eCDF base estimation of intensity scaling parameters
  Copyright (C) 2019  Filip Hroch, Masaryk University, Brno, CZ

  Fitspng is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  Fitspng is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Fitspng.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "fitspng.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <string.h>

#include <assert.h>

#define QSENSE 0.95
#define NMAX 32768
#define MAX(x,y) (x) > (y) ? (x) : (y);


int tone(long npix, const float *pic, float qblack, float rsense,
	 float *thresh, float *slope, int verb)
{
  int skip = MAX(npix / NMAX,1);
  long nmax = npix / skip;
  long nbytes = (nmax+1)*sizeof(float);

  float *d, *xcdf, *ycdf;
  if( (d = malloc(nbytes)) == NULL ||
      (xcdf = malloc(nbytes)) == NULL || (ycdf = malloc(nbytes)) == NULL ) {
    perror("There is no room for determination of an empirical CDF");
    return 0;
  }

  /* PART I. Estimate of background level */
  long m = 0;
  for(long i = 0; i < npix - skip && m < nmax; i += skip)
    d[m++] = pic[i];
  int ncdf = ecdf(m,d,xcdf,ycdf);

  float black = quantile(ncdf,xcdf,ycdf,qblack);

  float median = quantile(ncdf,xcdf,ycdf,0.5);
  float q75 = quantile(ncdf,xcdf,ycdf,0.75);
  float q25 = quantile(ncdf,xcdf,ycdf,0.25);
  float mad = (q75 - q25) / 2;
  float threshold = 3 * mad / 0.6745; /* 3-sigma rule */

  /* Grapes */
  /*
  FILE *file = fopen("b.dat","w");
  for(int i = 0; i < ncdf; i++)
    fprintf(file,"%f %f\n",xcdf[i],ycdf[i]);
  fclose(file);

  int nhist = 65536;
  int *hist;
  hist = malloc(nhist*sizeof(int));
  for(int i = 0; i < nhist; i++ ) hist[i] = 0;
  for(int i = 0; i < npix; i++) {
    int n = pic[i];
    if( 0 <= n && n < nhist )
      hist[n] = hist[n] + 1;
  }
  file = fopen("bh.dat","w");
  for(int i = 0; i < nhist; i++)
    fprintf(file,"%d %d\n",i,hist[i]);
  fclose(file);
  */

  /* PART II. Estimate of starlight */
  long n = 0;
  long side = 1;
  while( n < (NMAX / 10) && side <= 16 ) {
    side = 2*side;
    long skip2 = MAX(skip / side,1);
    long imax = npix - skip2;
    for(long i = 0; n < nmax && i < imax; i += skip2) {
      float r = pic[i] - median;
      if( r > threshold )
	d[n++] = r;
    }
  }

  float q1  = quantile(ncdf,xcdf,ycdf,1-QSENSE);
  float q99 = quantile(ncdf,xcdf,ycdf,QSENSE);
  float refsense = black > 0 ? sqrt(black) : 1;
  if( n > 0 ) {
    ncdf = ecdf(n,d,xcdf,ycdf);
    refsense = quantile(ncdf,xcdf,ycdf,QSENSE);
  }
  else if ( median > 100 * FLT_EPSILON )
    refsense = 3 * (q99 - q1);

  free(d);
  free(xcdf);
  free(ycdf);

  if( verb )
    fprintf(stderr,"Median=%f MAD=%f Black=%f Qlight=%f n=%ld\n",
	    median, mad, black, refsense, n);


  *thresh = black;
  *slope = refsense / rsense;

  return 1;
}