source: trunk/Utils/zscale.cc @ 3

Last change on this file since 3 was 3, checked in by Matthew Whiting, 18 years ago

This is the first full import of all working code to
the Duchamp repository.
Made three directories at top level:

branches/ tags/ trunk/

and trunk/ has the full set of code:
ATrous/ Cubes/ Detection/ InputComplete? InputExample? README Utils/ docs/ mainDuchamp.cc param.cc param.hh
File size: 3.1 KB
Line 
1#include <iostream>
2#include <Utils/utils.hh>
3const int nsample=1000;
4// const int width=300;
5const int width=150;
6const float contrast=0.25;
7
8void zscale(long imagesize, float *image, float &z1, float &z2)
9{
10
11  float *smallarray = new float[nsample];
12  float *ct = new float[nsample];
13  long size=0;
14
15  float mean=0.,sd=0.;
16  for(int i=0;i<imagesize;i++) mean+=image[i];
17  mean /= float(imagesize);
18  for(int i=0;i<imagesize;i++) sd+=(image[i]-mean)*(image[i]-mean);
19  sd /= float(imagesize);
20 
21
22  if(imagesize>nsample){
23    int step = (imagesize / nsample) + 1;
24    for(int i=0;i<imagesize;i++){
25      if((i%step)==0){
26        smallarray[size] = image[i];
27        ct[size]=(float)size+1.;
28        size++;
29      }
30    }
31  }
32  else{
33    for(int i=0;i<imagesize;i++){
34      smallarray[i] = image[i];
35      ct[i] = float(i+1);
36    }
37    size=imagesize;
38  }
39
40  sort(smallarray,0,size-1);
41 
42  /* fit a linear slope to the centre of the cumulative distribution */
43  long midpt = size/2;
44  long imin = midpt - width;
45  float slope,intercept,errSlope,errIntercept,r;
46  if(size<2*width)
47    linear_regression(size,ct,smallarray,0,size-1,slope,errSlope,
48                      intercept,errIntercept,r);
49  else linear_regression(size,ct,smallarray,imin,imin+2*width,slope,errSlope,
50                         intercept,errIntercept,r);
51
52  z1 = smallarray[midpt] + (slope/contrast)*(float)(1-midpt);
53  z2 = smallarray[midpt] + (slope/contrast)*(float)(nsample-midpt);
54
55  if(z1==z2){
56 
57    if(z1!=0) z2 = z1 * 1.05;
58    else z2 = z1+0.01;
59
60  }
61
62  delete [] smallarray;
63  delete [] ct;
64           
65}
66
67
68void zscale(long imagesize, float *image, float &z1, float &z2, float blankVal)
69{
70  float *newimage = new float[imagesize];
71  int newsize=0;
72  for(int i=0;i<imagesize;i++) if(image[i]!=blankVal) newimage[newsize++] = image[i];
73
74//   cerr<<"Sizes: "<<imagesize<<"   "<<newsize<<endl;
75
76  float *smallarray = new float[nsample];
77  float *ct = new float[nsample];
78  long size=0;
79
80  float mean=0.,sd=0.;
81  for(int i=0;i<newsize;i++) mean+=newimage[i];
82  mean /= float(newsize);
83  for(int i=0;i<newsize;i++) sd+=(newimage[i]-mean)*(newimage[i]-mean);
84  sd /= float(newsize);
85 
86
87  if(newsize>nsample){
88    int step = (newsize / nsample) + 1;
89    for(int i=0;i<newsize;i++){
90      if((i%step)==0){
91        smallarray[size] = newimage[i];
92        ct[size]=(float)size+1.;
93        size++;
94      }
95    }
96  }
97  else{
98    for(int i=0;i<newsize;i++){
99      smallarray[i] = newimage[i];
100      ct[i] = float(i+1);
101    }
102    size=newsize;
103  }
104
105  sort(smallarray,0,size-1);
106 
107  /* fit a linear slope to the centre of the cumulative distribution */
108  long midpt = size/2;
109  long imin = midpt - width;
110  float slope,intercept,errSlope,errIntercept,r;
111  if(size<2*width)
112    linear_regression(size,ct,smallarray,0,size-1,slope,errSlope,
113                      intercept,errIntercept,r);
114  else linear_regression(size,ct,smallarray,imin,imin+2*width,slope,errSlope,
115                    intercept,errIntercept,r);
116
117  z1 = smallarray[midpt] + (slope/contrast)*(float)(1-midpt);
118  z2 = smallarray[midpt] + (slope/contrast)*(float)(nsample-midpt);
119
120  if(z1==z2){
121 
122    if(z1!=0) z2 = z1 * 1.05;
123    else z2 = z1+0.01;
124
125  }
126
127  delete [] newimage;
128  delete [] smallarray;
129  delete [] ct;
130
131           
132}
133
134
Note: See TracBrowser for help on using the repository browser.