☰

渐近线：根据 3D 中的数据点绘制曲面

tag-icon • tag-icon diagrams 3d asymptote
渐近线：根据 3D 中的数据点绘制曲面

在下图中，您可以看到 25 条曲线。每条曲线都是根据文件中的点列表绘制的mypoints_i.dat。每个文件都包含完全相同数量的点。

在此处输入图片描述

你会如何绘制曲线所描述的表面？提出了一个解决方案在这个答案中使用三角测量技术，但这里不需要使用三角测量技术，因为所有曲线都是从相同数量的点绘制的。此外，

  • 表面的颜色应从蓝色变为红色（如下面曲线所示）；
  • 表面应尽可能光滑。

编辑我设法绘制了表面，但无法使用平滑运算符operator..，因此表面非常粗糙，如下所示。这是因为我逐个绘制了平面四边形。当然，我可以细化点，但这会导致文件很大，并且 3D 操作会耗费资源。

一定有一种方法可以顺利地绘制表面，对吗？

请注意，尽管并不明显，但表面实际上对应于上面的曲线（一半）。 在此处输入图片描述

表面的代码如下：

//////// SURFACE PART I ////////

int increment=1;
for (int i=1; i<51-increment; i=i+increment)
{
string filename1 = filebasename + string(i) + "PartI.dat";
string filename2 = filebasename + string(i+increment) + "PartI.dat";
file in1=input(filename1).line().csv();
file in2=input(filename2).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
real[][] a2=in2.dimension(0,0);
a1=transpose(a1);
a2=transpose(a2);

real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];

real[] x2=a2[0];
real[] y2=a2[1];
real[] z2=a2[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

int step=1;
for (int j=0; j<21-step; j=j+step) 
{
    path3 orbit1=graph(x,y,z,operator--);
    path3 orbit2=reverse(graph(x2,y2,z2));

    triple pointA1=(x[j],y[j],z[j]);
    triple pointA2=(x[j+step],y[j+step],z[j+step]);
    triple pointB1=(x2[j],y2[j],z2[j]);
    triple pointB2=(x2[j+step],y2[j+step],z2[j+step]);

    draw(surface(pointA1..pointA2--pointB2..pointB1--cycle),red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(1),currentlight);
    }
}

编辑2完整的独立示例，包括数据文件，可用这里。

答案1

我又想了想，并设法进行了相当大的简化。此解决方案和前一个解决方案都遵循用户符号 1 建议的相同原理：如果您仔细查看光滑表面绘图器的操作，它实际上只绘制 uv 空间中矩形网格上的那些点。因此，即使您只有这些点（例如，在 s 矩阵中给出triple），您也可以“欺骗”绘图器相信您有一个完整的函数——它永远不会知道差异，因为它只会在网格点处评估函数。我相信这也是 OG 的“简短解决方案”所做的。

我之前的大部分解决方案都是为了弄清楚如何组合以任意顺序给出的“矩形”块而设计的。但由于这里的点是按常规顺序给出的，所以这样做有点过头了。（我主要是这样写的，因为我已经写好了代码，而 OP 已经在自己的代码中将网格排列成矩形块。）

关于调色板的注释：draw当该方法传递一组笔时，它将按常规（可预测）顺序排列“矩形”块，并用第一支笔绘制第一个块，用第二支笔绘制第二个块，等等。因此，从技术上讲，这幅图中绘制的不是平滑渐变；但在这个特定示例中，它足够接近。相比之下，surface.colors()符号 1 和 OG 的答案中使用的方法对顶点而不是块进行着色（允许平滑渐变，尽管不是prc）并将此信息存储在surface对象本身中。

// Global Asymptote definitions can be put here.
import three;
import grid3;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);


real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
    // PRC TRUE
    settings.prc=true;
    settings.embed=true;
}
else {
    // RASTERIZE
    settings.outformat="png";
    settings.prc=false;
    settings.render=3;
}

