source: trunk/src/Utils/zscale.cc @ 258

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

Improved the setCubeStats function and changed the sorting function over to the STL sort and stable_sort function. Summary:

  • Utils/getStats.cc : changed calls to sort to std::sort
  • Utils/utils.hh : made swap function an inline one.
  • Utils/zscale.cc : changed calls to sort to std::sort
  • ATrous/ReconSearch.cc : minor change from setCubeStats testing
  • Cubes/cubes.cc : Some tweaking to improve performance, and changing calls to sort to std::sort.
  • Detection/sorting.cc : Define a new class "Pair" to sort pairs of arrays in the same sense -- removes need for second sort function in sort.cc. Make use of stable_sort to preserve order of matching elements.
File size: 3.1 KB
Line 
1#include <iostream>
2#include <algorithm>
3#include <Utils/utils.hh>
4const int nsample=1000;
5// const int width=300;
6const int width=150;
7const float contrast=0.25;
8
9void zscale(long imagesize, float *image, float &z1, float &z2)
10{
11
12  float *smallarray = new float[nsample];
13  float *ct = new float[nsample];
14  long size=0;
15
16  float mean=0.,sd=0.;
17  for(int i=0;i<imagesize;i++) mean+=image[i];
18  mean /= float(imagesize);
19  for(int i=0;i<imagesize;i++) sd+=(image[i]-mean)*(image[i]-mean);
20  sd /= float(imagesize);
21 
22
23  if(imagesize>nsample){
24    int step = (imagesize / nsample) + 1;
25    for(int i=0;i<imagesize;i++){
26      if((i%step)==0){
27        smallarray[size] = image[i];
28        ct[size]=(float)size+1.;
29        size++;
30      }
31    }
32  }
33  else{
34    for(int i=0;i<imagesize;i++){
35      smallarray[i] = image[i];
36      ct[i] = float(i+1);
37    }
38    size=imagesize;
39  }
40
41  std::sort(smallarray,smallarray+size);
42 
43  /* fit a linear slope to the centre of the cumulative distribution */
44  long midpt = size/2;
45  long imin = midpt - width;
46  float slope,intercept,errSlope,errIntercept,r;
47  if(size<2*width)
48    linear_regression(size,ct,smallarray,0,size-1,slope,errSlope,
49                      intercept,errIntercept,r);
50  else linear_regression(size,ct,smallarray,imin,imin+2*width,slope,errSlope,
51                         intercept,errIntercept,r);
52
53  z1 = smallarray[midpt] + (slope/contrast)*(float)(1-midpt);
54  z2 = smallarray[midpt] + (slope/contrast)*(float)(nsample-midpt);
55
56  if(z1==z2){
57 
58    if(z1!=0) z2 = z1 * 1.05;
59    else z2 = z1+0.01;
60
61  }
62
63  delete [] smallarray;
64  delete [] ct;
65           
66}
67
68
69void zscale(long imagesize, float *image, float &z1, float &z2, float blankVal)
70{
71  float *newimage = new float[imagesize];
72  int newsize=0;
73  for(int i=0;i<imagesize;i++) if(image[i]!=blankVal) newimage[newsize++] = image[i];
74
75  float *smallarray = new float[nsample];
76  float *ct = new float[nsample];
77  long size=0;
78
79  float mean=0.,sd=0.;
80  for(int i=0;i<newsize;i++) mean+=newimage[i];
81  mean /= float(newsize);
82  for(int i=0;i<newsize;i++) sd+=(newimage[i]-mean)*(newimage[i]-mean);
83  sd /= float(newsize);
84 
85
86  if(newsize>nsample){
87    int step = (newsize / nsample) + 1;
88    for(int i=0;i<newsize;i++){
89      if((i%step)==0){
90        smallarray[size] = newimage[i];
91        ct[size]=(float)size+1.;
92        size++;
93      }
94    }
95  }
96  else{
97    for(int i=0;i<newsize;i++){
98      smallarray[i] = newimage[i];
99      ct[i] = float(i+1);
100    }
101    size=newsize;
102  }
103
104  std::sort(smallarray,smallarray+size);
105 
106  /* fit a linear slope to the centre of the cumulative distribution */
107  long midpt = size/2;
108  long imin = midpt - width;
109  float slope,intercept,errSlope,errIntercept,r;
110  if(size<2*width)
111    linear_regression(size,ct,smallarray,0,size-1,slope,errSlope,
112                      intercept,errIntercept,r);
113  else linear_regression(size,ct,smallarray,imin,imin+2*width,slope,errSlope,
114                    intercept,errIntercept,r);
115
116  z1 = smallarray[midpt] + (slope/contrast)*(float)(1-midpt);
117  z2 = smallarray[midpt] + (slope/contrast)*(float)(nsample-midpt);
118
119  if(z1==z2){
120 
121    if(z1!=0) z2 = z1 * 1.05;
122    else z2 = z1+0.01;
123
124  }
125
126  delete [] newimage;
127  delete [] smallarray;
128  delete [] ct;
129
130           
131}
132
133
Note: See TracBrowser for help on using the repository browser.