source: trunk/src/STGrid.cpp@ 2367

Last change on this file since 2367 was 2365, 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...

Warning message modified and plot style is changed.


File size: 21.9 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 <Scantable.h>
24#include <MathUtils.h>
25
26#include "STGrid.h"
27
28using namespace std ;
29using namespace casa ;
30using namespace asap ;
31
32namespace asap {
33
34// constructor
35STGrid::STGrid()
36{
37 init() ;
38}
39
40STGrid::STGrid( const string infile )
41{
42 init() ;
43
44 setFileIn( infile ) ;
45}
46
47void STGrid::init()
48{
49 ifno_ = -1 ;
50 nx_ = -1 ;
51 ny_ = -1 ;
52 npol_ = 0 ;
53 nchan_ = 0 ;
54 nrow_ = 0 ;
55 cellx_ = 0.0 ;
56 celly_ = 0.0 ;
57 center_ = Vector<Double> ( 2, 0.0 ) ;
58 convType_ = "BOX" ;
59 wtype_ = "UNIFORM" ;
60 convSupport_ = -1 ;
61 userSupport_ = -1 ;
62 convSampling_ = 100 ;
63}
64
65void STGrid::setFileIn( const string infile )
66{
67 String name( infile ) ;
68 if ( infile_.compare( name ) != 0 ) {
69 infile_ = String( infile ) ;
70 tab_ = Table( infile_ ) ;
71 }
72}
73
74void STGrid::setPolList( vector<unsigned int> pols )
75{
76 pollist_.assign( Vector<uInt>( pols ) ) ;
77 cout << "pollist_ = " << pollist_ << endl ;
78}
79
80void STGrid::setScanList( vector<unsigned int> scans )
81{
82 scanlist_.assign( Vector<uInt>( scans ) ) ;
83 cout << "scanlist_ = " << scanlist_ << endl ;
84}
85
86void STGrid::setWeight( const string wType )
87{
88 wtype_ = String( wType ) ;
89 wtype_.upcase() ;
90 cout << "wtype_ = " << wtype_ << endl ;
91}
92
93void STGrid::defineImage( int nx,
94 int ny,
95 string scellx,
96 string scelly,
97 string scenter )
98{
99 ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
100 Matrix<Double> direction = dirCol.getColumn() ;
101 Double rmax, rmin, dmax, dmin ;
102 minMax( rmin, rmax, direction.row( 0 ) ) ;
103 minMax( dmin, dmax, direction.row( 1 ) ) ;
104
105 Int npx = (Int)nx ;
106 Int npy = (Int)ny ;
107 String cellx( scellx ) ;
108 String celly( scelly ) ;
109 String center( scenter ) ;
110 setupGrid( npx, npy,
111 cellx, celly,
112 rmin, rmax,
113 dmin, dmax,
114 center ) ;
115}
116
117void STGrid::setFunc( string convType,
118 int convSupport )
119{
120 convType_ = String( convType ) ;
121 convType_.upcase() ;
122 userSupport_ = (Int)convSupport ;
123}
124
125#define NEED_UNDERSCORES
126#if defined(NEED_UNDERSCORES)
127#define ggridsd ggridsd_
128#endif
129extern "C" {
130 void ggridsd(Double*,
131 const Complex*,
132 Int*,
133 Int*,
134 Int*,
135 const Int*,
136 const Int*,
137 const Float*,
138 Int*,
139 Int*,
140 Complex*,
141 Float*,
142 Int*,
143 Int*,
144 Int *,
145 Int *,
146 Int*,
147 Int*,
148 Float*,
149 Int*,
150 Int*,
151 Double*);
152}
153void STGrid::grid()
154{
155 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
156
157 // retrieve data
158 Cube<Float> spectra ;
159 Matrix<Double> direction ;
160 Cube<uChar> flagtra ;
161 Matrix<uInt> rflag ;
162 Matrix<Float> weight ;
163 double t0, t1 ;
164 t0 = mathutil::gettimeofday_sec() ;
165 getData( spectra, direction, flagtra, rflag, weight ) ;
166 t1 = mathutil::gettimeofday_sec() ;
167 os << "getData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
168 IPosition sshape = spectra.shape() ;
169 //os << "spectra.shape()=" << spectra.shape() << LogIO::POST ;
170 //os << "max(spectra) = " << max(spectra) << LogIO::POST ;
171 //os << "weight = " << weight << LogIO::POST ;
172
173 // flagtra: uChar -> Int
174 // rflag: uInt -> Int
175 Cube<Int> flagI ;
176 Matrix<Int> rflagI ;
177 t0 = mathutil::gettimeofday_sec() ;
178 toInt( &flagtra, &flagI ) ;
179 toInt( &rflag, &rflagI ) ;
180 t1 = mathutil::gettimeofday_sec() ;
181 os << "toInt: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
182
183 // grid parameter
184 os << LogIO::DEBUGGING ;
185 os << "----------" << endl ;
186 os << "Grid parameter summary" << endl ;
187 os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
188 os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
189 os << " center = " << center_ << endl ;
190 os << "----------" << LogIO::POST ;
191 os << LogIO::NORMAL ;
192
193 // convolution kernel
194 Vector<Float> convFunc ;
195 t0 = mathutil::gettimeofday_sec() ;
196 setConvFunc( convFunc ) ;
197 t1 = mathutil::gettimeofday_sec() ;
198 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
199 //cout << "convSupport=" << convSupport_ << endl ;
200 //cout << "convFunc=" << convFunc << endl ;
201
202 // world -> pixel
203 Matrix<Double> xypos( direction.shape(), 0.0 ) ;
204 t0 = mathutil::gettimeofday_sec() ;
205 toPixel( direction, xypos ) ;
206 t1 = mathutil::gettimeofday_sec() ;
207 os << "toPixel: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
208
209 // call ggridsd
210 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
211 Double *xypos_p = xypos.getStorage( deletePos ) ;
212 Cube<Complex> dataC( spectra.shape(), 0.0 ) ;
213 setReal( dataC, spectra ) ;
214 const Complex *data_p = dataC.getStorage( deleteData ) ;
215 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
216 const Int *flag_p = flagI.getStorage( deleteFlag ) ;
217 const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
218 Float *conv_p = convFunc.getStorage( deleteConv ) ;
219 // Extend grid plane with convSupport_
220 Int gnx = nx_ ;
221 Int gny = ny_ ;
222// Int gnx = nx_+convSupport_*2 ;
223// Int gny = ny_+convSupport_*2 ;
224 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
225 Array<Complex> gdataArrC( gshape, 0.0 ) ;
226 Array<Float> gwgtArr( gshape, 0.0 ) ;
227 Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
228 Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
229 Int idopsf = 0 ;
230 Int *chanMap = new Int[nchan_] ;
231 {
232 Int *work_p = chanMap ;
233 for ( Int i = 0 ; i < nchan_ ; i++ ) {
234 *work_p = i ;
235 work_p++ ;
236 }
237 }
238 Int *polMap = new Int[npol_] ;
239 {
240 Int *work_p = polMap ;
241 for ( Int i = 0 ; i < npol_ ; i++ ) {
242 *work_p = i ;
243 work_p++ ;
244 }
245 }
246 Double *sumw_p = new Double[npol_*nchan_] ;
247 {
248 Double *work_p = sumw_p ;
249 for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
250 *work_p = 0.0 ;
251 work_p++ ;
252 }
253 }
254 t0 = mathutil::gettimeofday_sec() ;
255 Int irow = -1 ;
256 ggridsd( xypos_p,
257 data_p,
258 &npol_,
259 &nchan_,
260 &idopsf,
261 flag_p,
262 rflag_p,
263 wgt_p,
264 &nrow_,
265 &irow,
266 gdata_p,
267 wdata_p,
268 &gnx,
269 &gny,
270 &npol_,
271 &nchan_,
272 &convSupport_,
273 &convSampling_,
274 conv_p,
275 chanMap,
276 polMap,
277 sumw_p ) ;
278 t1 = mathutil::gettimeofday_sec() ;
279 os << "ggridsd: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
280 xypos.putStorage( xypos_p, deletePos ) ;
281 dataC.freeStorage( data_p, deleteData ) ;
282 weight.freeStorage( wgt_p, deleteWgt ) ;
283 flagI.freeStorage( flag_p, deleteFlag ) ;
284 rflagI.freeStorage( rflag_p, deleteFlagR ) ;
285 convFunc.putStorage( conv_p, deleteConv ) ;
286 delete polMap ;
287 delete chanMap ;
288 gdataArrC.putStorage( gdata_p, deleteDataG ) ;
289 gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
290 Array<Float> gdataArr = real( gdataArrC ) ;
291 t0 = mathutil::gettimeofday_sec() ;
292 data_.resize( gdataArr.shape() ) ;
293 data_ = 0.0 ;
294 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
295 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
296 for ( Int ip = 0 ; ip < npol_ ; ip++ ) {
297 for ( Int ic = 0 ; ic < nchan_ ; ic++ ) {
298 IPosition pos( 4, ix, iy, ip, ic ) ;
299 IPosition gpos( 4, ix, iy, ip, ic ) ;
300// IPosition gpos( 4, ix+convSupport_, iy+convSupport_, ip, ic ) ;
301 if ( gwgtArr( gpos ) > 0.0 )
302 data_( pos ) = gdataArr( gpos ) / gwgtArr( gpos ) ;
303 }
304 }
305 }
306 }
307 t1 = mathutil::gettimeofday_sec() ;
308 os << "set data: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
309 //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
310 delete sumw_p ;
311 //cout << "sumWeight = " << sumWeight << endl ;
312// os << "gdataArr = " << gdataArr << LogIO::POST ;
313// os << "gwgtArr = " << gwgtArr << LogIO::POST ;
314// os << "data_ " << data_ << LogIO::POST ;
315}
316
317void STGrid::setupGrid( Int &nx,
318 Int &ny,
319 String &cellx,
320 String &celly,
321 Double &xmin,
322 Double &xmax,
323 Double &ymin,
324 Double &ymax,
325 String &center )
326{
327 //cout << "nx=" << nx << ", ny=" << ny << endl ;
328
329 // center position
330 if ( center.size() == 0 ) {
331 center_(0) = 0.5 * ( xmin + xmax ) ;
332 center_(1) = 0.5 * ( ymin + ymax ) ;
333 }
334 else {
335 String::size_type pos0 = center.find( " " ) ;
336 if ( pos0 == String::npos ) {
337 throw AipsError( "bad string format in parameter center" ) ;
338 }
339 String::size_type pos1 = center.find( " ", pos0+1 ) ;
340 String typestr, xstr, ystr ;
341 if ( pos1 != String::npos ) {
342 typestr = center.substr( 0, pos0 ) ;
343 xstr = center.substr( pos0+1, pos1-pos0 ) ;
344 ystr = center.substr( pos1+1 ) ;
345 // todo: convert to J2000 (or direction ref for DIRECTION column)
346 }
347 else {
348 typestr = "J2000" ;
349 xstr = center.substr( 0, pos0 ) ;
350 ystr = center.substr( pos0+1 ) ;
351 }
352 QuantumHolder qh ;
353 String err ;
354 qh.fromString( err, xstr ) ;
355 Quantum<Double> xcen = qh.asQuantumDouble() ;
356 qh.fromString( err, ystr ) ;
357 Quantum<Double> ycen = qh.asQuantumDouble() ;
358 center_(0) = xcen.getValue( "rad" ) ;
359 center_(1) = ycen.getValue( "rad" ) ;
360 }
361
362
363 //Double wx = xmax - xmin ;
364 //Double wy = ymax - ymin ;
365 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
366 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
367 // take 10% margin
368 wx *= 1.10 ;
369 wy *= 1.10 ;
370 Quantum<Double> qcellx ;
371 Quantum<Double> qcelly ;
372 nx_ = nx ;
373 ny_ = ny ;
374 if ( nx < 0 && ny > 0 ) {
375 nx_ = ny ;
376 ny_ = ny ;
377 }
378 if ( ny < 0 && nx > 0 ) {
379 nx_ = nx ;
380 ny_ = nx ;
381 }
382 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
383 if ( cellx.size() != 0 && celly.size() != 0 ) {
384 readQuantity( qcellx, cellx ) ;
385 readQuantity( qcelly, celly ) ;
386 }
387 else if ( celly.size() != 0 ) {
388 cout << "Using celly to x-axis..." << endl ;
389 readQuantity( qcelly, celly ) ;
390 qcellx = qcelly ;
391 }
392 else if ( cellx.size() != 0 ) {
393 cout << "Using cellx to y-axis..." << endl ;
394 readQuantity( qcellx, cellx ) ;
395 qcelly = qcellx ;
396 }
397 else {
398 if ( nx_ < 0 ) {
399 cout << "No user preference in grid setting. Using default..." << endl ;
400 readQuantity( qcellx, "1.0arcmin" ) ;
401 qcelly = qcellx ;
402 }
403 else {
404 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
405 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
406 }
407 }
408 cellx_ = qcellx.getValue( "rad" ) ;
409 celly_ = qcelly.getValue( "rad" ) ;
410 if ( nx_ < 0 ) {
411 nx_ = Int( ceil( wx/cellx_ ) ) ;
412 ny_ = Int( ceil( wy/celly_ ) ) ;
413 }
414}
415
416void STGrid::selectData( Table &tab )
417{
418 Int ifno = ifno_ ;
419 Table taborg( infile_ ) ;
420 if ( ifno == -1 ) {
421 LogIO os( LogOrigin("STGrid","getData",WHERE) ) ;
422// os << LogIO::SEVERE
423// << "Please set IFNO before actual gridding"
424// << LogIO::EXCEPTION ;
425 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
426 ifno = ifnoCol( 0 ) ;
427 os << LogIO::WARN
428 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
429 }
430// tab = taborg( taborg.col("IFNO") == ifno ) ;
431 TableExprNode node ;
432 node = taborg.col("IFNO") == ifno ;
433 if ( scanlist_.size() > 0 ) {
434 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
435 }
436 tab = taborg( node ) ;
437 if ( tab.nrow() == 0 ) {
438 LogIO os( LogOrigin("STGrid","getData",WHERE) ) ;
439 os << LogIO::SEVERE
440 << "No corresponding rows for given selection: IFNO " << ifno
441 << " SCANNO " << scanlist_
442 << LogIO::EXCEPTION ;
443 }
444}
445
446void STGrid::getData( Cube<Float> &spectra,
447 Matrix<Double> &direction,
448 Cube<uChar> &flagtra,
449 Matrix<uInt> &rflag,
450 Matrix<Float> &weight )
451{
452 Table tab ;
453 selectData( tab ) ;
454 ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
455 Vector<uInt> pols = polnoCol.getColumn() ;
456 Vector<uInt> pollistOrg ;
457 uInt npolOrg = 0 ;
458 for ( uInt i = 0 ; i < pols.size() ; i++ ) {
459 if ( allNE( pollistOrg, pols[i] ) ) {
460 pollistOrg.resize( npolOrg+1, True ) ;
461 pollistOrg[npolOrg] = pols[i] ;
462 npolOrg++ ;
463 }
464 }
465 if ( pollist_.size() == 0 )
466 pollist_ = pollistOrg ;
467 else {
468 Vector<uInt> newlist ;
469 uInt newsize = 0 ;
470 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
471 if ( anyEQ( pols, pollist_[i] ) ) {
472 newlist.resize( newsize+1, True ) ;
473 newlist[newsize] = pollist_[i] ;
474 newsize++ ;
475 }
476 }
477 pollist_ = newlist ;
478 }
479 npol_ = pollist_.size() ;
480 ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
481 nchan_ = tmpCol( 0 ).nelements() ;
482 nrow_ = tab.nrow() / npolOrg ;
483// cout << "npol_ = " << npol_ << endl ;
484// cout << "nchan_ = " << nchan_ << endl ;
485// cout << "nrow_ = " << nrow_ << endl ;
486 spectra.resize( npol_, nchan_, nrow_ ) ;
487 flagtra.resize( npol_, nchan_, nrow_ ) ;
488 rflag.resize( npol_, nrow_ ) ;
489 Cube<Float> tsys( npol_, nchan_, nrow_ ) ;
490 Matrix<Double> tint( npol_, nrow_ ) ;
491 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
492 Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
493 ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
494 ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
495 ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
496 ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
497 ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
498 ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
499 Matrix<Float> tmpF = spectra.yzPlane( ipol ) ;
500 Matrix<uChar> tmpUC = flagtra.yzPlane( ipol ) ;
501 Vector<uInt> tmpUI = rflag.row( ipol ) ;
502 spectraCol.getColumn( tmpF ) ;
503 flagtraCol.getColumn( tmpUC ) ;
504 rflagCol.getColumn( tmpUI ) ;
505 if ( ipol == 0 )
506 directionCol.getColumn( direction ) ;
507 Matrix<Float> tmpF2 = tsysCol.getColumn() ;
508 Vector<Double> tmpD = tint.row( ipol ) ;
509 if ( tmpF2.shape()(0) == nchan_ ) {
510 tsys.yzPlane( ipol ) = tmpF2 ;
511 }
512 else {
513 tsys.yzPlane( ipol ) = tmpF2(0,0) ;
514 }
515 tintCol.getColumn( tmpD ) ;
516 }
517
518 getWeight( weight, tsys, tint ) ;
519}
520
521void STGrid::getWeight( Matrix<Float> &w,
522 Cube<Float> &tsys,
523 Matrix<Double> &tint )
524{
525 // resize
526 w.resize( nchan_, nrow_ ) ;
527
528 // set weight
529 w = 1.0 ;
530 Bool warn = False ;
531 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
532 // do nothing
533 }
534 else if ( wtype_.compare( "TINT" ) == 0 ) {
535 if ( npol_ > 1 ) warn = True ;
536 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
537 Float val = mean( tint.column( irow ) ) ;
538 w.column( irow ) = w.column( irow ) * val ;
539 }
540 }
541 else if ( wtype_.compare( "TSYS" ) == 0 ) {
542 if ( npol_ > 1 ) warn = True ;
543 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
544 Matrix<Float> arr = tsys.xyPlane( irow ) ;
545 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
546 Float val = mean( arr.column( ichan ) ) ;
547 w(ichan,irow) = w(ichan,irow) / ( val * val ) ;
548 }
549 }
550 }
551 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
552 if ( npol_ > 1 ) warn = True ;
553 for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
554 Float interval = mean( tint.column( irow ) ) ;
555 Matrix<Float> arr = tsys.xyPlane( irow ) ;
556 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
557 Float temp = mean( arr.column( ichan ) ) ;
558 w(ichan,irow) = w(ichan,irow) * interval / ( temp * temp ) ;
559 }
560 }
561 }
562 else {
563 LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
564 os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
565 }
566
567 if ( npol_ > 1 ) {
568 LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
569 os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
570 }
571}
572
573void STGrid::toInt( Array<uChar> *u, Array<Int> *v )
574{
575 uInt len = u->nelements() ;
576 Int *int_p = new Int[len] ;
577 Bool deleteIt ;
578 const uChar *data_p = u->getStorage( deleteIt ) ;
579 Int *i_p = int_p ;
580 const uChar *u_p = data_p ;
581 for ( uInt i = 0 ; i < len ; i++ ) {
582 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
583 i_p++ ;
584 u_p++ ;
585 }
586 u->freeStorage( data_p, deleteIt ) ;
587 v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
588}
589
590void STGrid::toInt( Array<uInt> *u, Array<Int> *v )
591{
592 uInt len = u->nelements() ;
593 Int *int_p = new Int[len] ;
594 Bool deleteIt ;
595 const uInt *data_p = u->getStorage( deleteIt ) ;
596 Int *i_p = int_p ;
597 const uInt *u_p = data_p ;
598 for ( uInt i = 0 ; i < len ; i++ ) {
599 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
600 i_p++ ;
601 u_p++ ;
602 }
603 u->freeStorage( data_p, deleteIt ) ;
604 v->takeStorage( u->shape(), int_p, TAKE_OVER ) ;
605}
606
607void STGrid::toPixel( Matrix<Double> &world, Matrix<Double> &pixel )
608{
609 // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
610 // grid plane to avoid unexpected behavior on grid edge
611 Vector<Double> pixc( 2 ) ;
612 pixc(0) = Double( nx_-1 ) * 0.5 ;
613 pixc(1) = Double( ny_-1 ) * 0.5 ;
614// pixc(0) = Double( nx_+2*convSupport_-1 ) * 0.5 ;
615// pixc(1) = Double( ny_+2*convSupport_-1 ) * 0.5 ;
616 uInt nrow = world.shape()[1] ;
617 Vector<Double> cell( 2 ) ;
618 cell(0) = cellx_ ;
619 cell(1) = celly_ ;
620 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
621 for ( uInt i = 0 ; i < 2 ; i++ ) {
622 pixel( i, irow ) = pixc(i) + ( world(i, irow) - center_(i) ) / cell(i) ;
623 }
624 }
625// String gridfile = "grid."+convType_+"."+String::toString(convSupport_)+".dat" ;
626// ofstream ofs( gridfile.c_str(), ios::out ) ;
627// ofs << "center " << center_(0) << " " << pixc(0)
628// << " " << center_(1) << " " << pixc(1) << endl ;
629// for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
630// ofs << irow ;
631// for ( uInt i = 0 ; i < 2 ; i++ ) {
632// ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
633// }
634// ofs << endl ;
635// }
636// ofs.close() ;
637}
638
639void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
640{
641 convFunc = 0.0 ;
642 for ( Int i = 0 ; i < convSize/2 ; i++ )
643 convFunc(i) = 1.0 ;
644}
645
646#define NEED_UNDERSCORES
647#if defined(NEED_UNDERSCORES)
648#define grdsf grdsf_
649#endif
650extern "C" {
651 void grdsf(Double*, Double*);
652}
653void STGrid::spheroidalFunc( Vector<Float> &convFunc )
654{
655 convFunc = 0.0 ;
656 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
657 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
658 Double val ;
659 grdsf( &nu, &val ) ;
660 convFunc(i) = ( 1.0 - nu * nu ) * val ;
661 }
662}
663
664void STGrid::gaussFunc( Vector<Float> &convFunc )
665{
666 convFunc = 0.0 ;
667 // HWHM of the Gaussian is convSupport_ / 4
668 // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
669 Int len = convSampling_ * convSupport_ ;
670 Double hwhm = len * 0.25 ;
671 for ( Int i = 0 ; i < len ; i++ ) {
672 Double val = Double(i) / hwhm ;
673 convFunc(i) = exp( -log(2)*val*val ) ;
674 }
675}
676
677void STGrid::pbFunc( Vector<Float> &convFunc )
678{
679 convFunc = 0.0 ;
680}
681
682void STGrid::setConvFunc( Vector<Float> &convFunc )
683{
684 convSupport_ = userSupport_ ;
685 if ( convType_ == "BOX" ) {
686 if ( convSupport_ < 0 )
687 convSupport_ = 0 ;
688 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
689 convFunc.resize( convSize ) ;
690 boxFunc( convFunc, convSize ) ;
691 }
692 else if ( convType_ == "SF" ) {
693 if ( convSupport_ < 0 )
694 convSupport_ = 3 ;
695 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
696 convFunc.resize( convSize ) ;
697 spheroidalFunc( convFunc ) ;
698 }
699 else if ( convType_ == "GAUSS" ) {
700 // to take into account Gaussian tail
701 if ( convSupport_ < 0 )
702 convSupport_ = 12 ; // 3 * 4
703 else {
704 convSupport_ = userSupport_ * 4 ;
705 }
706 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
707 convFunc.resize( convSize ) ;
708 gaussFunc( convFunc ) ;
709 }
710 else if ( convType_ == "PB" ) {
711 if ( convSupport_ < 0 )
712 convSupport_ = 0 ;
713 pbFunc( convFunc ) ;
714 }
715 else {
716 throw AipsError( "Unsupported convolution function" ) ;
717 }
718}
719
720string STGrid::saveData( string outfile )
721{
722 //Int polno = 0 ;
723 string outfile_ ;
724 if ( outfile.size() == 0 ) {
725 if ( infile_.lastchar() == '/' ) {
726 outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
727 }
728 else {
729 outfile_ = infile_ ;
730 }
731 outfile_ += ".grid" ;
732 }
733 else {
734 outfile_ = outfile ;
735 }
736 CountedPtr<Scantable> ref( new Scantable( infile_, Table::Memory ) ) ;
737 //cout << "ref->nchan()=" << ref->nchan() << endl ;
738 CountedPtr<Scantable> out( new Scantable( *ref, True ) ) ;
739 Table tab = out->table() ;
740 IPosition dshape = data_.shape() ;
741 Int nrow = nx_ * ny_ * npol_ ;
742 tab.rwKeywordSet().define( "nPol", npol_ ) ;
743 tab.addRow( nrow ) ;
744 Vector<Double> cpix( 2 ) ;
745 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
746 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
747 Vector<Double> dir( 2 ) ;
748 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
749 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
750 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
751 Int irow = 0 ;
752 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
753 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
754 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
755 IPosition start( 4, ix, iy, ipol, 0 ) ;
756 IPosition end( 4, ix, iy, ipol, nchan_-1 ) ;
757 IPosition inc( 4, 1, 1, 1, 1 ) ;
758 Vector<Float> sp = data_( start, end, inc ) ;
759 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
760 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
761 spectraCol.put( irow, sp ) ;
762 directionCol.put( irow, dir ) ;
763 polnoCol.put( irow, pollist_[ipol] ) ;
764 irow++ ;
765 }
766 }
767 }
768 //cout << "outfile_=" << outfile_ << endl ;
769 out->makePersistent( outfile_ ) ;
770
771 return outfile_ ;
772}
773
774}
Note: See TracBrowser for help on using the repository browser.