source: tags/release-1.0.7/src/Detection/lutz_detect.cc @ 1441

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

Large commit. The addition of a new Statistics namespace & class, plus a number of changes to the code to make the cube-wide statistics calculation work. The FDR calculations are incorporated into the new class, and a number of functions have been made into templates to ease the calculations. Details follow:

  • New namespace and class (StatsContainer? -- templated) in Statistics.hh, that holds mean,median,stddev & madfm, and provide accessor and calculator functions for these. It also holds the threshold values for sigma-clipping and FDR methods, as well as the PValue evaluation functions
  • The correctionFactor incorporated into the namespace, and given a conversion function that other functions can call (eg. the atrous_Xd_reconstruct functions).
  • Templated the statistics functions in getStats.cc.
  • Templated the sort functions, and made swap an inline one defined in utils.hh.
  • A number of changes to cubes.cc and cubes.hh:
    • DataArray? gains a StatsContainer? object, to store stats info.
    • Image has lost its pValue array (not needed as we calculate on the fly) and mask array (not used).
    • Image has also lost all its stats members, but keeps minPix.
    • Functions to go are Image::maskObject, Image::findStats. Removed calls to the former. Latter never used.
    • Cube::setCubeStats does the cube-wide stats calculations, including setupFDR (now a Cube function).
    • Cube loses the specMean etc arrays.
  • The Search functions (ReconSearch? and CubicSearch?) changed to accommodate the exchange of StatsContainer? objects. This changed the prototypes as well.
  • The growObject function incorporates the new StatsContainer? object.
  • Minor change to Merger, in the preparation for calling growObject.
  • A new par introduced: flagUserThreshold -- set to true when the user enters a value for the threshold.
  • Removed thresholding_functions.cc and incorporated its functions into cubes.cc and cubes.hh.
