C++ Matrix::Clo方法代码示例

void MpDensityPlot::Draw (Scene& scene, const Matrix& z,
                          const ColorMap &cmap, double cmin, double cmax,
                          int mode)
{
  if (! scene.IsOpen()) 
    Matpack.Error("MpDensityPlot: scene is not open");

  int outlined = (mode & MpDensityPlot::Outlined) ? (Fill | Outline) : Fill;

  ColorB BlendFrom(0,0,0), BlendTo(255,255,255); // black to white
        
  // scale factor for color mapping
  int ncmap;
  double cscale = 0, 
         dz = cmax - cmin;
  if (cmap) {
    ncmap = cmap.Size()-1;
    cscale = (dz == 0.0) ? 0.0 : ncmap/dz;
  } else 
    ncmap = 0;
    
  // plot all density cells 
  for (int x = z.Clo(); x <= z.Chi(); ++x) {

    double dx1 = (x != z.Clo()) ? 0.5 : 0,
           dx2 = (x != z.Chi()) ? 0.5 : 0;
        
    for (int y = z.Rlo(); y <= z.Rhi(); ++y) {
                
      // calculate color
      long color;
      if (ncmap) {
        int c = int((z[y][x]-cmin)*cscale);
        MpForceRange(c,0,ncmap);
        color = scene.NewColor(cmap[c]);
              
      } else {    
        // linear map to color range black - white
        double intens = (dz == 0.0) ? 0.0 : (z[y][x] - cmin)/dz;
        MpForceRange(intens,0.0,1.0);
        color = scene.NewColor(((1.0-intens)*BlendFrom.red+intens*BlendTo.red)/255.0,
                               ((1.0-intens)*BlendFrom.green+intens*BlendTo.green)/255.0,
                               ((1.0-intens)*BlendFrom.blue+intens*BlendTo.blue)/255.0, RGB);
      }
            
      // draw squares
      double dy1 = (y != z.Rlo()) ? 0.5 : 0,
             dy2 = (y != z.Rhi()) ? 0.5 : 0;
      double px[4] = {y-dy1,y+dy2,y+dy2,y-dy1},
             py[4] = {x-dx1,x-dx1,x+dx2,x+dx2};
      scene.Polygon(4,px,py,color,outlined);
    }
  }
}

static void MakeWater (Matrix& M, float waterlevel = 0.0,
		float scale_x = 80, float scale_y = 120, float scale_z = 0.1)
{
  for (int i = M.Rlo(); i <= M.Rhi(); i++)
    for (int j = M.Clo(); j <= M.Chi(); j++)
      if (M(i,j) < waterlevel) {		   // apply noise to values below
	float x = scale_x * (i-M.Rlo())/M.Rows(),  // the water level to simulate
	      y = scale_y * (j-M.Clo())/M.Cols();	      // waves on the sea
	M(i,j) =  scale_z * Noise(x,y,0);
     }
}

int main (void)
{
  enum { N = 20, L = 10, resolution = 600};

  Matrix M = TestMatrix(N,N);

  Vector level(1,L);
  double mi = Min(M), ma = Max(M);
  for (int i = 1; i <= L; ++i) 
    level[i] = mi + (i-1)*(ma-mi)/(L-1);
  level[4]  = -3; 
  level[5]  = -2.5;
  cerr << mi << " " << ma << endl <<  level << endl;

  ColorMap cmap("spectral-256");
  MpImage image(resolution,resolution);
  Scene scene(image);

  MpContourSurface(scene,M, M.Rlo(),M.Rhi(), M.Clo(),M.Chi(), 
		   level,cmap,Min(M),Max(M));

  scene.LookFrom(Vector3D(100,50,-30));    // define the camera position
  scene.SetColoring(Scene::PerVertex);     // coloring (GLobal/PerVertex/PerFacet)
  scene.SetGlobalShading(Facet::Flat);     // shading (Constant/Flat/Gouraud/Phong)
  scene.SetGlobalEdgeLines(true);	   // edge lines (true/false)

  scene.Show();				   // now start rendering of the surface

  MpFrame display("Contour",resolution+4,resolution+4);
  new MpScrollImageWindow(display,image,resolution+4,resolution+4);
  new MpQuitButton(display,"Quit",80,25,0,0);
  display.EventLoop();
}

