source: trunk/src/STGrid.cpp@ 2356

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

New Development: Yes

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...

Defined new class STGrid and its python interface.
The STGrid class performs single dish gridding.
At the moment, only single polarization data is supported.


File size: 15.1 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/Inputs/Input.h>
10#include <casa/Quanta/Quantum.h>
11#include <casa/Quanta/QuantumHolder.h>
12#include <casa/Utilities/CountedPtr.h>
13
14#include <tables/Tables/Table.h>
15#include <tables/Tables/ScalarColumn.h>
16#include <tables/Tables/ArrayColumn.h>
17
18#include <measures/Measures/MDirection.h>
19
20#include <Scantable.h>
21#include "STGrid.h"
22
23using namespace std ;
24using namespace casa ;
25using namespace asap ;
26
27namespace asap {
28
29// constructor
30STGrid::STGrid()
31{
32 init() ;
33}
34
35STGrid::STGrid( const string infile )
36{
37 init() ;
38
39 setFileIn( infile ) ;
40}
41
42void STGrid::init()
43{
44 nx_ = -1 ;
45 ny_ = -1 ;
46 cellx_ = 0.0 ;
47 celly_ = 0.0 ;
48 center_ = Vector<Double> ( 2, 0.0 ) ;
49 convType_ = "BOX" ;
50 convSupport_ = -1 ;
51 userSupport_ = -1 ;
52 convSampling_ = 100 ;
53}
54
55void STGrid::setFileIn( const string infile )
56{
57 String name( infile ) ;
58 if ( infile_.compare( name ) != 0 ) {
59 infile_ = String( infile ) ;
60 tab_ = Table( infile_ ) ;
61 }
62}
63
64void STGrid::defineImage( int nx,
65 int ny,
66 string scellx,
67 string scelly,
68 string scenter )
69{
70 ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;
71 Matrix<Double> direction = dirCol.getColumn() ;
72 Double rmax, rmin, dmax, dmin ;
73 minMax( rmin, rmax, direction.row( 0 ) ) ;
74 minMax( dmin, dmax, direction.row( 1 ) ) ;
75
76 Int npx = (Int)nx ;
77 Int npy = (Int)ny ;
78 String cellx( scellx ) ;
79 String celly( scelly ) ;
80 String center( scenter ) ;
81 setupGrid( npx, npy,
82 cellx, celly,
83 rmin, rmax,
84 dmin, dmax,
85 center ) ;
86}
87
88void STGrid::setOption( string convType,
89 int convSupport )
90{
91 convType_ = String( convType ) ;
92 convType_.upcase() ;
93 userSupport_ = (Int)convSupport ;
94}
95
96#define NEED_UNDERSCORES
97#if defined(NEED_UNDERSCORES)
98#define ggridsd ggridsd_
99#endif
100extern "C" {
101 void ggridsd(Double*,
102 const Complex*,
103 Int*,
104 Int*,
105 Int*,
106 const Int*,
107 const Int*,
108 const Float*,
109 Int*,
110 Int*,
111 Complex*,
112 Float*,
113 Int*,
114 Int*,
115 Int *,
116 Int *,
117 Int*,
118 Int*,
119 Float*,
120 Int*,
121 Int*,
122 Double*);
123}
124void STGrid::grid()
125{
126
127 // retrieve data
128 Matrix<Float> spectra ;
129 Matrix<Double> direction ;
130 Matrix<uChar> flagtra ;
131 Vector<uInt> rflag ;
132 getData( infile_, spectra, direction, flagtra, rflag ) ;
133 IPosition sshape = spectra.shape() ;
134 Int nchan = sshape[0] ;
135 Int nrow = sshape[1] ;
136 //cout << "data.shape()=" << data.shape() << endl ;
137
138 // flagtra: uChar -> Int
139 // rflag: uInt -> Int
140 Matrix<Int> flagI ;
141 Vector<Int> rflagI ;
142 toInt( flagtra, flagI ) ;
143 toInt( rflag, rflagI ) ;
144
145 //cout << "flagI.shape() = " << flagI.shape() << endl ;
146 //cout << "rflagI.shape() = " << rflagI.shape() << endl ;
147
148 // grid parameter
149 cout << "----------" << endl ;
150 cout << "Grid parameter summary" << endl ;
151 cout << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
152 cout << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
153 cout << " center = " << center_ << endl ;
154 cout << "----------" << endl ;
155
156 // world -> pixel
157 Matrix<Double> xypos( direction.shape(), 0.0 ) ;
158 toPixel( direction, xypos ) ;
159 //cout << "max(xypos.row(0))=" << max(xypos.row(0)) << endl ;
160 //cout << "min(xypos.row(0))=" << min(xypos.row(0)) << endl ;
161 //cout << "max(xypos.row(1))=" << max(xypos.row(1)) << endl ;
162 //cout << "min(xypos.row(1))=" << min(xypos.row(1)) << endl ;
163// for ( Int irow = 0 ; irow < nrow ; irow++ ) {
164// cout << irow << ": xypos=" << xypos.column( irow )
165// << " data = " << spectra.column( irow ) << endl ;
166// }
167
168 // convolution kernel
169 Vector<Float> convFunc ;
170 setConvFunc( convFunc ) ;
171 //cout << "convSupport=" << convSupport_ << endl ;
172 //cout << "convFunc=" << convFunc << endl ;
173
174 // weighting factor
175 Matrix<Float> weight( IPosition( 2, nchan, nrow ), 1.0 ) ;
176
177 // call ggridsd
178 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
179 Double *xypos_p = xypos.getStorage( deletePos ) ;
180 Matrix<Complex> dataC( spectra.shape(), 0.0 ) ;
181 setReal( dataC, spectra ) ;
182 const Complex *data_p = dataC.getStorage( deleteData ) ;
183 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
184 const Int *flag_p = flagI.getStorage( deleteFlag ) ;
185 const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
186 Float *conv_p = convFunc.getStorage( deleteConv ) ;
187 Int npol = 1 ;
188 IPosition gshape( 4, nx_, ny_, npol, nchan ) ;
189 Array<Complex> gdataArrC( gshape, 0.0 ) ;
190 Array<Float> gwgtArr( gshape, 0.0 ) ;
191 Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
192 Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
193 Int idopsf = 0 ;
194 Int irow = -1 ;
195 Int *chanMap = new Int[nchan] ;
196 {
197 Int *work_p = chanMap ;
198 for ( Int i = 0 ; i < nchan ; i++ ) {
199 *work_p = i ;
200 work_p++ ;
201 }
202 }
203 Int *polMap = new Int[npol] ;
204 {
205 Int *work_p = polMap ;
206 for ( Int i = 0 ; i < npol ; i++ ) {
207 *work_p = i ;
208 work_p++ ;
209 }
210 }
211 Double *sumw_p = new Double[npol*nchan] ;
212 {
213 Double *work_p = sumw_p ;
214 for ( Int i = 0 ; i < npol*nchan ; i++ ) {
215 *work_p = 0.0 ;
216 work_p++ ;
217 }
218 }
219 ggridsd( xypos_p,
220 data_p,
221 &npol,
222 &nchan,
223 &idopsf,
224 flag_p,
225 rflag_p,
226 wgt_p,
227 &nrow,
228 &irow,
229 gdata_p,
230 wdata_p,
231 &nx_,
232 &ny_,
233 &npol,
234 &nchan,
235 &convSupport_,
236 &convSampling_,
237 conv_p,
238 chanMap,
239 polMap,
240 sumw_p ) ;
241 xypos.putStorage( xypos_p, deletePos ) ;
242 dataC.freeStorage( data_p, deleteData ) ;
243 weight.freeStorage( wgt_p, deleteWgt ) ;
244 flagI.freeStorage( flag_p, deleteFlag ) ;
245 rflagI.freeStorage( rflag_p, deleteFlagR ) ;
246 convFunc.putStorage( conv_p, deleteConv ) ;
247 delete polMap ;
248 delete chanMap ;
249 gdataArrC.putStorage( gdata_p, deleteDataG ) ;
250 gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
251 Array<Float> gdataArr = real( gdataArrC ) ;
252 //Array<Float> gdataArrN( gdataArr.shape(), 0.0 ) ;
253 data_.resize( gdataArr.shape() ) ;
254 data_ = 0.0 ;
255 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
256 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
257 for ( Int ip = 0 ; ip < npol ; ip++ ) {
258 for ( Int ic = 0 ; ic < nchan ; ic++ ) {
259 IPosition pos( 4, ix, iy, ip, ic ) ;
260 if ( gwgtArr( pos ) > 0.0 )
261 //gdataArrN( pos ) = gdataArr( pos ) / gwgtArr( pos ) ;
262 data_( pos ) = gdataArr( pos ) / gwgtArr( pos ) ;
263 }
264 }
265 }
266 }
267 Matrix<Double> sumWeight( IPosition( 2, npol, nchan ), sumw_p, TAKE_OVER ) ;
268 //cout << "sumWeight = " << sumWeight << endl ;
269 //cout << "gdataArr = " << gdataArr << endl ;
270 //cout << "gwgtArr = " << gwgtArr << endl ;
271 //cout << "gdataArr/gwgtArr = " << gdataArrN << endl ;
272}
273
274void STGrid::setupGrid( Int &nx,
275 Int &ny,
276 String &cellx,
277 String &celly,
278 Double &xmin,
279 Double &xmax,
280 Double &ymin,
281 Double &ymax,
282 String &center )
283{
284 //cout << "nx=" << nx << ", ny=" << ny << endl ;
285 Double wx = xmax - xmin ;
286 Double wy = ymax - ymin ;
287 // take some margin
288 wx *= 1.10 ;
289 wy *= 1.10 ;
290 Quantum<Double> qcellx ;
291 Quantum<Double> qcelly ;
292 nx_ = nx ;
293 ny_ = ny ;
294 if ( nx < 0 && ny > 0 ) {
295 nx_ = ny ;
296 ny_ = ny ;
297 }
298 if ( ny < 0 && nx > 0 ) {
299 nx_ = nx ;
300 ny_ = nx ;
301 }
302 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
303 if ( cellx.size() != 0 && celly.size() != 0 ) {
304 readQuantity( qcellx, cellx ) ;
305 readQuantity( qcelly, celly ) ;
306 }
307 else if ( celly.size() != 0 ) {
308 cout << "Using celly to x-axis..." << endl ;
309 readQuantity( qcelly, celly ) ;
310 qcellx = qcelly ;
311 }
312 else if ( cellx.size() != 0 ) {
313 cout << "Using cellx to y-axis..." << endl ;
314 readQuantity( qcellx, cellx ) ;
315 qcelly = qcellx ;
316 }
317 else {
318 if ( nx_ < 0 ) {
319 cout << "No user preference in grid setting. Using default..." << endl ;
320 readQuantity( qcellx, "1.0arcmin" ) ;
321 qcelly = qcellx ;
322 }
323 else {
324 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
325 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
326 }
327 }
328 cellx_ = qcellx.getValue( "rad" ) ;
329 celly_ = qcelly.getValue( "rad" ) ;
330 if ( nx_ < 0 ) {
331 nx_ = Int( ceil( wx/cellx_ ) ) ;
332 ny_ = Int( ceil( wy/celly_ ) ) ;
333 }
334
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
368void STGrid::getData( String &infile,
369 Matrix<Float> &spectra,
370 Matrix<Double> &direction,
371 Matrix<uChar> &flagtra,
372 Vector<uInt> &rflag )
373{
374 Table tab( infile ) ;
375 ROArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
376 ROArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
377 ROArrayColumn<uChar> flagtraCol( tab, "FLAGTRA" ) ;
378 ROScalarColumn<uInt> rflagCol( tab, "FLAGROW" ) ;
379 spectraCol.getColumn( spectra ) ;
380 directionCol.getColumn( direction ) ;
381 flagtraCol.getColumn( flagtra ) ;
382 rflagCol.getColumn( rflag ) ;
383}
384
385void STGrid::toInt( Matrix<uChar> &u, Matrix<Int> &v )
386{
387 uInt len = u.nelements() ;
388 Int *int_p = new Int[len] ;
389 Bool deleteIt ;
390 const uChar *data_p = u.getStorage( deleteIt ) ;
391 Int *i_p = int_p ;
392 const uChar *u_p = data_p ;
393 for ( uInt i = 0 ; i < len ; i++ ) {
394 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
395 i_p++ ;
396 u_p++ ;
397 }
398 u.freeStorage( data_p, deleteIt ) ;
399 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
400}
401
402void STGrid::toInt( Vector<uInt> &u, Vector<Int> &v )
403{
404 uInt len = u.nelements() ;
405 Int *int_p = new Int[len] ;
406 Bool deleteIt ;
407 const uInt *data_p = u.getStorage( deleteIt ) ;
408 Int *i_p = int_p ;
409 const uInt *u_p = data_p ;
410 for ( uInt i = 0 ; i < len ; i++ ) {
411 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
412 i_p++ ;
413 u_p++ ;
414 }
415 u.freeStorage( data_p, deleteIt ) ;
416 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
417}
418
419void STGrid::toPixel( Matrix<Double> &world, Matrix<Double> &pixel )
420{
421 Vector<Double> pixc( 2 ) ;
422 pixc(0) = Double( nx_-1 ) * 0.5 ;
423 pixc(1) = Double( ny_-1 ) * 0.5 ;
424 uInt nrow = world.shape()[1] ;
425 Vector<Double> cell( 2 ) ;
426 cell(0) = cellx_ ;
427 cell(1) = celly_ ;
428 //ofstream ofs( "grid.dat", ios::out ) ;
429 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
430 //ofs << irow ;
431 for ( uInt i = 0 ; i < 2 ; i++ ) {
432 pixel( i, irow ) = pixc(i) + ( world(i, irow) - center_(i) ) / cell(i) ;
433 //ofs << " " << world(i, irow) << " " << pixel(i, irow) ;
434 }
435 //ofs << endl ;
436 }
437 //ofs.close() ;
438}
439
440void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
441{
442 convFunc = 0.0 ;
443 for ( Int i = 0 ; i < convSize/2 ; i++ )
444 convFunc(i) = 1.0 ;
445}
446
447#define NEED_UNDERSCORES
448#if defined(NEED_UNDERSCORES)
449#define grdsf grdsf_
450#endif
451extern "C" {
452 void grdsf(Double*, Double*);
453}
454void STGrid::spheroidalFunc( Vector<Float> &convFunc )
455{
456 convFunc = 0.0 ;
457 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
458 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
459 Double val ;
460 grdsf( &nu, &val ) ;
461 convFunc(i) = ( 1.0 - nu * nu ) * val ;
462 }
463}
464
465void STGrid::gaussFunc( Vector<Float> &convFunc )
466{
467 convFunc = 0.0 ;
468 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
469 Double hwhm = convSampling_ * convSupport_ * 0.25 ;
470 Double val = Double(i) / hwhm ;
471 convFunc(i) = exp( -log(2)*val*val ) ;
472 }
473}
474
475void STGrid::pbFunc( Vector<Float> &convFunc )
476{
477 convFunc = 0.0 ;
478}
479
480void STGrid::setConvFunc( Vector<Float> &convFunc )
481{
482 convSupport_ = userSupport_ ;
483 if ( convType_ == "BOX" ) {
484 if ( convSupport_ < 0 )
485 convSupport_ = 0 ;
486 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
487 convFunc.resize( convSize ) ;
488 boxFunc( convFunc, convSize ) ;
489 }
490 else if ( convType_ == "SF" ) {
491 if ( convSupport_ < 0 )
492 convSupport_ = 3 ;
493 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
494 convFunc.resize( convSize ) ;
495 spheroidalFunc( convFunc ) ;
496 }
497 else if ( convType_ == "GAUSS" ) {
498 if ( convSupport_ < 0 )
499 convSupport_ = 3 ;
500 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
501 convFunc.resize( convSize ) ;
502 gaussFunc( convFunc ) ;
503 }
504 else if ( convType_ == "PB" )
505 pbFunc( convFunc ) ;
506 else {
507 throw AipsError( "Unsupported convolution function" ) ;
508 }
509}
510
511string STGrid::saveData( string outfile )
512{
513 Int polno = 0 ;
514 string outfile_ ;
515 if ( outfile.size() == 0 ) {
516 if ( infile_.lastchar() == '/' ) {
517 outfile_ = infile_.substr( 0, infile_.size()-1 ) ;
518 }
519 else {
520 outfile_ = infile_ ;
521 }
522 outfile_ += ".grid" ;
523 }
524 else {
525 outfile_ = outfile ;
526 }
527 CountedPtr<Scantable> ref( new Scantable( infile_, Table::Memory ) ) ;
528 //cout << "ref->nchan()=" << ref->nchan() << endl ;
529 CountedPtr<Scantable> out( new Scantable( *ref, True ) ) ;
530 Table tab = out->table() ;
531 Int nrow = nx_ * ny_ ;
532 tab.addRow( nrow ) ;
533 IPosition dshape = data_.shape() ;
534 Int npol = dshape[2] ;
535 Int nchan = dshape[3] ;
536 Vector<Double> cpix( 2 ) ;
537 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
538 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
539 Vector<Double> dir( 2 ) ;
540 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
541 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
542 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
543 Int irow = 0 ;
544 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
545 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
546 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
547 IPosition start( 4, ix, iy, ipol, 0 ) ;
548 IPosition end( 4, ix, iy, ipol, nchan-1 ) ;
549 IPosition inc( 4, 1, 1, 1, 1 ) ;
550 Vector<Float> sp = data_( start, end, inc ) ;
551 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
552 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
553 spectraCol.put( irow, sp ) ;
554 directionCol.put( irow, dir ) ;
555 polnoCol.put( irow, polno ) ;
556 irow++ ;
557 }
558 }
559 }
560 //cout << "outfile_=" << outfile_ << endl ;
561 out->makePersistent( outfile_ ) ;
562
563 return outfile_ ;
564}
565
566}
Note: See TracBrowser for help on using the repository browser.