source: tags/release-1.0.7/src/Cubes/trimImage.cc

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

A suite of changes aimed at calculating a cube-wide threshold, rather than individual thresholds for each channel map & spectrum. Summary of changes:

  • New Cube::setCubeStats() function, that stores the cube's stats in Cube member variables (mean, median, stddev, madfm).
  • New Parameter threshold, which can be entered by the user.
  • Make use of the cube statistics in the search functions, so that the Image class receives the correct information. Some tidying up of the CubicSearch? and ReconSearch? functions as well.
  • Changing the thresholding functions to be able to use the new threshold parameter.

Warning though -- the FDR setup doesn't quite work at this stage!

File size: 7.3 KB
Line 
1#include <iostream>
2#include <param.hh>
3#include <Cubes/cubes.hh>
4
5void Cube::trimCube()
6{
7  /**
8   *  Cube::trimCube()
9   *
10   *   If the blankPix flag has been set, this routine trims excess blank
11   *    pixels from the edges of the spatial image.
12   *   It uses as its template the first channel, assuming that its non-BLANK
13   *    size is representative of the rest of the channels.
14   *   The edges are moved in until the first non-BLANK pixel is encountered.
15   *   All other arrays are similarly edited, and the amount of trimming is
16   *    recorded.
17   */
18
19  if(this->par.getFlagBlankPix()) {
20
21    long xdim = this->axisDim[0];
22    long ydim = this->axisDim[1];
23    long zdim = this->axisDim[2];
24    const long ZCHAN = 0;
25
26    for(int z = 0; z < zdim; z++){
27
28      for(int row=0;row<ydim;row++){
29        int borderLeft=0;
30        int borderRight=0;
31        while((borderLeft<xdim)&&
32              (this->par.isBlank(this->array[row*xdim+borderLeft+ZCHAN*xdim*ydim])) ) {
33          borderLeft++;
34        }
35        while((borderRight<xdim)&&
36              (this->par.isBlank(this->array[row*xdim+(xdim-1-borderRight)+ZCHAN*xdim*ydim]))){
37          borderRight++;
38        }
39
40        if( ((z==0)&&(row==0)) || (borderLeft < this->par.getBorderLeft()) )
41          this->par.setBorderLeft(borderLeft);
42        if( ((z==0)&&(row==0)) || (borderRight < this->par.getBorderRight()) )
43          this->par.setBorderRight(borderRight);
44
45      }
46   
47      for(int col=0;col<xdim;col++){
48        int borderBottom=0;
49        int borderTop=0;
50        while((borderBottom<ydim)&&
51              (this->par.isBlank(this->array[col+borderBottom*xdim+ZCHAN*xdim*ydim])) ) borderBottom++;
52
53        while((borderTop<ydim)&&
54              (this->par.isBlank(this->array[col+(ydim-1-borderTop)*xdim+ZCHAN*xdim*ydim])) ) borderTop++;
55
56        if( ((z==0)&&(col==0)) || (borderBottom < this->par.getBorderBottom()) )
57          this->par.setBorderBottom(borderBottom);
58        if( ((z==0)&&(col==0)) || (borderTop < this->par.getBorderTop()) )
59          this->par.setBorderTop(borderTop);
60
61      }
62
63    }
64   
65    this->axisDim[0] = xdim - this->par.getBorderLeft() - this->par.getBorderRight();
66    this->axisDim[1] = ydim - this->par.getBorderBottom() - this->par.getBorderTop();
67    this->axisDim[2] = zdim;
68    this->numPixels = this->axisDim[0]*this->axisDim[1]*this->axisDim[2];
69
70    long oldpos,newpos;
71   
72    // Do the trimming, but only if we need to -- is there a border of Blanks?
73    if((this->par.getBorderLeft()>0)||(this->par.getBorderRight()>0)||
74       (this->par.getBorderBottom()>0)||(this->par.getBorderTop()>0)) {
75
76      // Trim the array of pixel values
77      float *newarray  = new float[this->numPixels];
78      for(int x = 0; x < axisDim[0]; x++){
79        for(int y = 0; y < axisDim[1]; y++){
80          for(int z = 0; z < axisDim[2]; z++){
81            oldpos = (x+this->par.getBorderLeft()) + (y+this->par.getBorderBottom())*xdim + z*xdim*ydim;
82            newpos = x + y*axisDim[0] + z*axisDim[0]*axisDim[1];
83            newarray[newpos]  = this->array[oldpos];
84          }
85        }
86      }
87      delete [] this->array;
88      this->array = newarray;
89
90      // Trim the 2-D detection map
91      short *newdetect = new short[this->axisDim[0]*this->axisDim[1]];
92      for(int x = 0; x < axisDim[0]; x++){
93        for(int y = 0; y < axisDim[1]; y++){
94          oldpos = (x+this->par.getBorderLeft()) + (y+this->par.getBorderBottom())*xdim;
95          newpos = x + y*axisDim[0];
96          newdetect[newpos] = this->detectMap[oldpos];
97        }
98      }
99      delete [] this->detectMap;
100      this->detectMap = newdetect;
101
102      if(this->par.getFlagATrous()){
103        // Trim the reconstructed array if we are going to do the reconstruction
104        float *newrecon  = new float[this->numPixels];
105        for(int x = 0; x < axisDim[0]; x++){
106          for(int y = 0; y < axisDim[1]; y++){
107            for(int z = 0; z < axisDim[2]; z++){
108              oldpos = (x+this->par.getBorderLeft()) + (y+this->par.getBorderBottom())*xdim + z*xdim*ydim;
109              newpos = x + y*axisDim[0] + z*axisDim[0]*axisDim[1];
110              newrecon[newpos] = this->recon[oldpos];     
111            }
112          }
113        }
114        delete [] this->recon;
115        this->recon = newrecon;
116      }
117
118      // Set the flag indicating trimming has taken place only if it has.
119      this->par.setFlagCubeTrimmed(true);
120
121    }
122
123  }
124
125}
126
127
128void Cube::unTrimCube()
129{
130  /**
131   *  Cube::unTrimCube()
132   *
133   *   If the cube has been trimmed by trimCube(), this task adds back the BLANK pixels on
134   *    the edges, so that it returns to its original size.
135   *   All arrays are similarly edited.
136   */
137
138  if(this->par.getFlagCubeTrimmed()){
139
140    long smallXDim = this->axisDim[0];
141    long smallYDim = this->axisDim[1];
142    long smallZDim = this->axisDim[2];
143
144    // first correct the dimension sizes
145    this->axisDim[0] = smallXDim + this->par.getBorderLeft() + this->par.getBorderRight();
146    this->axisDim[1] = smallYDim + this->par.getBorderBottom() + this->par.getBorderTop();
147    this->axisDim[2] = smallZDim;
148    this->numPixels = this->axisDim[0]*this->axisDim[1]*this->axisDim[2];
149
150    long pos,smlpos;
151    bool isDud;
152
153    // Correct the array of pixel values
154    float *newarray  = new float[this->numPixels];
155    for(int x = 0; x < this->axisDim[0]; x++){
156      for(int y = 0; y < this->axisDim[1]; y++){
157        isDud = (x<this->par.getBorderLeft()) || (x>=smallXDim+this->par.getBorderLeft()) ||
158          (y<this->par.getBorderBottom()) || (y>=smallYDim+this->par.getBorderBottom());
159       
160        for(int z = 0; z < this->axisDim[2]; z++){
161          pos = x + y*this->axisDim[0] + z*this->axisDim[0]*this->axisDim[1];
162          smlpos = (x-this->par.getBorderLeft()) + (y-this->par.getBorderBottom())*smallXDim
163            + z * smallXDim * smallYDim;
164          if(isDud) newarray[pos] = this->par.getBlankPixVal();
165          else      newarray[pos] = this->array[smlpos];
166        }
167      }
168    }
169    delete [] this->array;
170    this->array = newarray;
171
172    if(this->reconExists){
173      // Correct the reconstructed array
174      float *newrecon   = new float[this->numPixels];
175      for(int x = 0; x < this->axisDim[0]; x++){
176        for(int y = 0; y < this->axisDim[1]; y++){
177          isDud = (x<this->par.getBorderLeft()) || (x>=smallXDim+this->par.getBorderLeft()) ||
178            (y<this->par.getBorderBottom()) || (y>=smallYDim+this->par.getBorderBottom());
179         
180          for(int z = 0; z < this->axisDim[2]; z++){
181            pos = x + y*this->axisDim[0] + z*this->axisDim[0]*this->axisDim[1];
182            smlpos = (x-this->par.getBorderLeft()) + (y-this->par.getBorderBottom())*smallXDim
183              + z * smallXDim * smallYDim;
184            if(isDud) newrecon[pos] = this->par.getBlankPixVal();
185            else      newrecon[pos] = this->recon[smlpos];
186          }
187        }
188      }
189      delete [] this->recon;
190      this->recon = newrecon;
191    }
192
193    // Correct the 2-D detection map
194    short *newdetect = new short[this->axisDim[0]*this->axisDim[1]];
195    for(int x = 0; x < this->axisDim[0]; x++){
196      for(int y = 0; y < this->axisDim[1]; y++){
197        pos = x + y*this->axisDim[0];
198        smlpos = (x-this->par.getBorderLeft()) + (y-this->par.getBorderBottom())*smallXDim;
199        isDud = (x<this->par.getBorderLeft()) || (x>=smallXDim+this->par.getBorderLeft()) ||
200          (y<this->par.getBorderBottom()) || (y>=smallYDim+this->par.getBorderBottom());
201        if(isDud) newdetect[pos]=0;
202        else newdetect[pos] = this->detectMap[smlpos];
203      }
204    }
205    delete [] this->detectMap;
206    this->detectMap = newdetect;   
207
208 
209    // Now update the positions for all the detections
210 
211    for(int i=0;i<this->objectList.size();i++){
212      for(int pix=0;pix<objectList[i].getSize();pix++){
213        long x = objectList[i].getX(pix);
214        long y = objectList[i].getY(pix);
215        objectList[i].setX(pix,x+this->par.getBorderLeft());
216        objectList[i].setY(pix,y+this->par.getBorderBottom());
217      }
218      objectList[i].calcParams();
219    }
220
221  }
222
223}
Note: See TracBrowser for help on using the repository browser.