void MpViewGraphData::DrawPolesPlot (Scene& scene, const Matrix& Z, 
				     const Trafo& P)
{
  int identical_colors = (Three_Dim_Poles_Color == Three_Dim_Poles_Head_Color);
  int r = Three_Dim_Poles_Head_Radius;

  // set color only once
  if (identical_colors) scene.SetColor(Three_Dim_Poles_Color);

  // set lines style for poles
  scene.LineStyle(Solid,Three_Dim_Poles_Width);

  for (int x = Z.Rlo(); x <= Z.Rhi(); x++)
    for (int y = Z.Clo(); y <= Z.Chi(); y++) {

      // draw line of pole
      if ( ! identical_colors ) scene.SetColor(Three_Dim_Poles_Color);
      if (Three_Dim_Poles_Width)
	scene.Line( P * Vector3D(y,Three_Dim_Poles_Base,x), 
		    P * Vector3D(y,Z[x][y],x) );

      // draw head of pole
      if (Three_Dim_Poles_Head_Radius > 0) {
	if ( ! identical_colors ) scene.SetColor(Three_Dim_Poles_Head_Color);
	scene.Arc( scene.Map(P * Vector3D(y,Z[x][y],x)), 
		   r,r, 0, 360*64, Fill);
      }
    }
}

void MpSurfacePrism (Scene &S, const Matrix &Z, int nsegments,
                     double wx, double wy, 
                     const ColorMap &cmap,double cmin, double cmax)
{
  // nothing that can be drawn
  if ( Z.Empty() ) return;

  // 4-segments faster special case
  if (nsegments == 4) {
    MpSurfacePrism4(S,Z,wx,wy,cmap,cmin,cmax);
    return;
  }

  // find the extrema
  double zmin = Min(Z);

  // check for valid nsegments values
  MpForceRange(nsegments, 3, MaxVertices);

  // check if a colormap should be used
  double dz,cscale = 0;
  bool usecolormap;
  int ncmap = cmap.Size()-1;
  if (cmap) { 
    usecolormap = true;
    dz = cmax - cmin; 
    cscale = (dz == 0.0) ? 0 : ncmap/dz;
  } else
    usecolormap = false;

  // add prism
  for (int x = Z.Rlo(); x <= Z.Rhi(); x++)
    for (int y = Z.Clo(); y <= Z.Chi(); y++) {

      double z = Z(x,y), 
             height = 0.5 * (z - zmin);
      Trafo T = trans(y,height+zmin,x)
                * scale(wy,height,wx)
                * rot(Trafo::XAxis,M_PI/2);

      // coloring
      const ColorB *color;
      if (usecolormap) { 
        int c = int((z-cmin)*cscale);
        if (c > ncmap) 
          c = ncmap;
        else if (c < 0) 
          c = 0;
        color = &cmap[c];
      } else
        color = (ColorB*)NULL;

      cylinder(S,nsegments,T,color);
  }
}

