source: trunk/src/STGrid.cpp@ 2377

Last change on this file since 2377 was 2377, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: Yes

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Disabled table sort.


File size: 27.3 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include <casa/BasicSL/String.h>
5#include <casa/Arrays/Vector.h>
6#include <casa/Arrays/Matrix.h>
7#include <casa/Arrays/Cube.h>
8#include <casa/Arrays/ArrayMath.h>
9#include <casa/Arrays/ArrayIter.h>
10#include <casa/Quanta/Quantum.h>
11#include <casa/Quanta/QuantumHolder.h>
12#include <casa/Utilities/CountedPtr.h>
13#include <casa/Logging/LogIO.h>
14
15#include <tables/Tables/Table.h>
16#include <tables/Tables/TableRecord.h>
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/ScalarColumn.h>
19#include <tables/Tables/ArrayColumn.h>
20
21#include <measures/Measures/MDirection.h>
22
23#include <MathUtils.h>
24
25#include "STGrid.h"
26
27using namespace std ;
28using namespace casa ;
29using namespace asap ;
30
31namespace asap {
32
33// constructor
34STGrid::STGrid()
35{
36 init() ;
37}
38
39STGrid::STGrid( const string infile )
40{
41 init() ;
42
43 setFileIn( infile ) ;
44}
45
46void STGrid::init()
47{
48 ifno_ = -1 ;
49 nx_ = -1 ;
50 ny_ = -1 ;
51 npol_ = 0 ;
52 nchan_ = 0 ;
53 nrow_ = 0 ;
54 cellx_ = 0.0 ;
55 celly_ = 0.0 ;
56 center_ = Vector<Double> ( 2, 0.0 ) ;
57 convType_ = "BOX" ;
58 wtype_ = "UNIFORM" ;
59 convSupport_ = -1 ;
60 userSupport_ = -1 ;
61 convSampling_ = 100 ;
62}
63
64void STGrid::setFileIn( const string infile )
65{
66 String name( infile ) ;
67 if ( infile_.compare( name ) != 0 ) {
68 infile_ = String( infile ) ;
69 tab_ = Table( infile_ ) ;
70 }
71}
72
73void STGrid::setPolList( vector<unsigned int> pols )
74{
75 pollist_.assign( Vector<uInt>( pols ) ) ;
76 cout << "pollist_ = " << pollist_ << endl ;
77}
78
79void STGrid::setScanList( vector<unsigned int> scans )
80{
81 scanlist_.assign( Vector<uInt>( scans ) ) ;
82 cout << "scanlist_ = " << scanlist_ << endl ;
83}
84
85void STGrid::setWeight( const string wType )
86{
87 wtype_ = String( wType ) ;
88 wtype_.upcase() ;
89 cout << "wtype_ = " << wtype_ << endl ;
90}
91
92void STGrid::defineImage( int nx,
93 int ny,
94 string scellx,
95 string scelly,
96 string scenter )
97{
98 ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
99 Matrix<Double> direction = dirCol.getColumn() ;
100 Double rmax, rmin, dmax, dmin ;
101 minMax( rmin, rmax, direction.row( 0 ) ) ;
102 minMax( dmin, dmax, direction.row( 1 ) ) ;
103
104 Int npx = (Int)nx ;
105 Int npy = (Int)ny ;
106 String cellx( scellx ) ;
107 String celly( scelly ) ;
108 String center( scenter ) ;
109 setupGrid( npx, npy,
110 cellx, celly,
111 rmin, rmax,
112 dmin, dmax,
113 center ) ;
114}
115
116void STGrid::setFunc( string convType,
117 int convSupport )
118{
119 convType_ = String( convType ) ;
120 convType_.upcase() ;
121 userSupport_ = (Int)convSupport ;
122}
123
124#define NEED_UNDERSCORES
125#if defined(NEED_UNDERSCORES)
126#define ggridsd ggridsd_
127#endif
128extern "C" {
129 void ggridsd(Double*,
130 const Complex*,
131 Int*,
132 Int*,
133 Int*,
134 const Int*,
135 const Int*,
136 const Float*,
137 Int*,
138 Int*,
139 Complex*,
140 Float*,
141 Int*,
142 Int*,
143 Int *,
144 Int *,
145 Int*,
146 Int*,
147 Float*,
148 Int*,
149 Int*,
150 Double*);
151}
152void STGrid::grid()
153{
154 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
155
156 // retrieve data
157 Array<Complex> spectra ;
158 Array<Double> direction ;
159 Array<Int> flagtra ;
160 Array<Int> rflag ;
161 Array<Float> weight ;
162 double t0, t1 ;
163 t0 = mathutil::gettimeofday_sec() ;
164 getData( spectra, direction, flagtra, rflag, weight ) ;
165 t1 = mathutil::gettimeofday_sec() ;
166 os << "getData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
167 IPosition sshape = spectra.shape() ;
168 //os << "spectra.shape()=" << spectra.shape() << LogIO::POST ;
169 //os << "max(spectra) = " << max(spectra) << LogIO::POST ;
170 //os << "weight = " << weight << LogIO::POST ;
171
172 // grid parameter
173 os << LogIO::DEBUGGING ;
174 os << "----------" << endl ;
175 os << "Grid parameter summary" << endl ;
176 os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
177 os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
178 os << " center = " << center_ << endl ;
179 os << "----------" << LogIO::POST ;
180 os << LogIO::NORMAL ;
181
182 // convolution kernel
183 Vector<Float> convFunc ;
184 t0 = mathutil::gettimeofday_sec() ;
185 setConvFunc( convFunc ) ;
186 t1 = mathutil::gettimeofday_sec() ;
187 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
188 //cout << "convSupport=" << convSupport_ << endl ;
189 //cout << "convFunc=" << convFunc << endl ;
190
191 // world -> pixel
192 Array<Double> xypos( direction.shape(), 0.0 ) ;
193 t0 = mathutil::gettimeofday_sec() ;
194 toPixel( direction, xypos ) ;
195 t1 = mathutil::gettimeofday_sec() ;
196 //os << "xypos=" << xypos << LogIO::POST ;
197 os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
198
199 // call ggridsd
200 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
201 Double *xypos_p = xypos.getStorage( deletePos ) ;
202 const Complex *data_p = spectra.getStorage( deleteData ) ;
203 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
204 //const Int *flag_p = flagI.getStorage( deleteFlag ) ;
205 //const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
206 const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
207 const Int *rflag_p = rflag.getStorage( deleteFlagR ) ;
208 Float *conv_p = convFunc.getStorage( deleteConv ) ;
209 // Extend grid plane with convSupport_
210 Int gnx = nx_ ;
211 Int gny = ny_ ;
212// Int gnx = nx_+convSupport_*2 ;
213// Int gny = ny_+convSupport_*2 ;
214 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
215 Array<Complex> gdataArrC( gshape, 0.0 ) ;
216 Array<Float> gwgtArr( gshape, 0.0 ) ;
217 Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
218 Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
219 Int idopsf = 0 ;
220 Int *chanMap = new Int[nchan_] ;
221 {
222 Int *work_p = chanMap ;
223 for ( Int i = 0 ; i < nchan_ ; i++ ) {
224 *work_p = i ;
225 work_p++ ;
226 }
227 }
228 Int *polMap = new Int[npol_] ;
229 {
230 Int *work_p = polMap ;
231 for ( Int i = 0 ; i < npol_ ; i++ ) {
232 *work_p = i ;
233 work_p++ ;
234 }
235 }
236 Double *sumw_p = new Double[npol_*nchan_] ;
237 {
238 Double *work_p = sumw_p ;
239 for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
240 *work_p = 0.0 ;
241 work_p++ ;
242 }
243 }
244 t0 = mathutil::gettimeofday_sec() ;
245 Int irow = -1 ;
246 ggridsd( xypos_p,
247 data_p,
248 &npol_,
249 &nchan_,
250 &idopsf,
251 flag_p,
252 rflag_p,
253 wgt_p,
254 &nrow_,
255 &irow,
256 gdata_p,
257 wdata_p,
258 &gnx,
259 &gny,
260 &npol_,
261 &nchan_,
262 &convSupport_,
263 &convSampling_,
264 conv_p,
265 chanMap,
266 polMap,
267 sumw_p ) ;
268 t1 = mathutil::gettimeofday_sec() ;
269 os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
270 xypos.putStorage( xypos_p, deletePos ) ;
271 spectra.freeStorage( data_p, deleteData ) ;
272 weight.freeStorage( wgt_p, deleteWgt ) ;
273 flagtra.freeStorage( flag_p, deleteFlag ) ;
274 rflag.freeStorage( rflag_p, deleteFlagR ) ;
275 convFunc.putStorage( conv_p, deleteConv ) ;
276 delete polMap ;
277 delete chanMap ;
278 gdataArrC.putStorage( gdata_p, deleteDataG ) ;
279 gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
280 setData( data_, gdataArrC, gwgtArr ) ;
281 //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
282 delete sumw_p ;
283 //cout << "sumWeight = " << sumWeight << endl ;
284// os << "gdataArr = " << gdataArr << LogIO::POST ;
285// os << "gwgtArr = " << gwgtArr << LogIO::POST ;
286// os << "data_ " << data_ << LogIO::POST ;
287}
288
289void STGrid::setData( Array<Float> &data,
290 Array<Complex> &gdata,
291 Array<Float> &gwgt )
292{
293 LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
294 double t0, t1 ;
295 t0 = mathutil::gettimeofday_sec() ;
296 data.resize( gdata.shape() ) ;
297 uInt len = data.nelements() ;
298 Float *w0_p ;
299 const Complex *w1_p ;
300 const Float *w2_p ;
301 Bool b0, b1, b2 ;
302 Float *data_p = data.getStorage( b0 ) ;
303 const Complex *gdata_p = gdata.getStorage( b1 ) ;
304 const Float *gwgt_p = gwgt.getStorage( b2 ) ;
305 w0_p = data_p ;
306 w1_p = gdata_p ;
307 w2_p = gwgt_p ;
308 for ( uInt i = 0 ; i < len ; i++ ) {
309 *w0_p = (*w2_p > 0.0) ? ((*w1_p).real() / *w2_p) : 0.0 ;
310 w0_p++ ;
311 w1_p++ ;
312 w2_p++ ;
313 }
314 data.putStorage( data_p, b0 ) ;
315 gdata.freeStorage( gdata_p, b1 ) ;
316 gwgt.freeStorage( gwgt_p, b2 ) ;
317 t1 = mathutil::gettimeofday_sec() ;
318 os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
319}
320
321void STGrid::setupGrid( Int &nx,
322 Int &ny,
323 String &cellx,
324 String &celly,
325 Double &xmin,
326 Double &xmax,
327 Double &ymin,
328 Double &ymax,
329 String &center )
330{
331 LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
332 //cout << "nx=" << nx << ", ny=" << ny << endl ;
333
334 // center position
335 if ( center.size() == 0 ) {
336 center_(0) = 0.5 * ( xmin + xmax ) ;
337 center_(1) = 0.5 * ( ymin + ymax ) ;
338 }
339 else {
340 String::size_type pos0 = center.find( " " ) ;
341 if ( pos0 == String::npos ) {
342 throw AipsError( "bad string format in parameter center" ) ;
343 }
344 String::size_type pos1 = center.find( " ", pos0+1 ) ;
345 String typestr, xstr, ystr ;
346 if ( pos1 != String::npos ) {
347 typestr = center.substr( 0, pos0 ) ;
348 xstr = center.substr( pos0+1, pos1-pos0 ) ;
349 ystr = center.substr( pos1+1 ) ;
350 // todo: convert to J2000 (or direction ref for DIRECTION column)
351 }
352 else {
353 typestr = "J2000" ;
354 xstr = center.substr( 0, pos0 ) ;
355 ystr = center.substr( pos0+1 ) ;
356 }
357 QuantumHolder qh ;
358 String err ;
359 qh.fromString( err, xstr ) ;
360 Quantum<Double> xcen = qh.asQuantumDouble() ;
361 qh.fromString( err, ystr ) ;
362 Quantum<Double> ycen = qh.asQuantumDouble() ;
363 center_(0) = xcen.getValue( "rad" ) ;
364 center_(1) = ycen.getValue( "rad" ) ;
365 }
366
367
368 nx_ = nx ;
369 ny_ = ny ;
370 if ( nx < 0 && ny > 0 ) {
371 nx_ = ny ;
372 ny_ = ny ;
373 }
374 if ( ny < 0 && nx > 0 ) {
375 nx_ = nx ;
376 ny_ = nx ;
377 }
378
379 //Double wx = xmax - xmin ;
380 //Double wy = ymax - ymin ;
381 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
382 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
383 // take 10% margin
384 wx *= 1.10 ;
385 wy *= 1.10 ;
386
387 Quantum<Double> qcellx ;
388 Quantum<Double> qcelly ;
389 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
390 if ( cellx.size() != 0 && celly.size() != 0 ) {
391 readQuantity( qcellx, cellx ) ;
392 readQuantity( qcelly, celly ) ;
393 }
394 else if ( celly.size() != 0 ) {
395 os << "Using celly to x-axis..." << LogIO::POST ;
396 readQuantity( qcelly, celly ) ;
397 qcellx = qcelly ;
398 }
399 else if ( cellx.size() != 0 ) {
400 os << "Using cellx to y-axis..." << LogIO::POST ;
401 readQuantity( qcellx, cellx ) ;
402 qcelly = qcellx ;
403 }
404 else {
405 if ( nx_ < 0 ) {
406 os << "No user preference in grid setting. Using default..." << LogIO::POST ;
407 readQuantity( qcellx, "1.0arcmin" ) ;
408 qcelly = qcellx ;
409 }
410 else {
411 if ( wx == 0.0 ) {
412 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
413 wx = 0.00290888 ;
414 }
415 if ( wy == 0.0 ) {
416 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
417 wy = 0.00290888 ;
418 }
419 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
420 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
421 }
422 }
423 cellx_ = qcellx.getValue( "rad" ) ;
424 celly_ = qcelly.getValue( "rad" ) ;
425 if ( nx_ < 0 ) {
426 if ( wx == 0.0 ) {
427 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
428 wx = 0.00290888 ;
429 }
430 if ( wy == 0.0 ) {
431 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
432 wy = 0.00290888 ;
433 }
434 nx_ = Int( ceil( wx/cellx_ ) ) ;
435 ny_ = Int( ceil( wy/celly_ ) ) ;
436 }
437}
438
439void STGrid::selectData( Table &tab )
440{
441 Int ifno = ifno_ ;
442 Table taborg( infile_ ) ;
443 if ( ifno == -1 ) {
444 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
445// os << LogIO::SEVERE
446// << "Please set IFNO before actual gridding"
447// << LogIO::EXCEPTION ;
448 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
449 ifno = ifnoCol( 0 ) ;
450 os << LogIO::WARN
451 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
452 }
453// tab = taborg( taborg.col("IFNO") == ifno ) ;
454 TableExprNode node ;
455 node = taborg.col("IFNO") == ifno ;
456 if ( scanlist_.size() > 0 ) {
457 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
458 }
459// Block<String> cols( 3 ) ;
460// cols[0] = "TIME" ;
461// cols[1] = "BEAMNO" ;
462// cols[2] = "POLNO" ;
463// Block<Int> order( 3, Sort::Ascending ) ;
464// tab = taborg( node ).sort( cols, order ) ;
465 tab = taborg( node ) ;
466 if ( tab.nrow() == 0 ) {
467 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
468 os << LogIO::SEVERE
469 << "No corresponding rows for given selection: IFNO " << ifno ;
470 if ( scanlist_.size() > 0 )
471 os << " SCANNO " << scanlist_ ;
472 os << LogIO::EXCEPTION ;
473 }
474}
475
476void STGrid::getData( Array<Complex> &spectra,
477 Array<Double> &direction,
478 Array<uChar> &flagtra,
479 Array<uInt> &rflag,
480 Array<Float> &weight )
481{
482 LogIO os( LogOrigin("STGrid","getData",WHERE) ) ;
483// os << "start" << LogIO::POST ;
484 Table tab ;
485 selectData( tab ) ;
486 setupArray( tab ) ;
487// os << "npol_ = " << npol_ << LogIO::POST ;
488// os << "nchan_ = " << nchan_ << LogIO::POST ;
489// os << "nrow_ = " << nrow_ << LogIO::POST ;
490 IPosition cshape( 3, npol_, nchan_, nrow_ ) ;
491 IPosition mshape( 2, npol_, nrow_ ) ;
492 spectra.resize( cshape ) ;
493 flagtra.resize( cshape ) ;
494 rflag.resize( mshape ) ;
495 direction.resize( IPosition(2,2,nrow_) ) ;
496 Array<Float> tsys( cshape ) ;
497 Array<Double> tint( mshape ) ;
498
499 ArrayIterator<uChar> fli( flagtra, IPosition(2,1,2) ) ;
500 ArrayIterator<uInt> fri( rflag, IPosition(1,1) ) ;
501 ArrayIterator<Float> tsi( tsys, IPosition(2,1,2) ) ;
502 ArrayIterator<Double> tii( tint, IPosition(1,1) ) ;
503
504 // boolean for pointer access
505 Bool bsp, bsps ;
506 // pointer to the data
507 Complex *sp_p = spectra.getStorage( bsp ) ;
508 // working pointer
509 Complex *wsp_p = sp_p ;
510 uInt len = nchan_ * nrow_ ;
511 IPosition mshape2( 2, nchan_, nrow_ ) ;
512 IPosition vshape( 1, nrow_ ) ;
513 Vector<Float> spSlice( nchan_ ) ;
514 const Float *sps_p = spSlice.getStorage( bsps ) ;
515 long cincr = npol_ ;
516 long rincr = npol_ * nchan_ ;
517 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
518 Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
519 ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
520 ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
521 ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
522 ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
523 ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
524 ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
525 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
526 spectraCol.get( irow, spSlice ) ;
527 const Float *wsps_p = sps_p ;
528 wsp_p = sp_p + (long)ipol + rincr * (long)irow ;
529 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
530 *wsp_p = *wsps_p ;
531 wsps_p++ ;
532 wsp_p += cincr ;
533 }
534 }
535 Array<uChar> flSlice = fli.array() ;
536 Vector<uInt> frSlice = fri.array() ;
537 flagtraCol.getColumn( flSlice ) ;
538 rflagCol.getColumn( frSlice ) ;
539 if ( ipol == 0 )
540 directionCol.getColumn( direction ) ;
541 Vector<Float> tmpF = tsysCol( 0 ) ;
542 Array<Float> tsSlice = tsi.array() ;
543 if ( tmpF.nelements() == (uInt)nchan_ ) {
544 tsysCol.getColumn( tsSlice ) ;
545 }
546 else {
547 tsSlice = tmpF( 0 ) ;
548 }
549 Vector<Double> tmpD = tii.array() ;
550 tintCol.getColumn( tmpD ) ;
551
552 wsp_p += len ;
553
554 fli.next() ;
555 fri.next() ;
556 tsi.next() ;
557 tii.next() ;
558 }
559 spSlice.freeStorage( sps_p, bsps ) ;
560 spectra.putStorage( sp_p, bsp ) ;
561
562// os << "spectra=" << spectra << LogIO::POST ;
563// os << "flagtra=" << flagtra << LogIO::POST ;
564// os << "rflag=" << rflag << LogIO::POST ;
565// os << "direction=" << direction << LogIO::POST ;
566
567 getWeight( weight, tsys, tint ) ;
568}
569
570void STGrid::getData( Array<Complex> &spectra,
571 Array<Double> &direction,
572 Array<Int> &flagtra,
573 Array<Int> &rflag,
574 Array<Float> &weight )
575{
576 LogIO os( LogOrigin("STGrid","getData",WHERE) ) ;
577 double t0, t1 ;
578
579 Array<uChar> flagUC ;
580 Array<uInt> rflagUI ;
581 getData( spectra, direction, flagUC, rflagUI, weight ) ;
582
583 t0 = mathutil::gettimeofday_sec() ;
584 toInt( flagUC, flagtra ) ;
585 toInt( rflagUI, rflag ) ;
586 t1 = mathutil::gettimeofday_sec() ;
587 os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
588}
589
590void STGrid::setupArray( Table &tab )
591{
592 LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
593 ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
594 Vector<uInt> pols = polnoCol.getColumn() ;
595 //os << pols << LogIO::POST ;
596 Vector<uInt> pollistOrg ;
597 uInt npolOrg = 0 ;
598 uInt polno ;
599 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
600 //polno = polnoCol( i ) ;
601 polno = pols( i ) ;
602 if ( allNE( pollistOrg, polno ) ) {
603 pollistOrg.resize( npolOrg+1, True ) ;
604 pollistOrg[npolOrg] = polno ;
605 npolOrg++ ;
606 }
607 }
608 if ( pollist_.size() == 0 )
609 pollist_ = pollistOrg ;
610 else {
611 Vector<uInt> newlist ;
612 uInt newsize = 0 ;
613 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
614 if ( anyEQ( pollistOrg, pollist_[i] ) ) {
615 newlist.resize( newsize+1, True ) ;
616 newlist[newsize] = pollist_[i] ;
617 newsize++ ;
618 }
619 }
620 pollist_.assign( newlist ) ;
621 }
622 npol_ = pollist_.size() ;
623 if ( npol_ == 0 ) {
624 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
625 }
626 nrow_ = tab.nrow() / npolOrg ;
627 ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
628 nchan_ = tmpCol( 0 ).nelements() ;
629// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
630// << "nchan_ = " << nchan_ << endl
631// << "nrow_ = " << nrow_ << LogIO::POST ;
632}
633
634void STGrid::getWeight( Array<Float> &w,
635 Array<Float> &tsys,
636 Array<Double> &tint )
637{
638 LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
639 double t0, t1 ;
640 t0 = mathutil::gettimeofday_sec() ;
641 // resize
642 IPosition refShape = tsys.shape() ;
643 Int nchan = refShape[1] ;
644 Int nrow = refShape[2] ;
645 w.resize( IPosition(2,nchan,nrow) ) ;
646
647 // set weight
648 Bool warn = False ;
649 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
650 w = 1.0 ;
651 }
652 else if ( wtype_.compare( "TINT" ) == 0 ) {
653 if ( npol_ > 1 ) warn = True ;
654 Bool b0, b1 ;
655 Float *w_p = w.getStorage( b0 ) ;
656 Float *w0_p = w_p ;
657 const Double *ti_p = tint.getStorage( b1 ) ;
658 const Double *w1_p = ti_p ;
659 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
660 Float val = (Float)(polMean( w1_p )) ;
661 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
662 *w0_p = val ;
663 w0_p++ ;
664 }
665 }
666 w.putStorage( w_p, b0 ) ;
667 tint.freeStorage( ti_p, b1 ) ;
668 }
669 else if ( wtype_.compare( "TSYS" ) == 0 ) {
670 if ( npol_ > 1 ) warn = True ;
671 Bool b0, b1 ;
672 Float *w_p = w.getStorage( b0 ) ;
673 Float *w0_p = w_p ;
674 const Float *ts_p = tsys.getStorage( b1 ) ;
675 const Float *w1_p = ts_p ;
676 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
677 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
678 Float val = polMean( w1_p ) ;
679 *w0_p = 1.0 / ( val * val ) ;
680 w0_p++ ;
681 }
682 }
683 w.putStorage( w_p, b0 ) ;
684 tsys.freeStorage( ts_p, b1 ) ;
685 }
686 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
687 if ( npol_ > 1 ) warn = True ;
688 Bool b0, b1, b2 ;
689 Float *w_p = w.getStorage( b0 ) ;
690 Float *w0_p = w_p ;
691 const Double *ti_p = tint.getStorage( b1 ) ;
692 const Double *w1_p = ti_p ;
693 const Float *ts_p = tsys.getStorage( b2 ) ;
694 const Float *w2_p = ts_p ;
695 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
696 Float interval = (Float)(polMean( w1_p )) ;
697 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
698 Float temp = polMean( w2_p ) ;
699 *w0_p = interval / ( temp * temp ) ;
700 w0_p++ ;
701 }
702 }
703 w.putStorage( w_p, b0 ) ;
704 tint.freeStorage( ti_p, b1 ) ;
705 tsys.freeStorage( ts_p, b2 ) ;
706 }
707 else {
708 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
709 os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
710 w = 1.0 ;
711 }
712
713 if ( npol_ > 1 ) {
714 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
715 os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
716 }
717 t1 = mathutil::gettimeofday_sec() ;
718 os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
719}
720
721Float STGrid::polMean( const Float *p )
722{
723 Float v = 0.0 ;
724 for ( Int i = 0 ; i < npol_ ; i++ ) {
725 v += *p ;
726 p++ ;
727 }
728 v /= npol_ ;
729 return v ;
730}
731
732Double STGrid::polMean( const Double *p )
733{
734 Double v = 0.0 ;
735 for ( Int i = 0 ; i < npol_ ; i++ ) {
736 v += *p ;
737 p++ ;
738 }
739 v /= npol_ ;
740 return v ;
741}
742
743void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
744{
745 uInt len = u.nelements() ;
746 Int *int_p = new Int[len] ;
747 Bool deleteIt ;
748 const uChar *data_p = u.getStorage( deleteIt ) ;
749 Int *i_p = int_p ;
750 const uChar *u_p = data_p ;
751 for ( uInt i = 0 ; i < len ; i++ ) {
752 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
753 i_p++ ;
754 u_p++ ;
755 }
756 u.freeStorage( data_p, deleteIt ) ;
757 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
758}
759
760void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
761{
762 uInt len = u.nelements() ;
763 Int *int_p = new Int[len] ;
764 Bool deleteIt ;
765 const uInt *data_p = u.getStorage( deleteIt ) ;
766 Int *i_p = int_p ;
767 const uInt *u_p = data_p ;
768 for ( uInt i = 0 ; i < len ; i++ ) {
769 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
770 i_p++ ;
771 u_p++ ;
772 }
773 u.freeStorage( data_p, deleteIt ) ;
774 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
775}
776
777void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
778{
779 // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
780 // grid plane to avoid unexpected behavior on grid edge
781 Block<Double> pixc( 2 ) ;
782 pixc[0] = Double( nx_-1 ) * 0.5 ;
783 pixc[1] = Double( ny_-1 ) * 0.5 ;
784// pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
785// pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
786 uInt nrow = world.shape()[1] ;
787 Bool bw, bp ;
788 const Double *w_p = world.getStorage( bw ) ;
789 Double *p_p = pixel.getStorage( bp ) ;
790 const Double *ww_p = w_p ;
791 Double *wp_p = p_p ;
792 for ( uInt i = 0 ; i < nrow ; i++ ) {
793 *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
794 wp_p++ ;
795 ww_p++ ;
796 *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
797 wp_p++ ;
798 ww_p++ ;
799 }
800 world.freeStorage( w_p, bw ) ;
801 pixel.putStorage( p_p, bp ) ;
802// String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
803// ofstream ofs( gridfile.c_str(), ios::out ) ;
804// ofs << "center " << center_(0) << " " << pixc(0)
805// << " " << center_(1) << " " << pixc(1) << endl ;
806// for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
807// ofs << irow ;
808// for ( uInt i = 0 ; i < 2 ; i++ ) {
809// ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
810// }
811// ofs << endl ;
812// }
813// ofs.close() ;
814}
815
816void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
817{
818 convFunc = 0.0 ;
819 for ( Int i = 0 ; i < convSize/2 ; i++ )
820 convFunc(i) = 1.0 ;
821}
822
823#define NEED_UNDERSCORES
824#if defined(NEED_UNDERSCORES)
825#define grdsf grdsf_
826#endif
827extern "C" {
828 void grdsf(Double*, Double*);
829}
830void STGrid::spheroidalFunc( Vector<Float> &convFunc )
831{
832 convFunc = 0.0 ;
833 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
834 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
835 Double val ;
836 grdsf( &nu, &val ) ;
837 convFunc(i) = ( 1.0 - nu * nu ) * val ;
838 }
839}
840
841void STGrid::gaussFunc( Vector<Float> &convFunc )
842{
843 convFunc = 0.0 ;
844 // HWHM of the Gaussian is convSupport_ / 4
845 // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
846 Int len = convSampling_ * convSupport_ ;
847 Double hwhm = len * 0.25 ;
848 for ( Int i = 0 ; i < len ; i++ ) {
849 Double val = Double(i) / hwhm ;
850 convFunc(i) = exp( -log(2)*val*val ) ;
851 }
852}
853
854void STGrid::pbFunc( Vector<Float> &convFunc )
855{
856 convFunc = 0.0 ;
857}
858
859void STGrid::setConvFunc( Vector<Float> &convFunc )
860{
861 convSupport_ = userSupport_ ;
862 if ( convType_ == "BOX" ) {
863 if ( convSupport_ < 0 )
864 convSupport_ = 0 ;
865 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
866 convFunc.resize( convSize ) ;
867 boxFunc( convFunc, convSize ) ;
868 }
869 else if ( convType_ == "SF" ) {
870 if ( convSupport_ < 0 )
871 convSupport_ = 3 ;
872 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
873 convFunc.resize( convSize ) ;
874 spheroidalFunc( convFunc ) ;
875 }
876 else if ( convType_ == "GAUSS" ) {
877 // to take into account Gaussian tail
878 if ( convSupport_ < 0 )
879 convSupport_ = 12 ; // 3 * 4
880 else {
881 convSupport_ = userSupport_ * 4 ;
882 }
883 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
884 convFunc.resize( convSize ) ;
885 gaussFunc( convFunc ) ;
886 }
887 else if ( convType_ == "PB" ) {
888 if ( convSupport_ < 0 )
889 convSupport_ = 0 ;
890 pbFunc( convFunc ) ;
891 }
892 else {
893 throw AipsError( "Unsupported convolution function" ) ;
894 }
895}
896
897string STGrid::saveData( string outfile )
898{
899 LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
900 double t0, t1 ;
901 t0 = mathutil::gettimeofday_sec() ;
902
903 //Int polno = 0 ;
904 String outfile_ ;
905 if ( outfile.size() == 0 ) {
906 if ( infile_.lastchar() == '/' ) {
907 outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
908 }
909 else {
910 outfile_ = infile_ ;
911 }
912 outfile_ += ".grid" ;
913 }
914 else {
915 outfile_ = outfile ;
916 }
917 Table tab ;
918 prepareTable( tab, outfile_ ) ;
919 IPosition dshape = data_.shape() ;
920 Int nrow = nx_ * ny_ * npol_ ;
921 tab.rwKeywordSet().define( "nPol", npol_ ) ;
922 tab.addRow( nrow ) ;
923 Vector<Double> cpix( 2 ) ;
924 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
925 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
926 Vector<Double> dir( 2 ) ;
927 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
928 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
929 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
930 Int irow = 0 ;
931 Vector<Float> sp( nchan_ ) ;
932 Bool bsp, bdata ;
933 const Float *data_p = data_.getStorage( bdata ) ;
934 Float *wsp_p, *sp_p ;
935 const Float *wdata_p = data_p ;
936 long step = nx_ * ny_ * npol_ ;
937 long offset ;
938 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
939 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
940 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
941 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
942 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
943 offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
944 //os << "offset = " << offset << LogIO::POST ;
945 sp_p = sp.getStorage( bsp ) ;
946 wsp_p = sp_p ;
947 wdata_p = data_p + offset ;
948 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
949 *wsp_p = *wdata_p ;
950 wsp_p++ ;
951 wdata_p += step ;
952 }
953 sp.putStorage( sp_p, bsp ) ;
954 spectraCol.put( irow, sp ) ;
955 directionCol.put( irow, dir ) ;
956 polnoCol.put( irow, pollist_[ipol] ) ;
957 irow++ ;
958 }
959 }
960 }
961 data_.freeStorage( data_p, bdata ) ;
962
963 t1 = mathutil::gettimeofday_sec() ;
964 os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
965
966 return outfile_ ;
967}
968
969void STGrid::prepareTable( Table &tab, String &name )
970{
971 Table t( infile_, Table::Old ) ;
972 t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
973 tab = Table( name, Table::Update ) ;
974}
975}
Note: See TracBrowser for help on using the repository browser.