/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
{
    string filename1 = filebasename + string(i) + "PartI.dat";
    file in1=input(filename1).line().csv();
    real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
    a1=transpose(a1);

    real[] x1=a1[0];
    real[] y1=a1[1];
    real[] z1=a1[3];

    pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
    pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

    path3 thepath1 = graph(x1, y1, z1, operator..);
    draw(thepath1,orbitpen,currentlight);
}



//////// SURFACE PART I ////////

surface smoothSurface(triple[][] points) {
  int nu = points.length - 1;
  if (nu <= 0) abort("Grid must have at least two rows to produce a surface.");
  int nv = points[0].length - 1;
  if (nv <= 0) abort("Grid must have at least two columns to produce a surface.");
  // Create a parametric function that is designed specifically for integer values.
  triple f(pair uv) { return points[floor(uv.x)][floor(uv.y)]; }
  // Now graph that parametric function:
  return surface(f, (0,0), (nu, nv),
                nu=nu, nv=nv,
                usplinetype=Spline, vsplinetype=Spline);
}

int increment=1;

material[] surfacecolors = new material[0];

triple[][] points = new triple[0][0];

for (int i=1; i<51; i=i+increment)
{
    string filename1 = filebasename + string(i) + "PartI.dat";
    file in1=input(filename1).line().csv();
    real[][] a1=in1.dimension(0,0); // 0 means up to the end of the file
    a1=transpose(a1);

    real[] x=a1[0];
    real[] y=a1[1];
    real[] z=a1[3];

    triple[] currentrow = new triple[0];
    int step = 1;
    for (int j = 0; j < 21; j += step) {
            currentrow.push((x[j],y[j],z[j]));

        // Remember the color in which this patch should be drawn.
        if (i < 51-increment && j < 21-step)
          surfacecolors.push(
             red*(1-sqrt(1-i/50)) + rgb(0,62/255,91/255)*(sqrt(1-i/50)) + opacity(1));
    }
    points.push(currentrow);
}


surface smoothsurface = smoothSurface(points);

draw(smoothsurface, surfacecolors);

////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

结果如下：

在此处输入图片描述

答案2

请寻找其他解决方案。抱歉来晚了，我还有其他任务……

由于数据组织得非常好，因此可以按照与surface parametrized渐近线程序相同的精神构建参数化的贝塞尔曲面片。事实上，由实函数在 上box(a,b)（仅在网格上，f从框到 R 定义）描述的表面是 Scilab 函数的改编。

对于由参数函数描述的曲面f，box(a,b)过程如下（f在 3D 中取值）：如果f取决于(t,s)变量，则在中进行点的Dt均匀细分，在 中进行点的均匀细分。网格上的值，，，，，值。使用前面的例程，我们构造了 3 个贝塞尔参数化曲面； ，，。最后我们得到了贝塞尔参数化曲面。这是一种非常简单的方法，在 2D 中众所周知，您可能会遇到一些不适当的行为（参见 De Boor）。如果有人有更好的参数化选择（而不是均匀的），请告诉我。tnuDssnvFfFxFyFzxyz(Fx,{1,...,nu},{1,...,nv})(Fy,{1,...,nu},{1,...,nv})(Fz,{1,...,nu},{1,...,nv})(Fx,Fy,Fz)

将其适应于您的情况并不是很困难。

// Global Asymptote definitions can be put here.
import three;
import grid3;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);


real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
  // PRC TRUE
  settings.prc=true;
  settings.embed=true;
}
else {
  // RASTERIZE
  settings.outformat="png";
  settings.prc=false;
  settings.render=3;
}


/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
  {
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

path3 thepath1 = graph(x1, y1, z1, operator..);
draw(thepath1,orbitpen,currentlight);
  }



//////// SURFACE PART I ////////


triple[][] v=new triple[50][21];


int increment=1;
for (int i=1; i<52-increment; i=i+increment)
  {
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);
real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];
int step=1;
for (int j=0; j<22-step; j=j+step) 
      {
        v[i-1][j]=(x[j],y[j],z[j]);

      }
  }