void FractalSurface (Matrix& X, int b, double H, long seed)
{
    // The dimension of the matrix is 2^b + 1. Resize if neccessary 
    // (don't do anything if size is ok), then reset to zero.
    int n = 1 << b;

    if (X.Rlo() != 0 || X.Clo() != 0 || X.Rhi() != n || X.Chi() != n) {
      X.Remove();
      X = Matrix(0,n,0,n);
    }

    // clear
    X = 0.0;
    
    // Create normal distributed deviates with zero mean and unit variance
    double mean = 0, sdev = 1;
    Ran002 rnd(seed);
    NormalDistribution N(mean, sdev, &rnd);
         
    // Voss random addition algorithm
    int stride = n/2, off = stride/2;

    // start
    for (int i = 0; i <= n; i += stride)
      for (int j = 0; j <= n; j += stride) X(i,j) = N();

    // subdivide
    for (int g = 1; g < b; g++) {
	int i,j;

	// interpolate
	for (i = off; i <= n-off; i += stride)
	  for (j = off; j <= n-off; j += stride)
	    X(i,j) = 0.25 * ( X(i-off,j-off) + X(i-off,j+off) 
			    + X(i+off,j+off) + X(i+off,j-off));	
	for (i = 0; i <= n; i += stride)
	  for (j = off; j <= n-off; j += stride)
	    X(i,j) = 0.5 * (X(i,j-off) + X(i,j+off));
	for (i = off; i <= n-off; i += stride)
	  for (j = 0; j <= n; j += stride)
	    X(i,j) = 0.5 * (X(i-off,j) + X(i+off,j));

	// add noise with reduced variance
	sdev *= pow(0.5,H);
	N = NormalDistribution(mean, sdev, &rnd);
	for (i = 0; i <= n; i += off) 
	  for (j = 0; j <= n; j += off) X(i,j) += N(); 

	// next subdivision
	stride = off; 
	off /= 2;
    } 
}

static void MpSurfacePrism4 (Scene &S, const Matrix &Z, 
                             double wx, double wy, 
                             const ColorMap &cmap, double cmin, double cmax)
{
  // nothing that can be drawn
  if ( Z.Empty() ) return;
  
  // find the minimum
  double zmin = Min(Z);
  
  // check if a colormap should be used
  double dz,cscale = 0;
  bool usecolormap;
  int ncmap = cmap.Size()-1;
  if (cmap) {
    usecolormap = true;
    dz = cmax - cmin; 
    cscale = (dz == 0.0) ? 0.0 : ncmap/dz;
  } else
    usecolormap = false;
  
  // add square prisms
  for (int x = Z.Rlo(); x <= Z.Rhi(); x++)
    for (int y = Z.Clo(); y <= Z.Chi(); y++) {

      double z = Z(x,y), 
             height = 0.5 * (z - zmin);
      Trafo T = trans(y,height+zmin,x) 
                * scale(wy,height,wx);

      // coloring
      const ColorB *color;
      if (usecolormap) { 
        int c = int((z-cmin)*cscale);
        if (c > ncmap) 
          c = ncmap;
        else if (c < 0) 
          c = 0;
        color = &cmap[c];
      } else
        color = (ColorB*)NULL;
      
      MpCube(S,T,color);  
    }
}

