source: tags/release-0.9/Utils/zscale.cc @ 813

Last change on this file since 813 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.