int nu=50-1;
int nv=21-1;
real[] ipt=sequence(nu+1);
real[] jpt=sequence(nv+1);
real[][] fx=new real[nu+1][nv+1];
real[][] fy=new real[nu+1][nv+1];
real[][] fz=new real[nu+1][nv+1];

splinetype[] usplinetype=Spline;
splinetype[] vsplinetype=Spline;

for(int i=0; i <nu+1; ++i) {
  real ui=i;
  real[] fxi=fx[i];
  real[] fyi=fy[i];
  real[] fzi=fz[i];
  for(int j=0; j < nv+1; ++j) {
    pair z=(ui,j);
    fxi[j]=v[i][j].x;
    fyi[j]=v[i][j].y;
    fzi[j]=v[i][j].z;
  }
}
real[][][] sx=bispline(fx,ipt,jpt);//,usplinetype[0],vsplinetype[0]);
real[][][] sy=bispline(fy,ipt,jpt);//,usplinetype[1],vsplinetype[1]);
real[][][] sz=bispline(fz,ipt,jpt);//,usplinetype[2],vsplinetype[2]);

surface s=surface(sx.length);
s.index=new int[nu][nv];
int k=-1;
for(int i=0; i < nu; ++i) {
  int[] indexi=s.index[i];
  for(int j=0; j < nv; ++j)
    indexi[j]=++k;
}
for(int k=0; k < sx.length; ++k) {
  triple[][] Q=new triple[4][];
  real[][] Px=sx[k];
  real[][] Py=sy[k];
  real[][] Pz=sz[k];
  for(int i=0; i < 4 ; ++i) {
    real[] Pxi=Px[i];
    real[] Pyi=Py[i];
    real[] Pzi=Pz[i];
    Q[i]=new triple[] {(Pxi[0],Pyi[0],Pzi[0]),
                       (Pxi[1],Pyi[1],Pzi[1]),
                       (Pxi[2],Pyi[2],Pzi[2]),
                       (Pxi[3],Pyi[3],Pzi[3])};
  }
  s.s[k]=patch(Q);
}

draw(s,red);
////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

代码没有优化。你获取图片。我不管理调色板。

在此处输入图片描述

可以得到一个非常短的解决方案（但可读性不强）。“贝塞尔参数化”过程（50 行）被 5 行取代。

// Global Asymptote definitions can be put here.
import three;
import grid3;
import palette;

usepackage("mathptmx");

// One can globally override the default toolbar settings here:
// settings.toolbar=true;

import graph3;

real xmin=-2, xmax=2;
real ymin=-2, ymax=1.5;
real zmin=-2.5, zmax=2.5;

limits((xmin,ymin,zmin),(xmax,ymax,zmax));

currentprojection=perspective(camera=(1.5,2,2.5));

unitsize(3cm,3cm,2cm);


real linewidth=1.1;
real linewidthprojections=.15;

string filebasename="./data/ForASY_mnl2ippGridBranche3T1_";

bool renderPRC = false;

if(renderPRC) {
  // PRC TRUE
  settings.prc=true;
  settings.embed=true;
}
else {
  // RASTERIZE
  settings.outformat="png";
  settings.prc=false;
  settings.render=3;
}


/////// ORBITS IN 3D, SMOOTH ////////
for (int i=1; i<=50; i=i+2)
  {
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen orbitpen=.7bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50));
pen projpen=.3bp+red*(1-sqrt(1-i/50))+rgb(0,62/255,91/255)*(sqrt(1-i/50))+opacity(0.3);

path3 thepath1 = graph(x1, y1, z1, operator..);
draw(thepath1,orbitpen,currentlight);
  }



//////// SURFACE PART I ////////


triple[][] v=new triple[50][21];


int increment=1;
for (int i=1; i<52-increment; i=i+increment)
  {
string filename1 = filebasename + string(i) + "PartI.dat";
file in1=input(filename1).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);
real[] x=a1[0];
real[] y=a1[1];
real[] z=a1[3];
int step=1;
for (int j=0; j<22-step; j=j+step) 
      {
        v[i-1][j]=(x[j],y[j],z[j]);

      }
  }
