simplify wcs code

4 files changed, 50 insertions, 13 deletions

diff --git a/tksao/frame/basecommand.C b/tksao/frame/basecommand.C
index 49e65c0..4871fba 100644
--- a/tksao/frame/basecommand.C
+++ b/tksao/frame/basecommand.C
@@ -562,7 +562,7 @@ void Base::crop3dCmd(double z0, double z1, Coord::CoordSystem sys,
     return;
 
   // center in IMAGE
-  Vector3d cc = Vector3d(ptr->center(),1) * Translate3d(-.5,-.5,-.5);
+  Vector3d cc = Vector3d(ptr->getWCSCRPIX(sys),1) * Translate3d(-.5,-.5,-.5);
   Vector3d wcc = ptr->mapFromRef(cc,sys,sky);
   Vector3d min = ptr->mapToRef(Vector3d(wcc[0],wcc[1],z0),sys,sky);
   Vector3d max = ptr->mapToRef(Vector3d(wcc[0],wcc[1],z1),sys,sky);
@@ -1420,10 +1420,10 @@ void Base::getCrop3dCmd(Coord::CoordSystem sys, Coord::SkyFrame sky)
 
   // need to move from edge to center of pixel
   Vector3d rmin =
-    Vector3d(ptr->center(),zparams->zmin) * Translate3d(.5,.5,.5);
+    Vector3d(ptr->getWCSCRPIX(sys),zparams->zmin) * Translate3d(.5,.5,.5);
   // need to move from edge to center of pixel
   Vector3d rmax =
-    Vector3d(ptr->center(),zparams->zmax) * Translate3d(-.5,-.5,-.5);
+    Vector3d(ptr->getWCSCRPIX(sys),zparams->zmax) * Translate3d(-.5,-.5,-.5);
   Vector3d min = ptr->mapFromRef(rmin,sys,sky);
   Vector3d max = ptr->mapFromRef(rmax,sys,sky);
 
@@ -1778,7 +1778,7 @@ void Base::getFitsSliceFromImageCmd(int ss, Coord::CoordSystem sys,
 {
   if (currentContext->cfits) {
     FitsImage* ptr = currentContext->fits;
-    Vector3d dd = Vector3d(ptr->center(),ss) * Translate3d(-.5,-.5,-.5);
+    Vector3d dd = Vector3d(ptr->getWCSCRPIX(sys),ss) * Translate3d(-.5,-.5,-.5);
     Vector3d out = ptr->mapFromRef(dd,sys,sky);
     printDouble(out[2]);
   }
@@ -1791,10 +1791,10 @@ void Base::getFitsSliceToImageCmd(double dd, Coord::CoordSystem sys,
 {
   if (currentContext->cfits) {
     FitsImage* ptr = currentContext->fits;
-    Vector3d cc = Vector3d(ptr->center(),1) * Translate3d(-.5,-.5,-.5);
+    Vector3d cc = Vector3d(ptr->getWCSCRPIX(sys),1) * Translate3d(-.5,-.5,-.5);
     Vector3d wcc = ptr->mapFromRef(cc,sys,sky);
-    Vector3d out = ptr->mapToRef(Vector3d(wcc[0],wcc[1],dd),sys,sky)
-      * Translate3d(.5,.5,.5);
+    Vector3d oo = ptr->mapToRef(Vector3d(wcc[0],wcc[1],dd),sys,sky);
+    Vector3d out = oo * Translate3d(.5,.5,.5);
     printInteger(out[2]);
   }
   else
@@ -2866,14 +2866,13 @@ void Base::sliceCmd(double dd, Coord::CoordSystem sys, Coord::SkyFrame sky)
     return;
   
   FitsImage* ptr = currentContext->fits;
-  Vector3d cc = Vector3d(ptr->center(),1) * Translate3d(-.5,-.5,-.5);
+  Vector3d cc = Vector3d(ptr->getWCSCRPIX(sys),1) * Translate3d(-.5,-.5,-.5);
   Vector3d wcc = ptr->mapFromRef(cc,sys,sky);
-  Vector3d out = ptr->mapToRef(Vector3d(wcc[0],wcc[1],dd),sys,sky)
-    * Translate3d(.5,.5,.5);
-  int ss = out[2];
+  Vector3d oo = ptr->mapToRef(Vector3d(wcc[0],wcc[1],dd),sys,sky);
+  Vector3d out = oo * Translate3d(.5,.5,.5);
   
   // IMAGE (ranges 1-n)
-  setSlice(2,ss);
+  setSlice(2,out[2]);
   updateMagnifier();
 }
 
diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 0b12fb4..bf4c8ba 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -89,6 +89,7 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
   wcsCel_ =NULL;
   wcsCelLon_ =NULL;
   wcsCelLat_ =NULL;
+  wcsCRPIX_ =NULL;
   wcsSize_ =NULL;
 
   wcsHPX_ =0;
@@ -168,6 +169,8 @@ FitsImage::~FitsImage()
     if (wcsCelLat_)
       delete [] wcsCelLat_;
 
+    if (wcsCRPIX_)
+      delete [] wcsCRPIX_;
     if (wcsSize_)
       delete [] wcsSize_;
 
@@ -1065,6 +1068,9 @@ void FitsImage::initWCS(FitsHead* hd)
       delete [] wcsCelLat_;
     wcsCelLat_ =NULL;
 
+    if (wcsCRPIX_)
+      delete [] wcsCRPIX_;
+    wcsCRPIX_ =NULL;
     if (wcsSize_)
       delete [] wcsSize_;
     wcsSize_ =NULL;
@@ -1092,6 +1098,7 @@ void FitsImage::initWCS(FitsHead* hd)
 	  wcsCel_ = ptr->wcsCel_;
 	  wcsCelLon_ = ptr->wcsCelLon_;
 	  wcsCelLat_ = ptr->wcsCelLat_;
+	  wcsCRPIX_ = ptr->wcsCRPIX_;
 	  wcsSize_ = ptr->wcsSize_;
 
 	  wcsHPX_ = ptr->wcsHPX_;
@@ -1129,6 +1136,7 @@ void FitsImage::initWCS(FitsHead* hd)
 
   wcsInit(hasWCSAST);
   wcsCelInit(hasWCSAST);
+  wcsCRPIXInit(hd);
   wcsHPXInit();
 
   // init wcsState
@@ -2583,6 +2591,14 @@ void FitsImage::updatePS(Matrix3d ps)
 
 // WCS
 
+Vector FitsImage::getWCSCRPIX(Coord::CoordSystem sys)
+{
+  if (!wcsCRPIX_ || sys<Coord::WCS)
+    return Vector();
+  else
+    return wcsCRPIX_[sys-Coord::WCS];
+}
+
 double FitsImage::getWCSSize(Coord::CoordSystem sys)
 {
   if (!wcsSize_ || sys<Coord::WCS)
@@ -3087,6 +3103,24 @@ void FitsImage::wcsHPXInit()
   }
 }
 
+void FitsImage::wcsCRPIXInit(FitsHead* hd)
+{
+  // init wcsCel_ array
+  if (wcsCRPIX_)
+    delete [] wcsCRPIX_;
+
+  wcsCRPIX_ = new Vector[MULTWCS];
+
+  Vector cc = center();
+  char crpix1[] = "CRPIX1 ";
+  char crpix2[] = "CRPIX2 ";
+  for (int ii=0; ii<MULTWCS; ii++) {
+    crpix1[6] = !ii ? ' ' : ii+'@';
+    crpix2[6] = !ii ? ' ' : ii+'@';
+    wcsCRPIX_[ii] = Vector(hd->getReal(crpix1,cc[0]),hd->getReal(crpix2,cc[1]));
+  }
+}
+
 void FitsImage::wcsSizeInit()
 {
   // init wcsSize_ array
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index c949167..7bd7b7d 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -115,6 +115,7 @@ class FitsImage {
   int* wcsCel_;
   int* wcsCelLon_;
   int* wcsCelLat_;
+  Vector* wcsCRPIX_;
   double* wcsSize_;
 
   WCSState* wcsState_;
@@ -150,6 +151,7 @@ class FitsImage {
 
   void wcsInit(int);
   void wcsCelInit(int);
+  void wcsCRPIXInit(FitsHead*);
   void wcsHPXInit();
   void wcsSizeInit();
   void wcsPhyInit();
@@ -405,6 +407,7 @@ class FitsImage {
   const char* getWCSAxisSymbol(Coord::CoordSystem, int);
 
   double getWCSSize(Coord::CoordSystem);
+  Vector getWCSCRPIX(Coord::CoordSystem);
   double calcWCSSize(Coord::CoordSystem);
 
   int hasWCS(Coord::CoordSystem);
diff --git a/tksao/frame/frblt.C b/tksao/frame/frblt.C
index 17361f4..6a6c2c8 100644
--- a/tksao/frame/frblt.C
+++ b/tksao/frame/frblt.C
@@ -284,7 +284,8 @@ int Base::markerAnalysisPlot3d(Marker* pp, double** x, double** y,
   SETSIGBUS
     for (int kk=0; kk<srcd; kk++) {
       double ss = kk+.5+.5+zparams->zmin;
-      Vector3d dd = Vector3d(ptr->center(),ss) * Translate3d(-.5,-.5,-.5);
+      Vector3d dd = Vector3d(ptr->getWCSCRPIX(sys),ss) *
+	Translate3d(-.5,-.5,-.5);
       Vector3d out = ptr->mapFromRef(dd,sys,sky);
       (*x)[kk] = out[2];