source: trunk/src/STGrid.cpp@ 2373

Last change on this file since 2373 was 2371, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: No

Interface Changes: 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...

Direct array data access in STGrid::getData and STGrid::saveData.
Modified duplicate table copy in STGrid::saveData.
File size: 25.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/ArrayPartMath.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 Cube<Float> spectra ;
158 Matrix<Double> direction ;
159 Cube<uChar> flagtra ;
160 Matrix<uInt> rflag ;
161 Matrix<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 // flagtra: uChar -> Int
173 // rflag: uInt -> Int
174 Cube<Int> flagI ;
175 Matrix<Int> rflagI ;
176 t0 = mathutil::gettimeofday_sec() ;
177 toInt( &flagtra, &flagI ) ;
178 toInt( &rflag, &rflagI ) ;
179 t1 = mathutil::gettimeofday_sec() ;
180 os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
181
182 // grid parameter
183 os << LogIO::DEBUGGING ;
184 os << "----------" << endl ;
185 os << "Grid parameter summary" << endl ;
186 os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
187 os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
188 os << " center = " << center_ << endl ;
189 os << "----------" << LogIO::POST ;
190 os << LogIO::NORMAL ;
191
192 // convolution kernel
193 Vector<Float> convFunc ;
194 t0 = mathutil::gettimeofday_sec() ;
195 setConvFunc( convFunc ) ;
196 t1 = mathutil::gettimeofday_sec() ;
197 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
198 //cout << "convSupport=" << convSupport_ << endl ;
199 //cout << "convFunc=" << convFunc << endl ;
200
201 // world -> pixel
202 Matrix<Double> xypos( direction.shape(), 0.0 ) ;
203 t0 = mathutil::gettimeofday_sec() ;
204 toPixel( direction, xypos ) ;
205 t1 = mathutil::gettimeofday_sec() ;
206 os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
207
208 // call ggridsd
209 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
210 Double *xypos_p = xypos.getStorage( deletePos ) ;
211 Cube<Complex> dataC( spectra.shape(), 0.0 ) ;
212 setReal( dataC, spectra ) ;
213 const Complex *data_p = dataC.getStorage( deleteData ) ;
214 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
215 const Int *flag_p = flagI.getStorage( deleteFlag ) ;
216 const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
217 Float *conv_p = convFunc.getStorage( deleteConv ) ;
218 // Extend grid plane with convSupport_
219 Int gnx = nx_ ;
220 Int gny = ny_ ;
221// Int gnx = nx_+convSupport_*2 ;
222// Int gny = ny_+convSupport_*2 ;
223 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
224 Array<Complex> gdataArrC( gshape, 0.0 ) ;
225 Array<Float> gwgtArr( gshape, 0.0 ) ;
226 Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
227 Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
228 Int idopsf = 0 ;
229 Int *chanMap = new Int[nchan_] ;
230 {
231 Int *work_p = chanMap ;
232 for ( Int i = 0 ; i < nchan_ ; i++ ) {
233 *work_p = i ;
234 work_p++ ;
235 }
236 }
237 Int *polMap = new Int[npol_] ;
238 {
239 Int *work_p = polMap ;
240 for ( Int i = 0 ; i < npol_ ; i++ ) {
241 *work_p = i ;
242 work_p++ ;
243 }
244 }
245 Double *sumw_p = new Double[npol_*nchan_] ;
246 {
247 Double *work_p = sumw_p ;
248 for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
249 *work_p = 0.0 ;
250 work_p++ ;
251 }
252 }
253 t0 = mathutil::gettimeofday_sec() ;
254 Int irow = -1 ;
255 ggridsd( xypos_p,
256 data_p,
257 &npol_,
258 &nchan_,
259 &idopsf,
260 flag_p,
261 rflag_p,
262 wgt_p,
263 &nrow_,
264 &irow,
265 gdata_p,
266 wdata_p,
267 &gnx,
268 &gny,
269 &npol_,
270 &nchan_,
271 &convSupport_,
272 &convSampling_,
273 conv_p,
274 chanMap,
275 polMap,
276 sumw_p ) ;
277 t1 = mathutil::gettimeofday_sec() ;
278 os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
279 xypos.putStorage( xypos_p, deletePos ) ;
280 dataC.freeStorage( data_p, deleteData ) ;
281 weight.freeStorage( wgt_p, deleteWgt ) ;
282 flagI.freeStorage( flag_p, deleteFlag ) ;
283 rflagI.freeStorage( rflag_p, deleteFlagR ) ;
284 convFunc.putStorage( conv_p, deleteConv ) ;
285 delete polMap ;
286 delete chanMap ;
287 gdataArrC.putStorage( gdata_p, deleteDataG ) ;
288 gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
289 Array<Float> gdataArr = real( gdataArrC ) ;
290 setData( data_, gdataArr, gwgtArr ) ;
291 //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
292 delete sumw_p ;
293 //cout << "sumWeight = " << sumWeight << endl ;
294// os << "gdataArr = " << gdataArr << LogIO::POST ;
295// os << "gwgtArr = " << gwgtArr << LogIO::POST ;
296// os << "data_ " << data_ << LogIO::POST ;
297}
298
299void STGrid::setData( Array<Float> &data,
300 Array<Float> &gdata,
301 Array<Float> &gwgt )
302{
303 LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
304 double t0, t1 ;
305 t0 = mathutil::gettimeofday_sec() ;
306 data.resize( gdata.shape() ) ;
307 uInt len = data.nelements() ;
308 Float *w0_p ;
309 const Float *w1_p, *w2_p ;
310 Bool b0, b1, b2 ;
311 Float *data_p = data.getStorage( b0 ) ;
312 const Float *gdata_p = gdata.getStorage( b1 ) ;
313 const Float *gwgt_p = gwgt.getStorage( b2 ) ;
314 w0_p = data_p ;
315 w1_p = gdata_p ;
316 w2_p = gwgt_p ;
317 for ( uInt i = 0 ; i < len ; i++ ) {
318 *w0_p = (*w2_p > 0.0) ? (*w1_p / *w2_p) : 0.0 ;
319 w0_p++ ;
320 w1_p++ ;
321 w2_p++ ;
322 }
323 data.putStorage( data_p, b0 ) ;
324 gdata.freeStorage( gdata_p, b1 ) ;
325 gwgt.freeStorage( gwgt_p, b2 ) ;
326 t1 = mathutil::gettimeofday_sec() ;
327 os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
328}
329
330void STGrid::setupGrid( Int &nx,
331 Int &ny,
332 String &cellx,
333 String &celly,
334 Double &xmin,
335 Double &xmax,
336 Double &ymin,
337 Double &ymax,
338 String &center )
339{
340 //cout << "nx=" << nx << ", ny=" << ny << endl ;
341
342 // center position
343 if ( center.size() == 0 ) {
344 center_(0) = 0.5 * ( xmin + xmax ) ;
345 center_(1) = 0.5 * ( ymin + ymax ) ;
346 }
347 else {
348 String::size_type pos0 = center.find( " " ) ;
349 if ( pos0 == String::npos ) {
350 throw AipsError( "bad string format in parameter center" ) ;
351 }
352 String::size_type pos1 = center.find( " ", pos0+1 ) ;
353 String typestr, xstr, ystr ;
354 if ( pos1 != String::npos ) {
355 typestr = center.substr( 0, pos0 ) ;
356 xstr = center.substr( pos0+1, pos1-pos0 ) ;
357 ystr = center.substr( pos1+1 ) ;
358 // todo: convert to J2000 (or direction ref for DIRECTION column)
359 }
360 else {
361 typestr = "J2000" ;
362 xstr = center.substr( 0, pos0 ) ;
363 ystr = center.substr( pos0+1 ) ;
364 }
365 QuantumHolder qh ;
366 String err ;
367 qh.fromString( err, xstr ) ;
368 Quantum<Double> xcen = qh.asQuantumDouble() ;
369 qh.fromString( err, ystr ) ;
370 Quantum<Double> ycen = qh.asQuantumDouble() ;
371 center_(0) = xcen.getValue( "rad" ) ;
372 center_(1) = ycen.getValue( "rad" ) ;
373 }
374
375
376 //Double wx = xmax - xmin ;
377 //Double wy = ymax - ymin ;
378 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
379 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
380 // take 10% margin
381 wx *= 1.10 ;
382 wy *= 1.10 ;
383 Quantum<Double> qcellx ;
384 Quantum<Double> qcelly ;
385 nx_ = nx ;
386 ny_ = ny ;
387 if ( nx < 0 && ny > 0 ) {
388 nx_ = ny ;
389 ny_ = ny ;
390 }
391 if ( ny < 0 && nx > 0 ) {
392 nx_ = nx ;
393 ny_ = nx ;
394 }
395 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
396 if ( cellx.size() != 0 && celly.size() != 0 ) {
397 readQuantity( qcellx, cellx ) ;
398 readQuantity( qcelly, celly ) ;
399 }
400 else if ( celly.size() != 0 ) {
401 cout << "Using celly to x-axis..." << endl ;
402 readQuantity( qcelly, celly ) ;
403 qcellx = qcelly ;
404 }
405 else if ( cellx.size() != 0 ) {
406 cout << "Using cellx to y-axis..." << endl ;
407 readQuantity( qcellx, cellx ) ;
408 qcelly = qcellx ;
409 }
410 else {
411 if ( nx_ < 0 ) {
412 cout << "No user preference in grid setting. Using default..." << endl ;
413 readQuantity( qcellx, "1.0arcmin" ) ;
414 qcelly = qcellx ;
415 }
416 else {
417 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
418 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
419 }
420 }
421 cellx_ = qcellx.getValue( "rad" ) ;
422 celly_ = qcelly.getValue( "rad" ) ;
423 if ( nx_ < 0 ) {
424 nx_ = Int( ceil( wx/cellx_ ) ) ;
425 ny_ = Int( ceil( wy/celly_ ) ) ;
426 }
427}
428
429void STGrid::selectData( Table &tab )
430{
431 Int ifno = ifno_ ;
432 Table taborg( infile_ ) ;
433 if ( ifno == -1 ) {
434 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
435// os << LogIO::SEVERE
436// << "Please set IFNO before actual gridding"
437// << LogIO::EXCEPTION ;
438 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
439 ifno = ifnoCol( 0 ) ;
440 os << LogIO::WARN
441 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
442 }
443// tab = taborg( taborg.col("IFNO") == ifno ) ;
444 TableExprNode node ;
445 node = taborg.col("IFNO") == ifno ;
446 if ( scanlist_.size() > 0 ) {
447 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
448 }
449 tab = taborg( node ) ;
450 if ( tab.nrow() == 0 ) {
451 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
452 os << LogIO::SEVERE
453 << "No corresponding rows for given selection: IFNO " << ifno
454 << " SCANNO " << scanlist_
455 << LogIO::EXCEPTION ;
456 }
457}
458
459void STGrid::getData( Cube<Float> &spectra,
460 Matrix<Double> &direction,
461 Cube<uChar> &flagtra,
462 Matrix<uInt> &rflag,
463 Matrix<Float> &weight )
464{
465 Table tab ;
466 selectData( tab ) ;
467 updatePolList( tab ) ;
468// cout << "npol_ = " << npol_ << endl ;
469// cout << "nchan_ = " << nchan_ << endl ;
470// cout << "nrow_ = " << nrow_ << endl ;
471 spectra.resize( npol_, nchan_, nrow_ ) ;
472 flagtra.resize( npol_, nchan_, nrow_ ) ;
473 rflag.resize( npol_, nrow_ ) ;
474 Cube<Float> tsys( npol_, nchan_, nrow_ ) ;
475 Matrix<Double> tint( npol_, nrow_ ) ;
476 // boolean for pointer access
477 Bool bsp, bfl, bfr, bts, bti ;
478 // pointer to the data
479 Float *sp_p = spectra.getStorage( bsp ) ;
480 uChar *fl_p = flagtra.getStorage( bfl ) ;
481 uInt *fr_p = rflag.getStorage( bfr ) ;
482 Float *ts_p = tsys.getStorage( bts ) ;
483 Double *ti_p = tint.getStorage( bti ) ;
484 // working pointer
485 Float *wsp_p = sp_p ;
486 uChar *wfl_p = fl_p ;
487 uInt *wfr_p = fr_p ;
488 Float *wts_p = ts_p ;
489 Double *wti_p = ti_p ;
490 uInt len = nchan_ * nrow_ ;
491 IPosition mshape( 2, nchan_, nrow_ ) ;
492 IPosition vshape( 1, nrow_ ) ;
493 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
494 Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
495 ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
496 ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
497 ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
498 ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
499 ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
500 ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
501 Matrix<Float> spSlice( mshape, wsp_p, SHARE ) ;
502 Matrix<uChar> flSlice( mshape, wfl_p, SHARE ) ;
503 Vector<uInt> frSlice( vshape, wfr_p, SHARE ) ;
504 spectraCol.getColumn( spSlice ) ;
505 flagtraCol.getColumn( flSlice ) ;
506 rflagCol.getColumn( frSlice ) ;
507 if ( ipol == 0 )
508 directionCol.getColumn( direction ) ;
509 Vector<Float> tmpF = tsysCol( 0 ) ;
510 Vector<Double> tmpD( vshape, wti_p, SHARE ) ;
511 Matrix<Float> tsSlice( mshape, wts_p, SHARE ) ;
512 if ( tmpF.nelements() == (uInt)nchan_ ) {
513 tsysCol.getColumn( tsSlice ) ;
514 }
515 else {
516 tsSlice = tmpF( 0 ) ;
517 }
518 tintCol.getColumn( tmpD ) ;
519
520 wsp_p += len ;
521 wfl_p += len ;
522 wfr_p += nrow_ ;
523 wts_p += len ;
524 wti_p += nrow_ ;
525 }
526 spectra.putStorage( sp_p, bsp ) ;
527 flagtra.putStorage( fl_p, bfl ) ;
528 rflag.putStorage( fr_p, bfr ) ;
529 tsys.putStorage( ts_p, bts ) ;
530 tint.putStorage( ti_p, bti ) ;
531
532 getWeight( weight, tsys, tint ) ;
533}
534
535void STGrid::updatePolList( Table &tab )
536{
537 ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
538 Vector<uInt> pols = polnoCol.getColumn() ;
539 Vector<uInt> pollistOrg ;
540 uInt npolOrg = 0 ;
541 uInt polno ;
542 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
543 //polno = polnoCol( i ) ;
544 polno = pols( i ) ;
545 if ( allNE( pollistOrg, polno ) ) {
546 pollistOrg.resize( npolOrg+1, True ) ;
547 pollistOrg[npolOrg] = polno ;
548 npolOrg++ ;
549 }
550 }
551 if ( pollist_.size() == 0 )
552 pollist_ = pollistOrg ;
553 else {
554 Vector<uInt> newlist ;
555 uInt newsize = 0 ;
556 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
557 if ( anyEQ( pollistOrg, pollist_[i] ) ) {
558 newlist.resize( newsize+1, True ) ;
559 newlist[newsize] = pollist_[i] ;
560 newsize++ ;
561 }
562 }
563 pollist_.assign( newlist ) ;
564 }
565 npol_ = pollist_.size() ;
566 if ( npol_ == 0 ) {
567 LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
568 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
569 }
570 nrow_ = tab.nrow() / npolOrg ;
571 ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
572 nchan_ = tmpCol( 0 ).nelements() ;
573// LogIO os( LogOrigin("STGrid","updatePolList",WHERE) ) ;
574// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
575// << "nchan_ = " << nchan_ << endl
576// << "nrow_ = " << nrow_ << LogIO::POST ;
577}
578
579void STGrid::getWeight( Matrix<Float> &w,
580 Cube<Float> &tsys,
581 Matrix<Double> &tint )
582{
583 LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
584 double t0, t1 ;
585 t0 = mathutil::gettimeofday_sec() ;
586 // resize
587 w.resize( nchan_, nrow_ ) ;
588
589 // set weight
590 Bool warn = False ;
591 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
592 w = 1.0 ;
593 }
594 else if ( wtype_.compare( "TINT" ) == 0 ) {
595 if ( npol_ > 1 ) warn = True ;
596 Bool b0, b1 ;
597 Float *w_p = w.getStorage( b0 ) ;
598 Float *w0_p = w_p ;
599 const Double *ti_p = tint.getStorage( b1 ) ;
600 const Double *w1_p = ti_p ;
601 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
602 Float val = (Float)(polMean( w1_p )) ;
603 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
604 *w0_p = val ;
605 w0_p++ ;
606 }
607 }
608 w.putStorage( w_p, b0 ) ;
609 tint.freeStorage( ti_p, b1 ) ;
610 }
611 else if ( wtype_.compare( "TSYS" ) == 0 ) {
612 if ( npol_ > 1 ) warn = True ;
613 Bool b0, b1 ;
614 Float *w_p = w.getStorage( b0 ) ;
615 Float *w0_p = w_p ;
616 const Float *ts_p = tsys.getStorage( b1 ) ;
617 const Float *w1_p = ts_p ;
618 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
619 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
620 Float val = polMean( w1_p ) ;
621 *w0_p = 1.0 / ( val * val ) ;
622 w0_p++ ;
623 }
624 }
625 w.putStorage( w_p, b0 ) ;
626 tsys.freeStorage( ts_p, b1 ) ;
627 }
628 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
629 if ( npol_ > 1 ) warn = True ;
630 Bool b0, b1, b2 ;
631 Float *w_p = w.getStorage( b0 ) ;
632 Float *w0_p = w_p ;
633 const Double *ti_p = tint.getStorage( b1 ) ;
634 const Double *w1_p = ti_p ;
635 const Float *ts_p = tsys.getStorage( b2 ) ;
636 const Float *w2_p = ts_p ;
637 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
638 Float interval = (Float)(polMean( w1_p )) ;
639 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
640 Float temp = polMean( w2_p ) ;
641 *w0_p = interval / ( temp * temp ) ;
642 w0_p++ ;
643 }
644 }
645 w.putStorage( w_p, b0 ) ;
646 tint.freeStorage( ti_p, b1 ) ;
647 tsys.freeStorage( ts_p, b2 ) ;
648 }
649 else {
650 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
651 os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
652 w = 1.0 ;
653 }
654
655 if ( npol_ > 1 ) {
656 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
657 os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
658 }
659 t1 = mathutil::gettimeofday_sec() ;
660 os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
661}
662
663Float STGrid::polMean( const Float *p )
664{
665 Float v = 0.0 ;
666 for ( Int i = 0 ; i < npol_ ; i++ ) {
667 v += *p ;
668 p++ ;
669 }
670 v /= npol_ ;
671 return v ;
672}
673
674Double STGrid::polMean( const Double *p )
675{
676 Double v = 0.0 ;
677 for ( Int i = 0 ; i < npol_ ; i++ ) {
678 v += *p ;
679 p++ ;
680 }
681 v /= npol_ ;
682 return v ;
683}
684
685void STGrid::toInt( Array<uChar> *u, Array<Int> *v )
686{
687 uInt len = u->nelements() ;
688 Int *int_p = new Int[len] ;
689 Bool deleteIt ;
690 const uChar *data_p = u->getStorage( deleteIt ) ;
691 Int *i_p = int_p ;
692 const uChar *u_p = data_p ;
693 for ( uInt i = 0 ; i < len ; i++ ) {
694 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
695 i_p++ ;
696 u_p++ ;
697 }
698 u->freeStorage( data_p, deleteIt ) ;
699 v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
700}
701
702void STGrid::toInt( Array<uInt> *u, Array<Int> *v )
703{
704 uInt len = u->nelements() ;
705 Int *int_p = new Int[len] ;
706 Bool deleteIt ;
707 const uInt *data_p = u->getStorage( deleteIt ) ;
708 Int *i_p = int_p ;
709 const uInt *u_p = data_p ;
710 for ( uInt i = 0 ; i < len ; i++ ) {
711 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
712 i_p++ ;
713 u_p++ ;
714 }
715 u->freeStorage( data_p, deleteIt ) ;
716 v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
717}
718
719void STGrid::toPixel( Matrix<Double> &world, Matrix<Double> &pixel )
720{
721 // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
722 // grid plane to avoid unexpected behavior on grid edge
723 Vector<Double> pixc( 2 ) ;
724 pixc(0) = Double( nx_-1 ) * 0.5 ;
725 pixc(1) = Double( ny_-1 ) * 0.5 ;
726// pixc(0) = Double( nx_+2*convSupport_-1 ) * 0.5 ;
727// pixc(1) = Double( ny_+2*convSupport_-1 ) * 0.5 ;
728 uInt nrow = world.shape()[1] ;
729 Vector<Double> cell( 2 ) ;
730 cell(0) = cellx_ ;
731 cell(1) = celly_ ;
732 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
733 for ( uInt i = 0 ; i < 2 ; i++ ) {
734 pixel( i, irow ) = pixc(i) + ( world(i, irow) - center_(i) ) / cell(i) ;
735 }
736 }
737// String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
738// ofstream ofs( gridfile.c_str(), ios::out ) ;
739// ofs << "center " << center_(0) << " " << pixc(0)
740// << " " << center_(1) << " " << pixc(1) << endl ;
741// for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
742// ofs << irow ;
743// for ( uInt i = 0 ; i < 2 ; i++ ) {
744// ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
745// }
746// ofs << endl ;
747// }
748// ofs.close() ;
749}
750
751void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
752{
753 convFunc = 0.0 ;
754 for ( Int i = 0 ; i < convSize/2 ; i++ )
755 convFunc(i) = 1.0 ;
756}
757
758#define NEED_UNDERSCORES
759#if defined(NEED_UNDERSCORES)
760#define grdsf grdsf_
761#endif
762extern "C" {
763 void grdsf(Double*, Double*);
764}
765void STGrid::spheroidalFunc( Vector<Float> &convFunc )
766{
767 convFunc = 0.0 ;
768 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
769 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
770 Double val ;
771 grdsf( &nu, &val ) ;
772 convFunc(i) = ( 1.0 - nu * nu ) * val ;
773 }
774}
775
776void STGrid::gaussFunc( Vector<Float> &convFunc )
777{
778 convFunc = 0.0 ;
779 // HWHM of the Gaussian is convSupport_ / 4
780 // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
781 Int len = convSampling_ * convSupport_ ;
782 Double hwhm = len * 0.25 ;
783 for ( Int i = 0 ; i < len ; i++ ) {
784 Double val = Double(i) / hwhm ;
785 convFunc(i) = exp( -log(2)*val*val ) ;
786 }
787}
788
789void STGrid::pbFunc( Vector<Float> &convFunc )
790{
791 convFunc = 0.0 ;
792}
793
794void STGrid::setConvFunc( Vector<Float> &convFunc )
795{
796 convSupport_ = userSupport_ ;
797 if ( convType_ == "BOX" ) {
798 if ( convSupport_ < 0 )
799 convSupport_ = 0 ;
800 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
801 convFunc.resize( convSize ) ;
802 boxFunc( convFunc, convSize ) ;
803 }
804 else if ( convType_ == "SF" ) {
805 if ( convSupport_ < 0 )
806 convSupport_ = 3 ;
807 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
808 convFunc.resize( convSize ) ;
809 spheroidalFunc( convFunc ) ;
810 }
811 else if ( convType_ == "GAUSS" ) {
812 // to take into account Gaussian tail
813 if ( convSupport_ < 0 )
814 convSupport_ = 12 ; // 3 * 4
815 else {
816 convSupport_ = userSupport_ * 4 ;
817 }
818 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
819 convFunc.resize( convSize ) ;
820 gaussFunc( convFunc ) ;
821 }
822 else if ( convType_ == "PB" ) {
823 if ( convSupport_ < 0 )
824 convSupport_ = 0 ;
825 pbFunc( convFunc ) ;
826 }
827 else {
828 throw AipsError( "Unsupported convolution function" ) ;
829 }
830}
831
832string STGrid::saveData( string outfile )
833{
834 LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
835 double t0, t1 ;
836 t0 = mathutil::gettimeofday_sec() ;
837
838 //Int polno = 0 ;
839 String outfile_ ;
840 if ( outfile.size() == 0 ) {
841 if ( infile_.lastchar() == '/' ) {
842 outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
843 }
844 else {
845 outfile_ = infile_ ;
846 }
847 outfile_ += ".grid" ;
848 }
849 else {
850 outfile_ = outfile ;
851 }
852 Table tab ;
853 prepareTable( tab, outfile_ ) ;
854 IPosition dshape = data_.shape() ;
855 Int nrow = nx_ * ny_ * npol_ ;
856 tab.rwKeywordSet().define( "nPol", npol_ ) ;
857 tab.addRow( nrow ) ;
858 Vector<Double> cpix( 2 ) ;
859 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
860 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
861 Vector<Double> dir( 2 ) ;
862 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
863 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
864 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
865 Int irow = 0 ;
866 Vector<Float> sp( nchan_ ) ;
867 Bool bsp, bdata ;
868 const Float *data_p = data_.getStorage( bdata ) ;
869 Float *wsp_p, *sp_p ;
870 const Float *wdata_p = data_p ;
871 long step = nx_ * ny_ * npol_ ;
872 long offset ;
873 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
874 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
875 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
876 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
877 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
878 offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
879 //os << "offset = " << offset << LogIO::POST ;
880 sp_p = sp.getStorage( bsp ) ;
881 wsp_p = sp_p ;
882 wdata_p = data_p + offset ;
883 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
884 *wsp_p = *wdata_p ;
885 wsp_p++ ;
886 wdata_p += step ;
887 }
888 sp.putStorage( sp_p, bsp ) ;
889 spectraCol.put( irow, sp ) ;
890 directionCol.put( irow, dir ) ;
891 polnoCol.put( irow, pollist_[ipol] ) ;
892 irow++ ;
893 }
894 }
895 }
896 data_.freeStorage( data_p, bdata ) ;
897
898 t1 = mathutil::gettimeofday_sec() ;
899 os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
900
901 return outfile_ ;
902}
903
904void STGrid::prepareTable( Table &tab, String &name )
905{
906 Table t( infile_, Table::Old ) ;
907 t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
908 tab = Table( name, Table::Update ) ;
909}
910}
Note: See TracBrowser for help on using the repository browser.