source: tags/release-1.0.5/src/Utils/zscale.cc @ 1455

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

Generally tidying up code and adding more useful warning/error messages.
Some improvement of the constructor and save array tasks in cubes.cc,
adding tests for negative sizes and changes in array size.
File size: 3.0 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);
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  float *smallarray = new float[nsample];
75  float *ct = new float[nsample];
76  long size=0;
77
78  float mean=0.,sd=0.;
79  for(int i=0;i<newsize;i++) mean+=newimage[i];
80  mean /= float(newsize);
81  for(int i=0;i<newsize;i++) sd+=(newimage[i]-mean)*(newimage[i]-mean);
82  sd /= float(newsize);
83 
84
85  if(newsize>nsample){
86    int step = (newsize / nsample) + 1;
87    for(int i=0;i<newsize;i++){
88      if((i%step)==0){
89        smallarray[size] = newimage[i];
90        ct[size]=(float)size+1.;
91        size++;
92      }
93    }
94  }
95  else{
96    for(int i=0;i<newsize;i++){
97      smallarray[i] = newimage[i];
98      ct[i] = float(i+1);
99    }
100    size=newsize;
101  }
102
103  sort(smallarray,0,size);
104 
105  /* fit a linear slope to the centre of the cumulative distribution */
106  long midpt = size/2;
107  long imin = midpt - width;
108  float slope,intercept,errSlope,errIntercept,r;
109  if(size<2*width)
110    linear_regression(size,ct,smallarray,0,size-1,slope,errSlope,
111                      intercept,errIntercept,r);
112  else linear_regression(size,ct,smallarray,imin,imin+2*width,slope,errSlope,
113                    intercept,errIntercept,r);
114
115  z1 = smallarray[midpt] + (slope/contrast)*(float)(1-midpt);
116  z2 = smallarray[midpt] + (slope/contrast)*(float)(nsample-midpt);
117
118  if(z1==z2){
119 
120    if(z1!=0) z2 = z1 * 1.05;
121    else z2 = z1+0.01;
122
123  }
124
125  delete [] newimage;
126  delete [] smallarray;
127  delete [] ct;
128
129           
130}
131
132
Note: See TracBrowser for help on using the repository browser.