int main (void)
{
  MpImage image(size,size);		     // allocate size x size raster image

  MpImage backimg;			                 // read background image
  if (back_image) {
    if ( ! backimg.ReadAnyFile(back_image) )
      Matpack.Error("Can't read \"%s\"",back_image);
  }

  // read potential matrix
  ifstream pot_stream("pot");
  if ( ! pot_stream) Matpack.Error("Can't open pot data file \"%s\"","pot");
  pot_stream >> pot;
  pot_stream.close();

  for (int i = 0; i <= 175; i += 5) {  			   // loop through frames

    fprintf(stdout,"\rframe %d...",i); fflush(stdout);         // what's going on

    // read complex wave function
    char file_name[200], pipe_name[200], output_name[200];
    sprintf(file_name,"psi%d.gz",i);		           // generate input data file name
    sprintf(pipe_name,"gunzip -c %s > psi.dat",file_name); // decompress data file 
    sprintf(output_name,"psi%03d.jpg",i);                  // generate output file name

    system(pipe_name);
    ifstream data_stream("psi.dat");
    if ( ! data_stream) Matpack.Error("Can't open data file psi.dat");
    data_stream >> psi;
    data_stream.close();
    unlink("psi.dat");

    // consistency checks
    if ( pot.Rlo() != psi.Rlo() || pot.Rhi() != psi.Rhi() || 
	 pot.Clo() != psi.Clo() || pot.Chi() != psi.Chi() )
      Matpack.Error("non-conformant index range of potential and wave function");

    Scene scene;				   // define scene for 3d drawing

    // create surface
    double u0 = psi.Clo(), u1 = psi.Chi(), 
           v0 = psi.Rlo(), v1 = psi.Rhi();
    int    nu = psi.Chi()-psi.Clo(), nv = psi.Rhi()-psi.Rlo(),  
           su = 0, sv = 0,  periodic = 0;
    ParametricSurface(scene, height_fcn, u0,u1,nu,su, v0,v1,nv,sv, 
		      periodic, Identity, color_fcn);

    scene.BoxRatios(1,z_aspect,1);               // scale scene to box with x:y:z

    float dist = 1.6;			                    // distance parameter
    scene.Look(Vector3D( 1.1*dist, 1.5*dist, 1.4*dist), // camera position vector
	       Vector3D(-1.1*dist,-1.65*dist,-1.4*dist), // look direction vector
	       FieldOfView(45), 		 // field of view angle in degree
	       0);		             // "twist your head" angle in degree

    scene.SetGlobalShading(Facet::Gouraud);   
                                       // shading (None, Flat, Gouraud, Phong...)
    scene.SetHiddenSurface(preview ? Scene::DepthSort : Scene::TopoSort);

    // edgeline drawing option (None,...)
    scene.SetGlobalEdgeLines(show_grid ? Facet::Individual : None);

    // Setting per-vertex coloring is the best choice to get smoothly changing
    // colors on the surface. This requires Gouraud or Phong shading to be
    // switched on. If you use per-facet coloring then the facet edges are still
    // slightly visible ("Scene::PerFacet")
    scene.SetColoring(Scene::PerVertex);

    // Define ambient light, RGB color, and specular light and exponent
    // If per-vertex coloring or per-facet coloring is set, then only the ambient
    // light factor and the specular light factor and exponent are in effect, and
    // the diffuse color is ignored (the vertex or facet color is used instead)
    Material material(0.6, ColorF(1.,0.8,0.), 0.3,1.8); 
    scene.SetFacetMaterial(material);		    // set all facets to material

    scene.Open(image);			        // direct drawing to raster image
    if (back_image)
      image.InsertTiled(backimg);
    else
      scene.SetBackground(back_color); 	        // draw background with RGB color
    scene.SetColor(edge_color);		           // define color for edge lines
    scene.Show();			  		    // render the surface
    scene.Close();				       // call always after close

    image.WriteJpegFile(output_name,jpeg_quality,jpeg_smooth);// write JPEG image

    // write GIF image, if you prefer that
    //image.WriteGifFile(output_name);     // write image as GIF

  }

  cout << "\nfinished" << endl;
}

void ContourLine (Scene& scene, const Matrix& z, 
                  double h, int mode, double hxy, const Trafo& trafo, 
                  char* label, double label_distance)
