source: trunk/src/GenericEdgeDetector.cpp@ 2896

Last change on this file since 2896 was 2639, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CAS-2825

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: not available

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Correct DEC correction on R.A.


File size: 12.3 KB
Line 
1//
2// C++ Implimentation: GenericEdgeDetector
3//
4// Description:
5//
6//
7// Author: Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2012
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <math.h>
13
14#include <casa/Arrays/Vector.h>
15#include <casa/Arrays/Matrix.h>
16#include <casa/Arrays/ArrayMath.h>
17#include <casa/Arrays/ArrayIO.h>
18#include <casa/Utilities/GenSort.h>
19
20#include <coordinates/Coordinates/DirectionCoordinate.h>
21
22#include "GenericEdgeDetector.h"
23
24using namespace std ;
25using namespace casa ;
26
27namespace asap {
28
29GenericEdgeDetector::GenericEdgeDetector()
30 : EdgeDetector()
31{}
32
33GenericEdgeDetector::~GenericEdgeDetector()
34{}
35
36Vector<uInt> GenericEdgeDetector::detect()
37{
38 os_.origin(LogOrigin( "GenericEdgeDetector", "detect", WHERE )) ;
39
40 initDetect() ;
41
42 topixel() ;
43 countup() ;
44 thresholding() ;
45 labeling() ;
46 trimming() ;
47 selection() ;
48 tuning() ;
49
50 os_ << LogIO::DEBUGGING
51 << "Detected " << off_.nelements() << " integrations as OFF" << LogIO::POST ;
52
53 return off_ ;
54}
55
56void GenericEdgeDetector::parseOption( const Record &option )
57{
58 os_.origin(LogOrigin( "GenericEdgeDetector", "parseOption", WHERE )) ;
59
60 String name = "fraction" ;
61 if ( option.isDefined( name ) ) {
62 if ( option.dataType( name ) == TpString ) {
63 // should be "xx%" format
64 String fstr = option.asString( name ) ;
65 fstr = fstr.substr(0,fstr.size()-1) ;
66 fraction_ = String::toFloat( fstr ) * 0.01 ;
67 }
68 else {
69 fraction_ = option.asFloat( name ) ;
70 }
71 }
72 else {
73 fraction_ = 0.1 ; // default is 10%
74 }
75
76 name = "width" ;
77 if ( option.isDefined( name ) ) {
78 width_ = option.asFloat( name ) ;
79 }
80 else {
81 width_ = 0.5 ; // default is half of median separation
82 }
83
84 name = "elongated" ;
85 if ( option.isDefined( name ) ) {
86 elongated_ = option.asBool( name ) ;
87 }
88 else {
89 elongated_ = False ; // default is two-dimensional processing
90 }
91
92 os_ << "OPTION SUMMARY: " << endl
93 << " fraction=" << fraction_ << endl
94 << " width=" << width_ << endl
95 << " elongated=" << (elongated_ ? "True" : "False") << LogIO::POST ;
96}
97
98void GenericEdgeDetector::topixel()
99{
100// os_.origin(LogOrigin( "GenericEdgeDetector", "topixel", WHERE )) ;
101
102 setup() ;
103 // using DirectionCoordinate
104 Matrix<Double> xform(2,2) ;
105 xform = 0.0 ;
106 xform.diagonal() = 1.0 ;
107 DirectionCoordinate coord( MDirection::J2000,
108 Projection( Projection::SIN ),
109 cenx_, ceny_,
110 dx_, dy_,
111 xform,
112 0.5*Double(nx_-1),
113 0.5*Double(ny_-1) ) ;
114 Double *pdir_p = new Double[dir_.nelements()] ;
115 pdir_.takeStorage( dir_.shape(), pdir_p, TAKE_OVER ) ;
116 uInt len = time_.nelements() ;
117 Bool b ;
118 Double *dir_p = dir_.getStorage( b ) ;
119 Double *wdir_p = dir_p ;
120 Vector<Double> world ;
121 Vector<Double> pixel ;
122 IPosition vshape( 1, 2 ) ;
123 for ( uInt i = 0 ; i < len ; i++ ) {
124 world.takeStorage( vshape, wdir_p, SHARE ) ;
125 pixel.takeStorage( vshape, pdir_p, SHARE ) ;
126 coord.toPixel( pixel, world ) ;
127 pdir_p += 2 ;
128 wdir_p += 2 ;
129 }
130 dir_.putStorage( dir_p, b ) ;
131}
132
133void GenericEdgeDetector::setup()
134{
135 os_.origin(LogOrigin( "GenericEdgeDetector", "setup", WHERE )) ;
136
137 Double xmax, xmin, ymax, ymin ;
138 minMax( xmin, xmax, dir_.row( 0 ) ) ;
139 minMax( ymin, ymax, dir_.row( 1 ) ) ;
140 Double wx = ( xmax - xmin ) * 1.1 ;
141 Double wy = ( ymax - ymin ) * 1.1 ;
142
143 cenx_ = 0.5 * ( xmin + xmax ) ;
144 ceny_ = 0.5 * ( ymin + ymax ) ;
145 Double decCorr = cos( ceny_ ) ;
146
147 uInt len = time_.nelements() ;
148 Matrix<Double> dd = dir_.copy() ;
149 for ( uInt i = len-1 ; i > 0 ; i-- ) {
150 //dd(0,i) = ( dd(0,i) - dd(0,i-1) ) * decCorr ;
151 dd(0,i) = ( dd(0,i) - dd(0,i-1) ) * cos( 0.5*(dd(1,i-1)+dd(1,i)) ) ;
152 dd(1,i) = dd(1,i) - dd(1,i-1) ;
153 }
154 Vector<Double> dr( len-1 ) ;
155 Bool b ;
156 const Double *dir_p = dd.getStorage( b ) ;
157 const Double *x_p = dir_p + 2 ;
158 const Double *y_p = dir_p + 3 ;
159 for ( uInt i = 0 ; i < len-1 ; i++ ) {
160 dr[i] = sqrt( (*x_p) * (*x_p) + (*y_p) * (*y_p) ) ;
161 x_p += 2 ;
162 y_p += 2 ;
163 }
164 dir_.freeStorage( dir_p, b ) ;
165 Double med = median( dr, False, True, True ) ;
166 dy_ = med * width_ ;
167 dx_ = dy_ / decCorr ;
168
169 nx_ = uInt( ceil( wx / dx_ ) ) ;
170 ny_ = uInt( ceil( wy / dy_ ) ) ;
171
172 pcenx_ = 0.5 * Double( nx_ - 1 ) ;
173 pceny_ = 0.5 * Double( ny_ - 1 ) ;
174
175 os_ << LogIO::DEBUGGING
176 << "rangex=(" << xmin << "," << xmax << ")" << endl
177 << "rangey=(" << ymin << "," << ymax << ")" << endl
178 << "median separation between pointings: " << med << endl
179 << "dx=" << dx_ << ", dy=" << dy_ << endl
180 << "wx=" << wx << ", wy=" << wy << endl
181 << "nx=" << nx_ << ", ny=" << ny_ << LogIO::POST ;
182}
183
184void GenericEdgeDetector::countup()
185{
186 os_.origin(LogOrigin( "GenericEdgeDetector", "countup", WHERE )) ;
187
188 uInt *a_p = new uInt[nx_*ny_] ;
189 apix_.takeStorage( IPosition(2,nx_,ny_), a_p, TAKE_OVER ) ;
190 apix_ = 0 ;
191
192 uInt len = time_.nelements() ;
193 uInt ix ;
194 uInt iy ;
195 // pdir_ is always contiguous
196 const Double *pdir_p = pdir_.data() ;
197 const Double *px_p = pdir_p ;
198 const Double *py_p = pdir_p + 1 ;
199 long offset ;
200 for ( uInt i = 0 ; i < len ; i++ ) {
201 ix = uInt(round( *px_p )) ;
202 iy = uInt(round( *py_p )) ;
203 offset = ix + iy * nx_ ;
204 *(a_p+offset) += 1 ;
205 px_p += 2 ;
206 py_p += 2 ;
207 }
208
209 os_ << LogIO::DEBUGGING
210 << "a.max()=" << max(apix_) << ",a.min()=" << min(apix_) << LogIO::POST ;
211}
212
213void GenericEdgeDetector::thresholding()
214{
215 uInt len = apix_.nelements() ;
216 uInt *a_p = apix_.data() ;
217 for ( uInt i = 0 ; i < len ; i++ ) {
218 *a_p = ((*a_p > 0) ? 1 : 0) ;
219 a_p++ ;
220 }
221}
222
223void GenericEdgeDetector::labeling()
224{
225 os_.origin(LogOrigin( "GenericEdgeDetector", "labeling", WHERE )) ;
226
227 uInt n = 1 ;
228 uInt niter = 0 ;
229 const uInt maxiter = 100 ;
230 while ( n > 0 && niter < maxiter ) {
231 n = _labeling() ;
232 os_ << LogIO::DEBUGGING
233 << "cycle " << niter << ": labeled " << n << " pixels" << LogIO::POST ;
234 niter++ ;
235 }
236 if ( niter == maxiter ) {
237 // WARN
238 os_ << LogIO::WARN << "labeling not converged before maxiter=" << maxiter << LogIO::POST ;
239 }
240}
241
242uInt GenericEdgeDetector::_labeling()
243{
244 uInt n = 0 ;
245 for ( uInt ix = 0 ; ix < nx_ ; ix++ ) {
246 Vector<uInt> v( apix_.row( ix ) ) ;
247 uInt nx = __labeling( v ) ;
248 n += nx ;
249 }
250 for ( uInt iy = 0 ; iy < ny_ ; iy++ ) {
251 Vector<uInt> v( apix_.column( iy ) ) ;
252 uInt ny = __labeling( v ) ;
253 n += ny ;
254 }
255 return n ;
256}
257
258uInt GenericEdgeDetector::__labeling( Vector<uInt> &a )
259{
260 uInt n = 0 ;
261 if ( allEQ( a, n ) ) {
262 return n ;
263 }
264
265 uInt start ;
266 uInt end ;
267 _search( start, end, a ) ;
268 for ( uInt i = start+1 ; i < end ; i++ ) {
269 if ( a[i] == 0 ) {
270 a[i] = 1 ;
271 n++ ;
272 }
273 }
274 return n ;
275}
276
277void GenericEdgeDetector::_search( uInt &start,
278 uInt &end,
279 const Vector<uInt> &a )
280{
281 uInt n = a.nelements() ;
282 start = 0 ;
283 while( a[start] == 0 ) {
284 start++ ;
285 }
286 end = n - 1 ;
287 while( a[end] == 0 ) {
288 end-- ;
289 }
290}
291
292void GenericEdgeDetector::trimming()
293{
294 os_.origin(LogOrigin( "GenericEdgeDetector", "trimming", WHERE )) ;
295
296 const uInt n1 = sum( apix_ ) ;
297 const uInt nTrim = uInt(ceil( n1 * fraction_ )) ;
298 os_ << LogIO::DEBUGGING
299 << "number of nonzero pixel: " << n1 << endl
300 << "number of pixels to be trimmed: " << nTrim << LogIO::POST ;
301 uInt n = 0 ;
302 uInt niter = 0 ;
303 const uInt maxiter = 100 ;
304 if ( !elongated_ ) {
305 while ( n < nTrim && niter < maxiter ) {
306 uInt m = _trimming() ;
307 os_ << LogIO::DEBUGGING
308 << "cycle " << niter << ": trimmed " << m << " pixels" << LogIO::POST ;
309 n += m ;
310 niter++ ;
311 }
312 }
313 else if ( nx_ > ny_ ) {
314 os_ << "1D triming along x-axis" << LogIO::POST ;
315 while ( n < nTrim && niter < maxiter ) {
316 uInt m = _trimming1DX() ;
317 os_ << LogIO::DEBUGGING
318 << "cycle " << niter << ": trimmed " << m << " pixels" << LogIO::POST ;
319 n += m ;
320 niter++ ;
321 }
322 }
323 else { // nx_ < ny_
324 os_ << "1D triming along y-axis" << LogIO::POST ;
325 while ( n < nTrim && niter < maxiter ) {
326 uInt m = _trimming1DY() ;
327 os_ << LogIO::DEBUGGING
328 << "cycle " << niter << ": trimmed " << m << " pixels" << LogIO::POST ;
329 n += m ;
330 niter++ ;
331 }
332 }
333 os_ << LogIO::DEBUGGING
334 << "number of pixels actually trimmed: " << n << LogIO::POST ;
335
336 if ( niter == maxiter ) {
337 // WARN
338 os_ << LogIO::WARN << "trimming not converged before maxiter=" << maxiter << LogIO::POST ;
339 }
340}
341
342uInt GenericEdgeDetector::_trimming()
343{
344 uInt n = 0 ;
345 Block<uInt> flatIdxList( apix_.nelements() ) ;
346 uInt start ;
347 uInt end ;
348 uInt flatIdx ;
349 for ( uInt ix = 0 ; ix < nx_ ; ix++ ) {
350 Vector<uInt> v( apix_.row( ix ) ) ;
351 if ( allEQ( v, (uInt)0 ) ) {
352 continue ;
353 }
354 _search( start, end, v ) ;
355 uInt offset = start * nx_ ;
356 flatIdx = offset + ix ;
357 flatIdxList[n++] = flatIdx ;
358 if ( start != end ) {
359 offset = end * nx_ ;
360 flatIdx = offset + ix ;
361 flatIdxList[n++] = flatIdx ;
362 }
363 }
364 for ( uInt iy = 0 ; iy < ny_ ; iy++ ) {
365 Vector<uInt> v( apix_.column( iy ) ) ;
366 if ( allEQ( v, (uInt)0 ) ) {
367 continue ;
368 }
369 uInt offset = iy * nx_ ;
370 _search( start, end, v ) ;
371 flatIdx = offset + start ;
372 flatIdxList[n++] = flatIdx ;
373 if ( start != end ) {
374 flatIdx = offset + end ;
375 flatIdxList[n++] = flatIdx ;
376 }
377 }
378 n = genSort( flatIdxList.storage(),
379 n,
380 Sort::Ascending,
381 Sort::QuickSort | Sort::NoDuplicates ) ;
382 Vector<uInt> v( IPosition(1,apix_.nelements()), apix_.data(), SHARE ) ;
383 const uInt *idx_p = flatIdxList.storage() ;
384 for ( uInt i = 0 ; i < n ; i++ ) {
385 v[*idx_p] = 0 ;
386 idx_p++ ;
387 }
388
389 return n ;
390}
391
392uInt GenericEdgeDetector::_trimming1DX()
393{
394 uInt n = 0 ;
395 const uInt nx = apix_.nrow() ;
396 Vector<uInt> v1, v2 ;
397 uInt ix, jx ;
398 for ( ix = 0 ; ix < nx_ ; ix++ ) {
399 v1.reference( apix_.row( ix ) ) ;
400 if ( anyNE( v1, n ) ) break ;
401 }
402 for ( jx = nx-1 ; jx > ix ; jx-- ) {
403 v2.reference( apix_.row( jx ) ) ;
404 if ( anyNE( v2, n ) ) break ;
405 }
406 n += _trimming1D( v1 ) ;
407 if ( ix != jx )
408 n+= _trimming1D( v2 ) ;
409
410 return n ;
411}
412
413uInt GenericEdgeDetector::_trimming1DY()
414{
415 uInt n = 0 ;
416 const uInt ny = apix_.ncolumn() ;
417 Vector<uInt> v1, v2 ;
418 uInt iy, jy ;
419 for ( iy = 0 ; iy < ny_ ; iy++ ) {
420 v1.reference( apix_.column( iy ) ) ;
421 if ( anyNE( v1, n ) ) break ;
422 }
423 for ( jy = ny-1 ; jy > iy ; jy-- ) {
424 v2.reference( apix_.column( jy ) ) ;
425 if ( anyNE( v2, n ) ) break ;
426 }
427 n += _trimming1D( v1 ) ;
428 if ( iy != jy )
429 n+= _trimming1D( v2 ) ;
430
431 return n ;
432}
433
434uInt GenericEdgeDetector::_trimming1D( Vector<uInt> &a )
435{
436 uInt len = a.nelements() ;
437 uInt n = 0 ;
438 for ( uInt i = 0 ; i < len ; i++ ) {
439 if ( a[i] == 1 ) {
440 a[i] = 0 ;
441 n++ ;
442 }
443 }
444
445 return n ;
446}
447
448void GenericEdgeDetector::selection()
449{
450// os_.origin(LogOrigin( "GenericEdgeDetector", "selection", WHERE )) ;
451
452 uInt nrow = pdir_.shape()[1] ;
453 const Double *px_p = pdir_.data() ;
454 const Double *py_p = px_p + 1 ;
455 Vector<uInt> v( IPosition(1,apix_.nelements()), apix_.data(), SHARE ) ;
456 uInt n = 0 ;
457 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
458 uInt idx = int(round(*px_p)) + int(round(*py_p)) * nx_ ;
459 if ( v[idx] == 0 ) {
460 tempuInt_[n++] = irow ;
461 }
462 px_p += 2 ;
463 py_p += 2 ;
464 }
465 off_ = vectorFromTempStorage( n ) ;
466}
467
468void GenericEdgeDetector::tuning()
469{
470 os_.origin(LogOrigin( "GenericEdgeDetector", "tuning", WHERE )) ;
471
472 const uInt len = off_.nelements() ;
473 if ( len == 0 )
474 return ;
475
476 Block<uInt> diff( len-1 ) ;
477 for ( uInt i = 0 ; i < len-1 ; i++ ) {
478 diff[i] = off_[i+1] - off_[i] ;
479 }
480 const uInt threshold = 3 ;
481 uInt n = 0 ;
482 for ( uInt i = 0 ; i < len ; i++ ) {
483 tempuInt_[n++] = off_[i] ;
484 }
485 for ( uInt i = 1 ; i < len ; i++ ) {
486 uInt ii = i - 1 ;
487 if ( diff[ii] != 1 && diff[ii] < threshold ) {
488 uInt t = off_[ii]+1 ;
489 uInt u = off_[i] ;
490 for ( uInt j = t ; j < u ; j++ ) {
491 os_ << LogIO::DEBUGGING
492 << "move " << j << " from ON to OFF" << LogIO::POST ;
493 tempuInt_[n++] = j ;
494 }
495 }
496 }
497 if ( n > len ) {
498 off_.resize() ;
499 off_ = vectorFromTempStorage( n ) ;
500 }
501}
502
503} // namespace asap
Note: See TracBrowser for help on using the repository browser.