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