Close #1701 : SURFACE obeys multiplot position when T3D is not enabled. (#1798)

OK for me. At least for some basics tests. 
Thanks @GillesDuvert

Author: GillesDuvert

src/gdlgstream.cpp

@@ -1398,7 +1398,24 @@ void GDLGStream::getSubpageRegion(PLFLT *sxmin, PLFLT *symin, PLFLT *sxmax, PLFL
    *szmax=*szmin+profund;
  }
 }
-//get region (3BPP data)
+
+void GDLGStream::getCurrentSubpageInfo(PLFLT &xratio, PLFLT &yratio, PLFLT &zratio, PLFLT* displacement) {
+  int p=thePage.curPage-1;
+  int nx = thePage.nx;
+  int ny = thePage.ny;
+  xratio=1.0/thePage.nx;
+  yratio=1.0/thePage.ny;
+  zratio=1.0/thePage.nz;
+ int k= p / (thePage.nx*thePage.ny);
+ int l= p - k*(thePage.nx*thePage.ny);
+ int j= l /thePage.nx ;
+ int i= (l - j*thePage.nx);
+ displacement[0]=i*xratio;
+ displacement[1]=1-yratio-j*yratio;
+ displacement[2]=k*zratio;
+}
+
+  //get region (3BPP data)
 
 bool GDLGStream::GetRegion(DLong& xoff, DLong& yoff, DLong& nx, DLong& ny) {
   long nxOrig,nyOrig;

src/gdlgstream.hpp

@@ -640,6 +640,7 @@ class GDLGStream: public plstream
   void ssub( PLINT nx, PLINT ny, PLINT nz=1);
   void adv(PLINT page);
   void getSubpageRegion(PLFLT *sxmin, PLFLT *symin, PLFLT *sxmax, PLFLT *symax, PLFLT *zmin=NULL, PLFLT *zmax=NULL);
+  void getCurrentSubpageInfo(PLFLT &xratio, PLFLT &yratio, PLFLT &zratio, PLFLT* displacement);
   void SetPageDPMM(float setPsCharFudge=1.0, float setPsSymFudge=1.0);
   void syncPageInfo();
   void updateBoxDeviceCoords();

src/plotting.cpp

@@ -3976,24 +3976,55 @@ void SelfNormLonLat(DDoubleGDL *lonlat) {
     gdlAxis(e, a, YAXIS, yStart, yEnd, yLog);
   }
 
-  void gdlBox3(EnvT *e, GDLGStream *a, DDouble xStart, DDouble xEnd, bool xLog, DDouble yStart, DDouble yEnd, bool yLog, DDouble zStart, DDouble zEnd, bool zLog, DDouble zValue) {
+  void gdlBox3(EnvT *e, GDLGStream *a, DDouble xStart, DDouble xEnd, bool xLog, DDouble yStart, DDouble yEnd, bool yLog, DDouble zStart, DDouble zEnd, bool zLog, DDouble zValue, bool zAxis, DLong zaxis_value) {
     gdlAxis(e, a, XAXIS, xStart, xEnd, xLog, 1); //only Bottom
     gdlAxis(e, a, YAXIS, yStart, yEnd, yLog, 1); //only left
-
+	
     PLFLT xnormmin, xnormmax, ynormmin, ynormmax, znormmin = 0, znormmax = 1;
     a->getCurrentNormBox(xnormmin, xnormmax, ynormmin, ynormmax);
     gdlGetCurrentAxisWindow(ZAXIS, znormmin, znormmax);
 
     //almost cut and paste of plotting_axis section on Z axis:
-    PLFLT yPos = ynormmax;
-    PLFLT xPos = xnormmin;
-    PLFLT viewportXSize = xnormmax - xPos;
-    if (viewportXSize < 0.001) viewportXSize = 0.01;
-    a->vpor(xPos, xPos + viewportXSize, znormmin, znormmax);
-    //special transform to use 'y' axis code, but with 'z' values and yz exch.
-    gdlSetZto3DDriverTransform(a, yPos);
-    gdlExchange3DDriverTransform(a);
-    gdlAxis(e, a, ZAXIS, zStart, zEnd, zLog, 1);
+      if (zAxis) { //no use drawing something unseen, but IDL does, and wee need values if ztick_get is asked for!!!
+	      PLFLT yPos ;
+	      PLFLT xPos ;
+	      PLFLT viewportXSize ;
+        DLong modifierCode=0;
+        switch (zaxis_value) {
+        case 1:
+          modifierCode = 2;
+          yPos = ynormmin;
+          xPos = xnormmax;
+          viewportXSize=xPos-xnormmin; if (viewportXSize<0.001) viewportXSize=0.01;
+          a->vpor(xPos - viewportXSize, xPos, znormmin, znormmax);
+          break;
+        case 4:
+          modifierCode = 2;
+          yPos = ynormmax;
+          xPos = xnormmax;
+          viewportXSize=xPos-xnormmin; if (viewportXSize<0.001) viewportXSize=0.01;
+          a->vpor(xPos - viewportXSize, xPos, znormmin, znormmax);
+          break;
+        case 2:
+          modifierCode = 1;
+          yPos = ynormmin;
+          xPos = xnormmin;
+          viewportXSize=xnormmax-xPos; if (viewportXSize<0.001) viewportXSize=0.01;
+          a->vpor(xPos, xPos + viewportXSize, znormmin, znormmax);
+          break;
+        default: //->3 = 0
+          modifierCode = 1;
+          yPos =ynormmax;
+          xPos =xnormmin;
+          viewportXSize=xnormmax-xPos; if (viewportXSize<0.001) viewportXSize=0.01;
+          a->vpor(xPos, xPos + viewportXSize, znormmin, znormmax);
+          break;
+        }
+        //special transform to use 'y' axis code, but with 'z' values and yz exch.
+		gdlSetZto3DDriverTransform(a, yPos);
+        gdlExchange3DDriverTransform(a);
+        gdlAxis(e, a, ZAXIS, zStart, zEnd, zLog, modifierCode);
+      }
 
     //restore (normally not useful?)
     gdlExchange3DDriverTransform(a);

src/plotting.hpp

@@ -234,7 +234,9 @@ namespace lib {
   void SelfPDotTTransformXYZ(DDoubleGDL *xt, DDoubleGDL *yt, DDoubleGDL *zt);
   void PDotTTransformXYZval(PLFLT x, PLFLT y, PLFLT *xt, PLFLT *yt, PLPointer data);
   DDoubleGDL* gdlDefinePlplotRotationMatrix(DDouble az, DDouble alt, DDouble *scale, bool save);
-  bool gdlInterpretT3DMatrixAsPlplotRotationMatrix(DDouble &az, DDouble &alt, DDouble &ay, DDouble *scale, T3DEXCHANGECODE &axisExchangeCode, bool &below);
+  void gdlMakeSubpageRotationMatrix3d(DDoubleGDL* me, PLFLT xratio, PLFLT yratio, PLFLT zratio, PLFLT* trans);
+  void gdlMakeSubpageRotationMatrix2d(DDoubleGDL* me, PLFLT xratio, PLFLT yratio, PLFLT zratio, PLFLT* trans);
+  bool gdlInterpretT3DMatrixAsPlplotRotationMatrix(DDouble &az, DDouble &alt, DDouble &ay, DDouble *scale, /* DDouble *trans,*/  T3DEXCHANGECODE &axisExchangeCode, bool &below);
   DDoubleGDL* gdlGetT3DMatrix();
   void gdlStartT3DMatrixDriverTransform(GDLGStream *a, DDouble zValue);
   void gdlStartSpecial3DDriverTransform( GDLGStream *a, GDL_3DTRANSFORMDEVICE &PlotDevice3D);
@@ -436,7 +438,7 @@ namespace lib {
 
   void gdlBox(EnvT *e, GDLGStream *a, DDouble xStart, DDouble xEnd, bool xLog, DDouble yStart, DDouble yEnd, bool yLog);
 
-  void gdlBox3(EnvT *e, GDLGStream *a, DDouble xStart, DDouble xEnd, bool xLog, DDouble yStart, DDouble yEnd, bool yLog,  DDouble zStart, DDouble zEnd, bool zLog, DDouble zValue);
+  void gdlBox3(EnvT *e, GDLGStream *a, DDouble xStart, DDouble xEnd, bool xLog, DDouble yStart, DDouble yEnd, bool yLog,  DDouble zStart, DDouble zEnd, bool zLog, DDouble zValue, bool drawZ=false, DLong zcode=0);
 
   //restore current clipbox, make another or remove it at all.
   bool gdlSwitchToClippedNormalizedCoordinates(EnvT *e, GDLGStream *actStream, bool invertedClipMeaning=false, bool commandHasCoordSys=true );

src/plotting_convert_coord.cpp

@@ -646,11 +646,13 @@ namespace lib {
   // create all 8 homogenous coordinates of the cubes points, 
   // compute their projected coord,
   // find min and max on both axes.
-    static DDoubleGDL* cube_coord = (new DDoubleGDL(dimension(8,4), BaseGDL::NOZERO));
-    static DDouble vals[32]={0,1,0,1,0,1,0,1,0,0,1,1,0,0,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1};
-    for (auto i=0; i< 32; ++i) (*cube_coord)[i]=vals[i];
+	static DDouble vals[32] = {0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+//    static DDoubleGDL cube_coord(dimension(8,4), BaseGDL::NOZERO);
+    static DDoubleGDL cube_coord(vals,32);
+	cube_coord.SetDim(dimension(8,4)); //needed as previous initializaton is 1D only.
+//    for (auto i=0; i< 32; ++i) (*(&cube_coord))[i]=vals[i];
     //    w=v#!p.t
-    DDoubleGDL* tmp = plplot3d->MatrixOp(cube_coord, false, true);
+    DDoubleGDL* tmp = plplot3d->MatrixOp(&cube_coord, false, true);
     DDoubleGDL* res = static_cast<DDoubleGDL*>(tmp->Transpose(NULL));
     GDLDelete(tmp);
 //    SelfPrint3d(res);
@@ -718,6 +720,22 @@ DDoubleGDL* gdlDoAsScale3(DDouble az, DDouble alt, DDouble *scalex, DDouble *sca
     return t3dMatrix; 
   }
 
+   //this displaces 3d elements made in NoSubpage mode to a particular subpage. Used for axis for example.
+  void gdlMakeSubpageRotationMatrix3d(DDoubleGDL* me, PLFLT xratio, PLFLT yratio, PLFLT zratio, PLFLT* trans) {
+	DDouble newscale[3]={ xratio, yratio, zratio};
+    SelfTranspose3d(me);
+	SelfScale3d(me,newscale);
+	SelfTranslate3d(me,trans);
+    SelfTranspose3d(me);
+  }
+  //this displaces 2D elements (images, already reprojected 3d) made in NoSubpage mode to a particluar subpage
+  void gdlMakeSubpageRotationMatrix2d(DDoubleGDL* me, PLFLT xratio, PLFLT yratio, PLFLT zratio, PLFLT* trans) {
+	DDouble newscale[3]={ xratio, yratio, zratio};
+	SelfReset3d(me);
+	SelfScale3d(me,newscale);
+	SelfTranslate3d(me,trans);
+    SelfTranspose3d(me);
+  } 
  //converts 3D values according to COORDSYS towards NORMAL coordinates and , logically, unset xLog,yLo,zLog and define code as NORMAL.
   void SelfConvertToNormXYZ(DDoubleGDL *x, bool &xLog, DDoubleGDL *y, bool &yLog, DDoubleGDL *z, bool &zLog, COORDSYS &code) {
 //    std::cerr<<"SelfConvertToNormXYZ(DDoubleGDL)\n";
@@ -981,7 +999,7 @@ bool isAxonometricRotation(DDoubleGDL* Matrix, DDouble &alt, DDouble &az, DDoubl
   // Note that if drawing of axes and hidden surfaces (by plplot) are not needed (e.g., PLOTS), any matrix can be used.
   // Returns a 'plplot-compatible' matrix that will be used in calls to plplot.
   // Retunrs NULL if conversion is impossible.
-  bool gdlInterpretT3DMatrixAsPlplotRotationMatrix(DDouble &az, DDouble &alt, DDouble &ay, DDouble *scale, T3DEXCHANGECODE &axisExchangeCode, bool &below) {
+  bool gdlInterpretT3DMatrixAsPlplotRotationMatrix(DDouble &az, DDouble &alt, DDouble &ay, DDouble *scale, /*DDouble *trans,*/ T3DEXCHANGECODE &axisExchangeCode, bool &below) {
 //    std::cerr<<"gdlInterpretT3DMatrixAsPlplotRotationMatrix(()\n";
     //returns NULL if error!
     DDoubleGDL* t3dMatrix = new DDoubleGDL(dimension(4, 4),BaseGDL::NOZERO);
@@ -991,6 +1009,12 @@ bool isAxonometricRotation(DDoubleGDL* Matrix, DDouble &alt, DDouble &az, DDoubl
     static unsigned tTag = pStruct->Desc()->TagIndex("T");
     for (int i = 0; i < t3dMatrix->N_Elements(); ++i)(*t3dMatrix)[i] = (*static_cast<DDoubleGDL*> (pStruct->GetTag(tTag, 0)))[i];
     SelfTranspose3d(t3dMatrix);
+//	//compute translation in vport space using the displacement of [0,0,0]:
+//	static DDouble center[4]={0.5,0.5,0.5,1.};
+//	static DDoubleGDL zeropos(center,4);
+//	DDoubleGDL* intermediary = t3dMatrix->MatrixOp(&zeropos, false, false);
+//	for (int i=0; i< 3; ++i) trans[i]= (*intermediary)[i]-0.5;
+//	GDLDelete(intermediary);	
     //check if valid, get rotation etc.
     if (!isAxonometricRotation(t3dMatrix, alt, az, ay, scale, axisExchangeCode, below)) return false;
     return true;

src/plotting_shade_surf.cpp

@@ -47,7 +47,8 @@ namespace lib
     bool below=false;
     DDouble* Current3DMatrix;
     DDoubleGDL* gdlBox3d;
-    
+	PLFLT xratio, yratio, zratio, trans[3];
+	PLINT nx, ny;    
  private:
     bool handle_args (EnvT* e)
     {
@@ -231,6 +232,14 @@ namespace lib
     {
       static int savet3dIx = e->KeywordIx("SAVE");
       bool saveT3d = e->KeywordSet(savet3dIx);
+	  bool zAxis=true;
+	  static int zaxisIx = e->KeywordIx("ZAXIS");
+	  DLong zaxis_value=0;
+      if (e->GetKW(zaxisIx) != NULL) {
+        e->AssureLongScalarKWIfPresent(zaxisIx, zaxis_value);
+		if (zaxis_value > 4) zaxis_value=0;
+	  }
+	  if (zaxis_value < 0) zAxis = false;
 
       // background BEFORE next plot since it is the only place plplot may redraw the background...
       gdlSetGraphicsBackgroundColorFromKw ( e, actStream );
@@ -266,9 +275,14 @@ namespace lib
           below=true;
           alt=-(360.-alt);
       }
-        //Compute special transformation matrix for the BOX and give it to the driver
+		//Compute special transformation matrix for the BOX and give it to the driver. Subpage info is important
         gdlBox3d=gdlDefinePlplotRotationMatrix( az, alt, scale, saveT3d);
-        Guard<BaseGDL> g(gdlBox3d);
+		actStream->getCurrentSubpageInfo(xratio, yratio, zratio, trans);
+		gdlMakeSubpageRotationMatrix3d(gdlBox3d, xratio, yratio, zratio,trans);
+		//now that 3D matrix is OK, we pass in 'No Sub'
+		actStream->NoSub();
+		zValue=gdlSetViewPortAndWorldCoordinates(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue);
+
         GDL_3DTRANSFORMDEVICE T3DForAXes;
         for (int i = 0; i < 16; ++i)T3DForAXes.T[i] =(*gdlBox3d)[i];
         T3DForAXes.zValue = (std::isfinite(zValue))?zValue:0;
@@ -290,7 +304,7 @@ namespace lib
       //Should draw 3d mesh before axes
       gdlSetGraphicsForegroundColorFromKw ( e, actStream ); //COLOR
       //write OUR box using our 3D PLESC tricks:
-      gdlBox3(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue);
+      gdlBox3(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue, zAxis, zaxis_value);
        //reset driver to 2D plotting routines, further 3D is just plplot drawing a mesh.
       gdlStop3DDriverTransform(actStream); 
 
@@ -320,7 +334,8 @@ namespace lib
 //This is the good version for shade_surf, shifting the shaded surface by som amount in the 3DDriverTransform of the driver, not the plplot library.     
       actStream->vpor(0,1,0,1);
       actStream->wind(-0.5/scale[0],0.5/scale[0],-0.5/scale[1],0.5/scale[1]); //mandatory: to center in (0,0,0) for 3D Matrix rotation.
-      if (below) { actStream->w3d(1,1,1,0,1,0,1,0,1, -alt, az);
+      if (below) {
+		actStream->w3d(1,1,1,0,1,0,1,0,1, -alt, az);
         DDouble xp = 0;
         DDouble yp1 = 0;
         DDouble yp2 = 0;
@@ -340,6 +355,14 @@ namespace lib
         Matrix3DTransformXYZval(0, 0, 0.5, &xp, &yp2,Current3DMatrix);
         gdlShiftYaxisUsing3DDriverTransform(actStream, yp1 - yp2, false);
       }
+		if (!doT3d) { //use a special matrix to shift and scale into current subpage
+		  gdlMakeSubpageRotationMatrix2d(gdlBox3d, xratio, yratio, zratio,trans);
+		  GDL_3DTRANSFORMDEVICE T3DForAXes;
+		  for (int i = 0; i < 16; ++i)T3DForAXes.T[i] = (*gdlBox3d)[i];
+		  T3DForAXes.zValue = (std::isfinite(zValue)) ? zValue : 0;
+		  gdlStartSpecial3DDriverTransform(actStream, T3DForAXes);
+          Guard<BaseGDL> g(gdlBox3d);
+		}
 
       return false;
       }

src/plotting_surface.cpp

@@ -49,7 +49,8 @@ namespace lib
     bool below=false;
     DDouble* Current3DMatrix;
     DDoubleGDL* gdlBox3d;
-    
+	PLFLT xratio, yratio, zratio, trans[3];
+	PLINT nx, ny;    
   private:
     bool handle_args (EnvT* e)
 {
@@ -233,7 +234,15 @@ namespace lib
     {
       static int savet3dIx = e->KeywordIx("SAVE");
       bool saveT3d = e->KeywordSet(savet3dIx);
-
+	  bool zAxis=true;
+	  static int zaxisIx = e->KeywordIx("ZAXIS");
+	  DLong zaxis_value=0;
+      if (e->GetKW(zaxisIx) != NULL) {
+        e->AssureLongScalarKWIfPresent(zaxisIx, zaxis_value);
+		if (zaxis_value > 4) zaxis_value=0;
+	  }
+	  if (zaxis_value < 0) zAxis = false;
+ 
       // background BEFORE next plot since it is the only place plplot may redraw the background...
       gdlSetGraphicsBackgroundColorFromKw ( e, actStream );
       //start a plot
@@ -267,10 +276,15 @@ namespace lib
         } else if (alt > 270) {
           below=true;
           alt=-(360.-alt);
-        }
-        //Compute special transformation matrix for the BOX and give it to the driver
+		}
+		//Compute special transformation matrix for the BOX and give it to the driver. Subpage info is important
         gdlBox3d=gdlDefinePlplotRotationMatrix( az, alt, scale, saveT3d);
-        Guard<BaseGDL> g(gdlBox3d);
+		actStream->getCurrentSubpageInfo(xratio, yratio, zratio, trans);
+		gdlMakeSubpageRotationMatrix3d(gdlBox3d, xratio, yratio, zratio,trans);
+		//now that 3D matrix is OK, we pass in 'No Sub'
+		actStream->NoSub();
+		zValue=gdlSetViewPortAndWorldCoordinates(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue);
+
         GDL_3DTRANSFORMDEVICE T3DForAXes;
         for (int i = 0; i < 16; ++i)T3DForAXes.T[i] =(*gdlBox3d)[i];
         T3DForAXes.zValue = (std::isfinite(zValue))?zValue:0;
@@ -291,7 +305,7 @@ namespace lib
       //Draw axes with normal color!
       gdlSetGraphicsForegroundColorFromKw ( e, actStream ); //COLOR
       //write OUR box using our 3D PLESC tricks:
-      gdlBox3(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue);
+      gdlBox3(e, actStream, xStart, xEnd, xLog, yStart, yEnd, yLog, zStart, zEnd, zLog, zValue, zAxis, zaxis_value);
        //reset driver to 2D plotting routines, further 3D is just plplot drawing a mesh.
       gdlStop3DDriverTransform(actStream); 
 
@@ -324,13 +338,21 @@ namespace lib
 // shade //      if (e->KeywordPresent(shadesIx)) doShade = true;
 // shade //      
 // shade //      if (!doShade) {
-        //This is the good version for surface without the shade argument.     
+        //This is the good version for surface without the shade argument.
         actStream->vpor(0, 1, 0, 1);
-        actStream->wind(-0.5 / scale[0], 0.5 / scale[0], -0.5 / scale[1], 0.5 / scale[1]);
-        if (below) {
-          actStream->w3d(1, 1, 1, 0, 1, 0, 1, 0.5, 1.5, -alt, az);
-          gdlFlipYPlotDirection(actStream); //special trick, not possible with plplot
-        } else actStream->w3d(1, 1, 1, 0, 1, 0, 1, 0.5, 1.5, alt, az);
+        actStream->wind(-0.5/scale[0],0.5/scale[0],-0.5/scale[1],0.5/scale[1]); //mandatory: to center in (0,0,0) for 3D Matrix rotation.
+		if (!doT3d) { //use a special matrix to shift and scale into current subpage
+		  gdlMakeSubpageRotationMatrix2d(gdlBox3d, xratio, yratio, zratio,trans);
+		  GDL_3DTRANSFORMDEVICE T3DForAXes;
+		  for (int i = 0; i < 16; ++i)T3DForAXes.T[i] = (*gdlBox3d)[i];
+		  T3DForAXes.zValue = (std::isfinite(zValue)) ? zValue : 0;
+		  gdlStartSpecial3DDriverTransform(actStream, T3DForAXes);
+          Guard<BaseGDL> g(gdlBox3d);
+		}
+		if (below) {
+		  actStream->w3d(1, 1, 1, 0, 1, 0, 1, 0.5, 1.5, -alt, az);
+		  gdlFlipYPlotDirection(actStream); //special trick, not possible with plplot
+		} else actStream->w3d(1, 1, 1, 0, 1, 0, 1, 0.5, 1.5, alt, az);
 // shade //      } else {
 // shade //        //This is the good version for shade_surf and surface with shade option
 // shade //        // (needs shifting the plplot plot by some amount in the 3DDriverTransform of the driver.)