Line data    Source code

       1             : #include <synthesis/Utilities/PointingDirectionProjector.h>
       2             : 
       3             : #include <casacore/casa/BasicSL/Constants.h>
       4             : #include <casacore/casa/Arrays/Matrix.h>
       5             : #include <casacore/casa/Arrays/Vector.h>
       6             : #include <casacore/casa/Arrays/ArrayMath.h>
       7             : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
       8             : 
       9             : using namespace casacore;
      10             : using namespace casacore;
      11             : 
      12             : using namespace casacore;
      13             : namespace casa {
      14          18 : Projector::Projector() :
      15          18 :     user_defined_center_(false), user_defined_pcenter_(false) {
      16          18 : }
      17             : 
      18          18 : void Projector::setDirection(const Matrix<Double> &dir) {
      19          18 :   dir_.reference(dir.copy());
      20          36 :   Vector<Double> ra(dir_.row(0));
      21          18 :   rotateRA(ra);
      22          18 : }
      23             : 
      24           9 : void Projector::setReferenceCoordinate(Double const lat, Double const lon) {
      25           9 :   cenx_user_ = lat;
      26           9 :   ceny_user_ = lon;
      27           9 :   user_defined_center_ = true;
      28           9 : }
      29           9 : void Projector::setReferencePixel(Double const refx, Double const refy) {
      30           9 :   pcenx_user_ = refx;
      31           9 :   pceny_user_ = refy;
      32           9 :   user_defined_pcenter_ = true;
      33           9 : }
      34             : 
      35           0 : void Projector::unsetReferenceCoordinate() {
      36           0 :   user_defined_center_ = false;
      37           0 : }
      38           0 : void Projector::unsetReferencePixel() {
      39           0 :   user_defined_pcenter_ = false;
      40           0 : }
      41             : 
      42          18 : void Projector::rotateRA(Vector<Double> &v) {
      43          18 :   uInt len = v.nelements();
      44          36 :   Vector<Double> work(len);
      45             : 
      46       48207 :   for (uInt i = 0; i < len; i++) {
      47       48189 :     work[i] = fmod(v[i], C::_2pi);
      48       48189 :     if (work[i] < 0.0) {
      49       47519 :       work[i] += C::_2pi;
      50             :     }
      51             :   }
      52             : 
      53          36 :   Vector<uInt> quad(len);
      54          36 :   Vector<uInt> nquad(4, 0);
      55       48207 :   for (uInt i = 0; i < len; i++) {
      56       48189 :     uInt q = uInt(work[i] / C::pi_2);
      57       48189 :     nquad[q]++;
      58       48189 :     quad[i] = q;
      59             :   }
      60             : 
      61          36 :   Vector<Bool> rot(4, false);
      62          18 :   if (nquad[0] > 0 && nquad[3] > 0 && (nquad[1] == 0 || nquad[2] == 0)) {
      63             :     //cout << "need rotation" << endl ;
      64           5 :     rot[3] = true;
      65           5 :     rot[2] = (nquad[1] == 0 && nquad[2] > 0);
      66             :   }
      67             : 
      68       48207 :   for (uInt i = 0; i < len; i++) {
      69       48189 :     if (rot[quad[i]]) {
      70         610 :       v[i] = work[i] - C::_2pi;
      71             :     } else {
      72       47579 :       v[i] = work[i];
      73             :     }
      74             :   }
      75          18 : }
      76             : 
      77          18 : OrthographicProjector::~OrthographicProjector() {
      78             :   // Do nothing
      79          18 : }
      80             : 
      81          18 : OrthographicProjector::OrthographicProjector(Float pixel_scale) :
      82          18 :     Projector(), pixel_scale_(pixel_scale), p_center_(2, 0.0), p_size_(2, 0.0)
      83             : 
      84             : {
      85          18 : }
      86             : 
      87          18 : const Matrix<Double>& OrthographicProjector::project() {
      88          18 :   scale_and_center();
      89             :   // using DirectionCoordinate
      90          36 :   Matrix<Double> identity(2, 2, Double(0.0));
      91          18 :   identity.diagonal() = 1.0;
      92           0 :   DirectionCoordinate coord(MDirection::J2000, Projection(Projection::SIN),
      93          36 :       cenx_, ceny_, dx_, dy_, identity, pcenx_, pceny_);
      94             : 
      95          18 :   pdir_ = Matrix<Double>(dir_.shape());
      96          18 :   double* pdir_p = pdir_.data();
      97          18 :   size_t len = dir_.ncolumn();
      98             :   Bool b;
      99          18 :   Double *dir_p = dir_.getStorage(b);
     100          18 :   Double *wdir_p = dir_p;
     101          36 :   Vector<Double> world;
     102          36 :   Vector<Double> pixel;
     103          18 :   IPosition vshape(1, 2);
     104       48207 :   for (uInt i = 0; i < len; i++) {
     105       48189 :     world.takeStorage(vshape, wdir_p, SHARE);
     106       48189 :     pixel.takeStorage(vshape, pdir_p, SHARE);
     107       48189 :     coord.toPixel(pixel, world);
     108       48189 :     pdir_p += 2;
     109       48189 :     wdir_p += 2;
     110             :   }
     111          18 :   dir_.putStorage(dir_p, b);
     112          36 :   return pdir_;
     113             : }
     114             : 
     115          18 : void OrthographicProjector::scale_and_center() {
     116          18 :   os_.origin(LogOrigin("OrthographicProjector", "scale_and_center", WHERE));
     117             : 
     118             :   Double xmax, xmin, ymax, ymin;
     119          18 :   minMax( xmin, xmax, dir_.row( 0 ) );
     120          18 :   minMax( ymin, ymax, dir_.row( 1 ) );
     121          18 :   Double wx = ( xmax - xmin ) * 1.1;
     122          18 :   Double wy = ( ymax - ymin ) * 1.1;
     123             : 
     124          18 :   if (isReferenceCoordinateSet()) {
     125           9 :     getUserDefinedReferenceCoordinate(cenx_, ceny_);
     126             :   } else {
     127           9 :     cenx_ = 0.5 * ( xmin + xmax );
     128           9 :     ceny_ = 0.5 * ( ymin + ymax );
     129             :   }
     130          18 :   Double decCorr = cos( ceny_ );
     131             : 
     132             :   // Renaud: uInt len = time_.nelements() ;
     133          18 :       uInt len = dir_.ncolumn();
     134          36 :       Matrix<Double> dd = dir_.copy();
     135       48189 :       for ( uInt i = len-1; i > 0; i-- ) {
     136             :         //dd(0,i) = ( dd(0,i) - dd(0,i-1) ) * decCorr ;
     137       48171 :         dd(0,i) = ( dd(0,i) - dd(0,i-1) ) * cos( 0.5*(dd(1,i-1)+dd(1,i)) );
     138       48171 :         dd(1,i) = dd(1,i) - dd(1,i-1);
     139             :       }
     140          36 :       Vector<Double> dr( len-1 );
     141             :       Bool b;
     142          18 :       const Double *dir_p = dd.getStorage( b );
     143          18 :       const Double *x_p = dir_p + 2;
     144          18 :       const Double *y_p = dir_p + 3;
     145       48189 :       for ( uInt i = 0; i < len-1; i++ ) {
     146       48171 :         dr[i] = sqrt( (*x_p) * (*x_p) + (*y_p) * (*y_p) );
     147       48171 :         x_p += 2;
     148       48171 :         y_p += 2;
     149             :       }
     150          18 :       dir_.freeStorage( dir_p, b );
     151          18 :       Double med = median( dr, false, true, true );
     152          18 :       dy_ = med * pixel_scale_;
     153          18 :       dx_ = dy_ / decCorr;
     154             : 
     155          18 :       Double nxTemp = ceil(wx / dx_);
     156          18 :       Double nyTemp = ceil(wy / dy_);
     157             : 
     158          18 :       os_ << LogIO::DEBUGGING
     159             :       << "len = " << len
     160             :       << "range x = (" << xmin << "," << xmax << ")" << endl
     161             :       << "range y = (" << ymin << "," << ymax << ")" << endl
     162             :       << "direction center = (" << cenx_ << "," << ceny_ << ")" << endl
     163             :       << "declination correction: cos(dir_center.y)=" << decCorr << endl
     164             :       << "median separation between pointings: " << med << endl
     165             :       << "dx=" << dx_ << ", dy=" << dy_ << endl
     166             :       << "wx=" << wx << ", wy=" << wy << endl
     167          18 :       << "nxTemp=" << nxTemp << ", nyTemp=" << nyTemp << LogIO::POST;
     168             : 
     169          18 :       if (nxTemp > (Double)UINT_MAX || nyTemp > (Double)UINT_MAX) {
     170           0 :         throw AipsError("Error in setup: Too large number of pixels.");
     171             :       }
     172          18 :       nx_ = uInt( nxTemp );
     173          18 :       ny_ = uInt( nyTemp );
     174             : 
     175             :       // Renaud debug
     176          18 :       p_size_[0] = nxTemp;
     177          18 :       p_size_[1] = nyTemp;
     178             : 
     179          18 :       if (isReferencePixelSet()) {
     180           9 :         getUserDefinedReferencePixel(pcenx_, pceny_);
     181             :       } else {
     182           9 :         pcenx_ = 0.5 * Double( nx_ - 1 );
     183           9 :         pceny_ = 0.5 * Double( ny_ - 1 );
     184             :       }
     185             : 
     186             :       // Renaud debug
     187          18 :       p_center_[0] = pcenx_;
     188          18 :       p_center_[1] = pceny_;
     189             : 
     190          18 :       os_ << LogIO::DEBUGGING
     191             :       << "pixel center = (" << pcenx_ << "," << pceny_ << ")" << endl
     192             :       << "nx=" << nx_ << ", ny=" << ny_
     193          18 :       << "n_pointings=" << len << " must be < n_pixels=" << nx_ * ny_ << LogIO::POST ;
     194          18 : }
     195             : 
     196             : }