triple f (pair t) {
  return(v[round(t.x)][round(t.y)]);
}

surface ns=surface(f,(0,0),(49,20),49,20,Spline);
ns.colors(palette(ns.map(zpart),Rainbow()));
draw(ns,render(merge=true));

////// CONTOUR OPENING //////

string fileImpactData = "./data/impactData.dat";
file in1=input(fileImpactData).line().csv();
real[][] a1=in1.dimension(0,0); // 0 pour dire jusqu'à la fin du fichier
a1=transpose(a1);

real[] x1=a1[0];
real[] y1=a1[1];
real[] z1=a1[3];

pen contourpen=green+1.5bp;

draw(graph(x1,y1,z1,operator--),contourpen,currentlight);

////// PLANES ///////
pen bg=gray(0.9)+opacity(0.2);
draw(surface((xmax,ymin,zmin)--(xmax,ymin,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmin,ymax,zmin)--(xmin,ymax,zmax)--(xmin,ymin,zmax)--(xmin,ymin,zmin)--cycle),bg);
draw(surface((xmax,ymax,zmin)--(xmax,ymin,zmin)--(xmin,ymin,zmin)--(xmin,ymax,zmin)--cycle),bg);

////// GRID LINES ///////
pen gridpen=.2bp+gray(0.7);

grid3(XYgrid,Step=1,gridpen);
grid3(YXgrid,Step=.5,gridpen);
grid3(XZgrid,Step=1,gridpen);
grid3(ZXgrid,Step=2,gridpen);
grid3(YZgrid,Step=.5,gridpen);
grid3(ZYgrid,Step=2,gridpen);

// No-go zone
draw((xmax,1,zmin)--(xmin,1,zmin)--(xmin,1,zmax),black+1bp);

xaxis3(Label("$x_1$",MidPoint,align=Y-Z),Bounds(Both,Min),InTicks(Step=1),p=black);
yaxis3(Label("$x_2$",MidPoint,align=X-Z),Bounds(Both,Min),InTicks(Step=.5),p=black);
zaxis3(Label("$\dot x_2$",MidPoint,align=X-Y),Bounds(Both,Min),InTicks(Step=2),p=black);

带有调色板的类似图片。 在此处输入图片描述

答案3

考虑这个例子。这个结果能说明你所需要的吗？

正如您在以下代码中看到的，2 数组P存储坐标，并且s是由构造的表面P。然后s由笔数组着色。

import three;
import palette;

size(12cm);
currentprojection=orthographic(1,1,1.5);
currentlight=(1,0,1);

triple P00=-X-Y+0.5*Z, P03=-X+Y, P33=X+Y, P30=X-Y;

triple[][] P={
  {P00,P00+(-0.5,0.5,-1),P03+(0,-0.5,1),P03},
  {P00+(0.5,-0.5,-1),(-0.5,-0.5,0.5),(-0.5,0.5,-1.5),P03+(0.5,0,1)},
  {P30+(-0.5,0,1),(0.5,-0.5,-1.5),(0.5,0.5,1),P33+(-0.5,0,1)},
  {P30,P30+(0,0.5,1),P33+(0,-0.5,1),P33}
};

surface s=surface(patch(P));
s.colors(palette(s.map(zpart),Gradient(yellow,red)));
// s.colors(palette(s.map(zpart),Rainbow()));

draw(s);
draw(sequence(new path3(int i){
      return s.s[i].external();},s.s.length), bp+orange);

if(!is3D())
  shipout(bbox(Fill(lightgrey)));

看起来不错。然而 另一个例子它警告

这里我们确定顶点的颜色（顶点着色）。由于 PRC 输出格式不支持贝塞尔曲面的顶点着色，因此 PRC 面片使用四种顶点颜色的平均值进行着色。

因此，只要您想要 PRC，就与渐近线无关。

警告：在我的电脑上，OpenGL 很难处理数千个贝塞尔曲面。因此，请在展示之前进行大量测试。

相关内容