//
//  Plots a single contour line at a given height 'h' in the
//  data set z. Depending on the mode flag
//  it uses a linear interpolation or Bezier polynomials to draw
//  smooth contour lines.
//
//  Possible values for the mode flag are:
//
//    ContourPlotClass::Lines
//                        to draw a contour line with straight line segments
//
//    ContourPlotClass::BezierSpline
//                        to draw a smooth contour line using cubic 
//                        Bezier polynomials (this is the default)
//
//    ContourPlotClass::CardinalSpline
//                        draw smooth lines using cubic Cardinal Splines
//
//    Project3d           to project the contour line in 3d space onto the 
//                        xy-plane at height hxy. (default is 2d).
//
//  The mode flags can be combined like 'Linear|Project3D'
//
{
    int x,y;

    if (! scene.IsOpen()) 
      Matpack.Error(Mat::UnspecifiedError,"ContourLine: scene is not open");

    // make some parameters global 
    S = &scene;
    P = trafo;
    height = h;
    xyheight = hxy;
    xdim = z.Chi()-z.Clo()+2;
    ydim = z.Rhi()-z.Rlo()+2;
    xoff = z.Clo()-1;
    yoff = z.Rlo()-1;
    lab = label;
    lab_dist = label_distance;

    // straight lines or Bezier polynomials are drawn 
    drawlinear = mode & ContourPlotClass::Lines; 
    drawbspline = mode & ContourPlotClass::CardinalSpline;

    // use 3D or 2D projection
    project3d  = (mode & Project3D);
    
    // set the resolution global variable - must at least be = 3 !!
    resol = scene.curve.resolution + 2; 

    // allocate a cut-array 
    cut = bytematrix(0,xdim,0,ydim);

    // clear cut-array boundaries 
    for (x = 0; x <= xdim; x++)
      cut[x][0] = cut[x][ydim] = lower;
    for (y = 0; y <= ydim; y++)
      cut[0][y] = cut[xdim][y] = lower;

    // determine cut-array 
    for (x = 1; x < xdim; x++)
      for (y = 1; y < ydim; y++)
        cut[x][y] = (z[y+yoff][x+xoff] >= height) ? higher : lower;

#ifdef DEBUG
    printf("\nBEFORE\n");
    for (x = 0; x <= xdim; x++) {
        for (y = 0; y <= ydim; y++) 
          printf("%c",cut[x][y]==0?'.':cut[x][y]==1?'+':'2');
        printf("\n");
    }
    printf("\n");
#endif    
    
    // find contour lines  
    // ------------------

    for (x = 1; x < xdim; x++)
      for (y = 1; y < ydim; y++)
        if (cut[x][y] == higher) 
          if (cut[x][y-1] == lower)      // initial search direction = 0
            search_contour(z,x,y-1,0); 
          else if (cut[x][y+1] == lower) // initial search direction = 2
            search_contour(z,x,y+1,2); 

#ifdef DEBUG
    printf("\nAFTER\n");
    for (x = 0; x <= xdim; x++) {
        for (y = 0; y <= ydim; y++) 
          printf("%c",cut[x][y]==0?'.':cut[x][y]==1?'+':'2');
        printf("\n");
    }
    printf("\n");
#endif    

    // flush graphics buffer 
//.........这里部分代码省略.........

int MpWarpMatrixToSphere (Scene& S,
                          Matrix& H,
                          Vector& Pol, Vector& Azm,
                          double r,
                          const Trafo& T,
                          int facet_culling,
                          bool periodic,
                          bool south_pole_cap,
                          double h_south,
                          bool north_pole_cap,
                          double h_north,
                          const ColorMap& cmap,
                          double cmin, double cmax,
                          int PolStride, int AzmStride)
{
  int a,p,ncmap = 0;
  double x,y,z,pol,hr,sp,cp;
  double dz,scale = 0.0;
  bool usecolormap;

  int p1 = H.Rlo(), p2 = H.Rhi(), 
      a1 = H.Clo(), a2 = H.Chi();

  // check matrix and vector dimensions
  if (p1 != Pol.Lo() || p2 != Pol.Hi() ||
      a1 != Azm.Lo() || a2 != Azm.Hi())
    return 1;

  // check edge line stride values
  if (PolStride <= 0 || AzmStride <= 0) return 4;

  // evaluate facet culling flag
  int addFront = (facet_culling == Scene::Front) || (facet_culling == false),
      addBack  = (facet_culling == Scene::Back)  || (facet_culling == false);

  // number of facets in polar and azimuthal direction
  int np = p2-p1+1,
      na = a2-a1+1;

  // nothing that can be drawn
  if (np < 2 || na < 2) return 3;
    
  // pointer to vertices in two consecutive slices
  Vertex **top,**bot=0,**act;

  // allocate auxilliary storage
  if ( !(top = new Vertex*[na]) || !(bot = new Vertex*[na]) )
    Matpack.Error("WarpMatrixToSphere: vector allocation failure");
 
  // check if a colormap should be used
  if (cmap) { 
    usecolormap = true;
    dz = cmax-cmin; 
    ncmap = cmap.Size()-1;
    scale = (dz == 0.0) ? 0 : ncmap/dz;
  } else
    usecolormap = false;
        
  // Precompute sines and cosines of the azimuthal angles for the
  // sake of high performance
  Vector s,c;
  s = Sin(Azm);
  c = Cos(Azm);

  // loop through latitudes (polar angle)
  for ( p = p1, act = top; p <= p2; p++ ) {

    pol = Pol(p);
    sp  = sin(pol);
    cp  = cos(pol);

    // loop through longitudes (azimuthal angle)
    for ( a = a1; a <= a2; a++ ) {

      // warp transformation to sphere
      hr  = H(p,a) + r;
      x =  hr * cp * s(a);
      y = -hr * cp * c(a);
      z =  hr * sp;

      // add vertices to scene storing a pointer to it
      // (and apply the global transformation)
      act[a-a1] = &S.AddVertex( T * Vertex((float)x,(float)z,(float)y) );  

      // set the vertex color
      if (usecolormap) {
        int cen = int((H(p,a)-cmin)*scale);
	MpForceRange(cen,0,ncmap);
        act[a-a1]->SetColor(cmap[cen]);
      }
    }

    // first slice or last slice - add triangles to create pole caps
    if ( (south_pole_cap && p == p1) || (north_pole_cap && p == p2) ) { 

      int addF = false,
          addB = false;
      double h = 0.0;
      if (p == p1 && south_pole_cap) {
        addF = addBack;
//.........这里部分代码省略.........

Matrix operator * (const Matrix& A, const Matrix& B)
{
    if (A.Clo() != B.Rlo() || A.Chi() != B.Rhi()) 
      Matpack.Error("Matrix operator * (const Matrix&, const Matrix&): "
                    "non conformant arguments\n");

    // allocate return matrix
    Matrix C(A.Rlo(),A.Rhi(),B.Clo(),B.Chi());
    
    //------------------------------------------------------------------------//
    // the BLAS version
    //------------------------------------------------------------------------//

#if defined ( _MATPACK_USE_BLAS_ )

    if ( LT(B) ) {                   // full matrix * lower triangle
#ifdef DEBUG
        cout << "GM*LT\n";
#endif
        checksquare(B);

        // copy A to C to protect from overwriting
        copyvec(C.Store(),A.Store(),A.Elements());

        charT   side('L'), uplo('U'), transc('N'), diag('N');
        intT    m(C.Cols()), n(C.Rows()),
                ldb(B.Cols()), ldc(C.Cols());
        doubleT alpha(1.0);
        
        F77NAME(dtrmm)(&side,&uplo,&transc,&diag,&m,&n,
                       &alpha,B.Store(),&ldb, C.Store(),&ldc);


    } else if ( UT(B) ) {             // full matrix * upper triangle
#ifdef DEBUG
        cout << "GM*UT\n";
#endif
        checksquare(B);

        // copy A to C to protect from overwriting
        copyvec(C.Store(),A.Store(),A.Elements());

        charT   side('L'), uplo('L'), transc('N'), diag('N');
        intT    m(C.Cols()), n(C.Rows()),
                ldb(B.Cols()), ldc(C.Cols());
        doubleT alpha(1.0);
        
        F77NAME(dtrmm)(&side,&uplo,&transc,&diag,&m,&n,
                       &alpha,B.Store(),&ldb, C.Store(),&ldc);


    } else if ( LT(A) ) {            // lower triangle * full matrix
#ifdef DEBUG
        cout << "LT*GM\n";
#endif

        checksquare(A);

        // copy B to C to protect from overwriting
        copyvec(C.Store(),B.Store(),B.Elements());

        charT   side('R'), uplo('U'), transc('N'), diag('N');
        intT    m(C.Cols()), n(C.Rows()),
                ldb(A.Cols()), ldc(C.Cols());
        doubleT alpha(1.0);
        
        F77NAME(dtrmm)(&side,&uplo,&transc,&diag,&m,&n,
                       &alpha,A.Store(),&ldb, C.Store(),&ldc);



    } else if ( UT(A) ) {            // upper triangle * full matrix
#ifdef DEBUG
        cout << "UT*GM\n";
#endif
        checksquare(A);

        // copy A to C to protect from overwriting
        copyvec(C.Store(),B.Store(),B.Elements());

        charT   side('R'), uplo('L'), transc('N'), diag('N');
        intT    m(C.Cols()), n(C.Rows()),
                ldb(A.Cols()), ldc(C.Cols());
        doubleT alpha(1.0);
        
        F77NAME(dtrmm)(&side,&uplo,&transc,&diag,&m,&n,
                       &alpha,A.Store(),&ldb, C.Store(),&ldc);

    } else /* GM(A) and GM(B) */ {   // GM*GM: full matrix * full matrix
#ifdef DEBUG
        cout << "GM*GM\n";
#endif

        charT   t('N');
        intT    m(B.Cols()), n(A.Rows()), k(B.Rows()),
                lda(A.Cols()), ldb(B.Cols()), ldc(C.Cols());
        doubleT alpha(1.0), beta(0.0);
        
        F77NAME(dgemm)(&t,&t, &m,&n,&k,
                       &alpha,B.Store(),&ldb, A.Store(),&lda, 
//.........这里部分代码省略.........

void MpSurfaceBands (Scene& S, const Matrix& z, 
                     int direction, double width, 
                     const ColorMap& cmap, double cmin, double cmax)
{
  Vertex *v0,*v1,*v2,*v3;
  double scale = 0.0;
  bool usecolormap;

  // the width of the bands
  MpForceRange(width,1e-3,1.0);

  // 0.5 is half the mesh unit
  width /= 2;
    
  int ncmap = 0,
    nx = z.Rows(), 
    ny = z.Cols();

  // check arguments
  if ( nx < 0 || ny < 0 ) return; // nothing to draw
   
  // check if a colormap should be used
  if (cmap) { 
    usecolormap = true;
    ncmap = cmap.Size()-1;
    double dz = cmax - cmin; 
    scale = (dz == 0.0) ? 0.0 : ncmap/dz;
  } else
    usecolormap = false;
        
  switch (direction) {

    case RowBandsGrid: 
    case RowBands:
        {

      for ( int y = z.Clo(); y <= z.Chi(); y++ ) { 
        int x = z.Rlo();            
        v0 = &S.AddVertex(y+width,z[x][y],x);
        v1 = &S.AddVertex(y-width,z[x][y],x);
        if (usecolormap) { 
          int c = int((z[x][y]-cmin)*scale);
          if (c > ncmap) 
            c = ncmap;
          else if (c < 0) 
            c = 0;
          const ColorB &col = cmap[c];
          v0->SetColor(col);
          v1->SetColor(col);
        }

        for ( ; x < z.Rhi(); x++) {                 
          // create facet
          Facet &f = S.AddFacet(4);
          v2 = &S.AddVertex(y-width,z[x+1][y],x+1);
          v3 = &S.AddVertex(y+width,z[x+1][y],x+1);
          f(0) = v0;
          f(1) = v1;
          f(2) = v2;
          f(3) = v3;
          v0 = v3;
          v1 = v2;

          // set edge line mask
          if (direction == RowBands)
            if (x == z.Rlo())
              f.SetEdgeLineVisibility(Facet::Individual,0xb); // (1011=0xb)
            else if (x == z.Rhi()-1)
              f.SetEdgeLineVisibility(Facet::Individual,0xe); // (1110=0xe)
            else
              f.SetEdgeLineVisibility(Facet::Individual,0xa); // (1010=0xa)

          // set color of facet according to colormap
          if (usecolormap) { 
            float cx,cy,cz;
            f.Center(cy,cz,cx);
            int c = int((cz-cmin)*scale);
            if (c > ncmap) 
              c = ncmap;
            else if (c < 0) 
              c = 0;
            const ColorB &center = cmap[c];
            f.SetColor(center);
            c = int((z[x+1][y]-cmin)*scale);
            if (c > ncmap) 
              c = ncmap;
            else if (c < 0) 
              c = 0;
            const ColorB &col = cmap[c];
            v2->SetColor(col);
            v3->SetColor(col);
          }
        }
                
      }
      break;
        }

    case ColumnBandsGrid:
    case ColumnBands:
//.........这里部分代码省略.........