File size: 6.5 KB
Line 
1/**
2 *  LUTZ_DETECT.CC
3 *
4 * A detection algorithm based on that of Lutz (1980).
5 *
6 * INPUTS:
7 *    image     -- an Image object, containing a 2D image that has been
8 *                 processed such that its pValue array is defined.
9 * OUTPUTS:
10 *   The detection array in image will be filled, according to
11 *   the location of the objects in the image.
12 *
13 */
14#include <Cubes/cubes.hh>
15#include <Detection/detection.hh>
16#include <vector>
17
18enum STATUS { NONOBJECT, OBJECT, COMPLETE, INCOMPLETE };
19enum ROW { PRIOR=0, CURRENT};
20enum NULLS { NULLSTART=-1, NULLMARKER=45}; //ASCII 45 = '-' --> eases printing in case of debugging
21
22class Object{
23public:
24  Object(){start=NULLSTART; end=NULLSTART;};
25  int start;
26  int end;
27  Detection info;
28};
29
30void Image::lutz_detect()
31{
32  // Allocate necessary arrays.
33  STATUS *status   = new STATUS[2];
34  Detection *store = new Detection[this->axisDim[0]+1];
35  char *marker     = new char[this->axisDim[0]+1];
36  for(int i=0; i<(this->axisDim[0]+1); i++) marker[i] = NULLMARKER;
37  vector <Object>  *oS    = new vector <Object>;
38  vector <STATUS> *psS    = new vector <STATUS>;
39
40  Pixel *pix = new Pixel;
41
42  for(int posY=0;posY<(this->axisDim[1]+1);posY++){
43    // Loop over each row -- consider rows one at a time
44   
45    status[PRIOR] = COMPLETE;
46    status[CURRENT] = NONOBJECT;
47
48    pix->setY(posY);
49
50    for(int posX=0;posX<(this->axisDim[0]+1);posX++){
51      // Now the loop for a given row, looking at each column individually.
52
53      char currentMarker = marker[posX];
54      marker[posX] = NULLMARKER;
55      pix->setX(posX);
56
57      bool isObject = false;
58      if((posX<this->axisDim[0])&&(posY<this->axisDim[1])){
59        // if we are in the original image
60        isObject = this->isDetection(posX,posY);
61        pix->setF( this->array[posY*this->axisDim[0] + posX] );
62      }
63      // else we're in the padding row/col and isObject=FALSE;
64
65      /**
66       * ------------------------------ START SEGMENT ------------------------
67       * If the current pixel is object and the previous pixel is not, then
68       * start a new segment.
69       * If the pixel touches an object on the prior row, the marker is either
70       * an S or an s, depending on whether the object has started yet.
71       * If it doesn't touch a prior object, this is the start of a completly
72       * new object on this row.
73       */
74      if ( (isObject) && (status[CURRENT] != OBJECT) ){
75
76        status[CURRENT] = OBJECT;
77        if(status[PRIOR] == OBJECT){
78         
79          if(oS->back().start==NULLSTART){
80            marker[posX] = 'S';
81            oS->back().start = posX;
82          }
83          else  marker[posX] = 's';
84        }
85        else{
86          psS->push_back(status[PRIOR]);  //PUSH PS onto PSSTACK;
87          marker[posX] = 'S';
88          oS->resize(oS->size()+1);        //PUSH OBSTACK;
89          oS->back().start = posX;
90
91          status[PRIOR] = COMPLETE;
92        }
93      }
94
95      /**
96       * ------------------------------ PROCESS MARKER -----------------------
97       * If the current marker is not blank, then we need to deal with it.
98       * Four cases:
99       *   S --> start of object on prior row. Push priorStatus onto PSSTACK
100       *         and set priorStatus to OBJECT
101       *   s --> start of a secondary segment of object on prior row.
102       *         If current object joined, pop PSSTACK and join the objects.
103       *         Set priorStatus to OBJECT.
104       *   f --> end of a secondary segment of object on prior row.
105       *         Set priorStatus to INCOMPLETE.
106       *   F --> end of object on prior row. If no more of the object is to
107       *          come (priorStatus=COMPLETE), then finish it and output data.
108       *         Add to list, but only if it has more than the minimum number
109       *           of pixels.
110       */
111      if(currentMarker != NULLMARKER){
112
113        if(currentMarker == 'S'){
114          psS->push_back(status[PRIOR]);      // PUSH PS onto PSSTACK
115          if(status[CURRENT] == NONOBJECT){
116            psS->push_back(COMPLETE);         // PUSH COMPLETE ONTO PSSTACK;
117            oS->resize(oS->size()+1);          // PUSH OBSTACK;
118            oS->back().info = store[posX];
119          }
120          else oS->back().info.addAnObject( store[posX] );
121         
122          status[PRIOR] = OBJECT;
123        }
124
125        /*---------*/
126        if(currentMarker == 's'){
127
128          if( (status[CURRENT] == OBJECT) && (status[PRIOR] == COMPLETE) ){
129            status[PRIOR] = psS->back();
130            psS->pop_back();                   //POP PSSTACK ONTO PS
131
132//          oS->at(oS->size()-2).info.addAnObject( oS->back().info );
133//          if(oS->at(oS->size()-2).start == NULLSTART) oS->at(oS->size()-2).start = oS->back().start;
134//          else marker[oS->back().start] = 's';
135            (*oS)[oS->size()-2].info.addAnObject( oS->back().info );
136            if((*oS)[oS->size()-2].start == NULLSTART)
137              (*oS)[oS->size()-2].start = oS->back().start;
138            else marker[oS->back().start] = 's';
139
140            oS->pop_back();
141          }
142
143          status[PRIOR] = OBJECT;
144        }
145
146        /*---------*/
147        if(currentMarker == 'f') status[PRIOR] = INCOMPLETE;
148
149        /*---------*/
150        if(currentMarker == 'F') {
151
152          status[PRIOR] = psS->back();
153          psS->pop_back();                    //POP PSSTACK ONTO PS
154
155          if( (status[CURRENT] == NONOBJECT) && (status[PRIOR] == COMPLETE) ){
156
157            if(oS->back().start == NULLSTART){ // object completed
158              if(oS->back().info.getSize() >= this->minSize){
159                // is it big enough?
160                oS->back().info.calcParams(); // work out midpoints, fluxes etc
161                this->addObject(oS->back().info);
162              }
163            }
164            else{
165              marker[ oS->back().end ] = 'F';
166              store[ oS->back().start ] = oS->back().info;
167            }
168
169            oS->pop_back();
170
171            status[PRIOR] = psS->back();
172            psS->pop_back();
173          }
174        }
175
176      } // end of PROCESSMARKER section ( if(currentMarker!=NULLMARKER) )
177
178      if (isObject){
179        oS->back().info.addPixel(*pix);
180      }
181      else{
182        /**
183         *------------------------------ END SEGMENT -------------------------
184         * If the current pixel is background and the previous pixel was an
185         * object, then finish the segment.
186         * If the prior status is COMPLETE, it's the end of the final segment
187         * of the object section.
188         * If not, it's end of the segment, but not necessarily the section.
189         */
190        if ( status[CURRENT] == OBJECT) {
191
192          status[CURRENT] = NONOBJECT;
193
194          if(status[PRIOR] != COMPLETE){
195            marker[posX] = 'f';
196            oS->back().end = posX;
197          }
198          else{
199            status[PRIOR] = psS->back();
200            psS->pop_back();                   //POP PSSTACK onto PS;
201            marker[posX] = 'F';
202            store[ oS->back().start ] = oS->back().info;
203            oS->pop_back();
204          }
205        }
206       
207      } //-> end of END SEGMENT else{ clause
208
209    }//-> end of loop over posX
210
211  }//-> end of loop over posY
212
213  // clean up and remove declared arrays
214  delete [] marker;
215  delete oS;
216  delete psS;
217  delete [] store;
218  delete [] status;
219  delete pix;
220
221}
Note: See TracBrowser for help on using the repository browser.