version update
[goodguy/cinelerra.git] / cinelerra-5.1 / db / utils / histeq.C
1 #include<stdio.h>
2 #include<stdint.h>
3 #include<stdlib.h>
4 #include<string.h>
5 #include<math.h>
6 #include "../guicast/linklist.h"
7
8
9
10 class HistEq {
11         void fit(int *dat, int n, double &a, double &b);
12 public:
13         void eq(unsigned char *dp, int len);
14         HistEq() {}
15         ~HistEq() {}
16 };
17
18 double mean(uint8_t *dat, int n)
19 {
20   int s = 0;
21   for( int i=0; i<n; ++i ) s += *dat++;
22   return (double)s / n;
23 }
24
25 double variance(uint8_t *dat, int n, double m=-1)
26 {
27   if( m < 0 ) m = mean(dat, n);
28   double ss = 0, s = 0;
29   for( int i=0; i<n; ++i ) {
30     double dx = dat[i]-m;
31     s += dx;  ss += dx*dx;
32   }
33   return (ss - s*s / n) / (n-1);
34 }
35
36 double std_dev(uint8_t *dat, int n, double m=-1)
37 {
38   return sqrt(variance(dat, n, m));
39 }
40
41 static inline int clip(int v, int mn, int mx)
42 {
43   return v<mn ? mn : v>mx ? mx : v;
44 }
45
46 void HistEq::
47 fit(int *dat, int n, double &a, double &b)
48 {
49   double st2 = 0, sb = 0;
50   int64_t sy = 0;
51   for( int i=0; i<n; ++i ) sy += dat[i];
52   double mx = (n-1)/2., my = (double)sy / n;
53   for( int i=0; i<n; ++i ) {
54     double t = i - mx;
55     st2 += t * t;
56     sb += t * dat[i];
57   }
58   b = sb / st2;
59   a = my - mx*b;
60 }
61
62 void HistEq::
63 eq(uint8_t *dp, int len)
64 {
65   int hist[256], wts[256], map[256];
66   for( int i=0; i<256; ++i ) hist[i] = 0;
67   uint8_t *bp = dp;
68   for( int i=len; --i>=0; ++bp ) ++hist[*bp];
69   int t = 0;
70   for( int i=0; i<256; ++i ) { t += hist[i];  wts[i] = t; }
71   int lmn = len/20, lmx = len-lmn;
72   int mn =   0;  while( mn < 256 && wts[mn] < lmn ) ++mn;
73   int mx = 255;  while( mx > mn && wts[mx] > lmx ) --mx;
74   double a, b;
75   fit(&wts[mn], mx-mn, a, b);
76   double r = 256./len;
77   double a1 = (a - b*mn) * r, b1 = b * r;
78   for( int i=0 ; i<256;  ++i ) map[i] = clip(a1 + i*b1, 0, 255);
79   bp = dp;
80   for( int i=len; --i>=0; ++bp ) *bp = map[*bp];
81 }
82
83
84 int main(int ac, char **av)
85 {
86   FILE *ifp = !strcmp(av[1],"-") ? stdin : fopen(av[1],"r");
87   FILE *ofp = !strcmp(av[2],"-") ? stdout : fopen(av[2],"w");
88   char line[120];
89   fgets(line,sizeof(line),ifp);  fputs(line,ofp);
90   fgets(line,sizeof(line),ifp);  fputs(line,ofp);
91   int w, h;  if( sscanf(line,"%d %d\n",&w,&h) != 2 ) exit(1);
92   fgets(line,sizeof(line),ifp);  fputs(line,ofp);
93   int len = w*h;
94   unsigned char data[len], *bp = data;
95   for( int ch, i=len; --i>=0 && (ch=getc(ifp)) >= 0; ++bp ) *bp = ch;
96   HistEq().eq(data,len);
97   bp = data;
98   for( int i=len; --i>=0; ++bp ) putc(*bp,ofp);
99
100   if( ifp != stdin ) fclose(ifp);
101   if( ofp != stdout ) fclose(ofp);
102   return 0;
103 }
104