☰

pgfplots：具有自由格式的数据的 3D 曲面图

tag-icon • tag-icon pgfplots external-files 3d gnuplot octave
pgfplots：具有自由格式的数据的 3D 曲面图

我想从文件中绘制 3D 数据（数据文件） 使用pgfplots。数据表示一个椭圆体。我想像下图这样绘制外壳，只是没有背景，并且是灰度的：

椭圆体

我尝试了几种方法这里和这里但我没有运气。要么我没有得到椭圆体的外壳四边形，要么 gnuplot 在 100% CPU 负载下运行了几分钟却没有任何结果。

是否可以使用pgfplots和格努普特或者GNU Octave？该怎么做？TeX 内存够用吗？

我的 MWE：

\documentclass{scrreprt}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[ngerman]{babel}
\usepackage{pgfplots}
\usepackage{tikz}
\usetikzlibrary{calc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\begin{figure}[htb]
\centering
\begin{tikzpicture}
\begin{axis}
\addplot3 [surf] gnuplot [raw gnuplot] {set dgrid3d 1152,1152 spline;splot 'criterion.txt';};
\end{axis}
\end{tikzpicture}
\end{figure}

\end{document}

答案1

您的数据似乎排列在一组矩形块中。因此关键是patch type=rectangle。使用opacity选项可以避免（大部分）当前限制产生的叠加问题（与轴有关）pgfplots。（要绝对控制 3D 对象和灯光，请改用 Asymptote，无论如何都能pgfplots让您达到 99% 的效果，并且比 有所改进gnuplot）。

\documentclass{article}
\usepackage{pgfplots}
\pgfplotsset{compat=1.9}
\begin{document}
\begin{tikzpicture}
\begin{axis}[colormap/greenyellow, view = {150}{20}, axis equal, axis line style={opacity=0.5}, axis lines=center, xlabel=\small{$\sigma_\parallel$}, ticks=none, ylabel=\small{$\sigma_\perp$}, zlabel=\small{T}, xtick={}]
\addplot3+[patch, patch type=rectangle, mark=none, opacity=0.5, faceted color=black] file {
criterion.txt
}; 
\end{axis}
\end{tikzpicture}
\end{document}

这将生成以下图表：

克里托

参考数据criterion.txt文件如下所示：

-1229.428   -137.007    0.0
-1214.681   -163.451    0.0
-1215.0 -159.764    10.003
-1229.428   -137.007    0.0
-1214.681   -163.451    0.0
-1175.463   -187.298    0.0
-1176.097   -179.989    19.834
-1215.0 -159.764    10.003
-1175.463   -187.298    0.0
-1112.445   -208.142    0.0
...

（总行数：1152）

答案2

您的问题有几种可能的解释。第一种解释是，您想绘制一个表示“外壳”（凸包？通过光滑表面插值？）的曲面，这些曲面是无特定顺序给出的某些数据点的。用于生成此类曲面的数学算法几乎总是会产生三角化表面。由于你坚持矩形的mesh，我一点也不惊讶你找不到任何你想要的东西。

但是，您提供的数据文件中的点并不是随机排序的。相反，正如 alfC 所认识到的，前四个连续的点给出一个四边形的角；接下来的四个连续的点给出第二个四边形的角；依此类推。如果您将数据文件解释为给出一组无特定顺序的四边形，并且目标是绘制它们的并集的表面，那么这项任务就容易得多。用户 alfC 已经展示了如何使用 来完成它pgfplots；这里是一个 Asymptote 版本，带有注释来解释正在做什么：

defaultpen(fontsize(10));
size(345.0pt,0);    //Set the width of the resulting image.


settings.outformat="png";
settings.render=16;
usepackage("lmodern");    //Vectorized fonts are easier to render in 3d
import three;            //For drawing 3d things.

// Set the camera angle. (These numbers were obtained by experimentation.)
currentprojection = orthographic(camera=(14,14,10));

//Input the data into a two-dimensional array of "real" numbers:
file datafile = input("criterion.txt");
real[][] data = datafile.dimension(0,3);
close(datafile);

surface ellipsoid;        // The surface we are building
surface ellipsoidFacing;    // The subset consisting of only those patches that face the camera.
triple[] currentpatch;    // The surface patch currently being built

/* There's always a bit of programming involved in translating from a file.
 * Iterate over all the rows (i.e., all the lines of the file):
  */
for (real[] row : data) {
    //Add the current row to the list of triples:
    currentpatch.push((row[0], row[1], row[2]));

    //If we've described an entire rectangular patch, then add it to the surface and start a new patch:
    if (currentpatch.length == 4) {
        patch toAdd = patch(currentpatch[0] -- currentpatch[1] -- currentpatch[2] -- currentpatch[3] -- cycle);
        ellipsoid.push(toAdd);
        // Transparent surfaces often look better if only the patches facing the camera are considered.
        if (dot(toAdd.normal(0.5,0.5), currentprojection.camera) >= 0)
            ellipsoidFacing.push(toAdd);
        currentpatch.delete();
    }
}

//Draw the ellipsoid we've just built:
draw(ellipsoidFacing, surfacepen = material(white + opacity(0.6), specularpen=black), meshpen=black + linewidth(0.2pt));

//Find appropriate values for the minimum and maximum of the axes:
triple min = 1.1*min(ellipsoid);
triple max = 1.1*max(ellipsoid);
//Further adjustments will be made based on actual experimentation.

//Create (but do not draw) the three axes:
path3 xaxis = (min.x, 0, 0) -- (max.x, 0, 0);
path3 yaxis = (0, min.y, 0)  -- (0, 1.5*max.y, 0);
path3 zaxis = (0, 0, 2*min.z) -- (0, 0, 2*max.z);

//Now, draw the axes, together with their labels:
draw(xaxis, arrow=Arrow3, L=Label("$\sigma_{\parallel}$", position=EndPoint));
draw(yaxis, arrow=Arrow3, L=Label("$\sigma_{\perp}$", position=EndPoint));
draw(zaxis, arrow=Arrow3, L=Label("$\tau_{\parallel \perp}$", position=EndPoint));

//Finally, find, draw, and label the intersection points:
triple[] temp = intersectionpoints(xaxis, ellipsoid, fuzz=.01);
dot(temp[0], L=Label("$R_{\parallel d}$", align=SE));
dot(temp[1], L=Label("$R_{\parallel z}$", align=NW));

temp = intersectionpoints(yaxis, ellipsoid, fuzz=.01);
dot(temp[0], L=Label("$R_{\perp d}$", align=3NW));
dot(temp[1], L=Label("$R_{\perp z}$", align=NE));

temp = intersectionpoints(zaxis, ellipsoid, fuzz=.01);
dot(temp[0], L=Label("$R_{\parallel \perp}$", align=2*SE));
dot(temp[1], L=Label("$R_{\parallel \perp}$", align=NE));

结果如下： 在此处输入图片描述

我还制作了一个替代方案，旨在添加几个附加功能：

  1. 输出是矢量图形而不是光栅化图形。
  2. 所示表面是光滑表面。
  3. 网格的密度可以调整，而不必基于实际给出的点数。

第二个标准尤其需要大量额外的编程，因为我假设四边形没有特定的顺序。本质上，我必须重建该顺序，然后告诉 Asymptote 使用样条插值来获得（大部分）平滑的表面。

在此处输入图片描述

代码如下：

settings.outformat="pdf";
settings.render=0;
settings.prc=false;
usepackage("lmodern");
size(20cm);
import graph3;

file datafile = input("criterion.txt");
real[][] data = datafile.dimension(0,3);
close(datafile);

typedef triple[] quadpatch;
triple[] topEdge(quadpatch p) { return p[1:3]; }
triple[] botEdge(quadpatch p) { return new triple[] {p[3], p[0]}; }
triple[] leftEdge(quadpatch p) { return p[0:2]; }
triple[] rightEdge(quadpatch p) { return p[2:4]; }
triple botleft(quadpatch p) { return p[0]; }
triple botright(quadpatch p) { return p[3]; }
triple topleft(quadpatch p) { return p[1]; }
triple topright(quadpatch p) { return p[2]; }

bool edgesMatch(triple[] a, triple[] b) {
  if (a.length != b.length) return false;
  b = reverse(b);
  for (int i = 0; i < a.length; ++i) {
    if (abs(a[i] - b[i]) > .0001) return false;
  }
  return true;
}

bool secondAbove(quadpatch a, quadpatch b) {
  return edgesMatch(topEdge(a), botEdge(b));
}
bool secondRight(quadpatch a, quadpatch b) {
  return edgesMatch(rightEdge(a), leftEdge(b));
}


quadpatch[][] matrix;
void addToMatrix(quadpatch p, int i, int j) {
  while (matrix.length - 1 < i)
    matrix.push(new quadpatch[]);
  quadpatch[] currentrow = matrix[i];
  if (currentrow.length - 1 < j)
    currentrow.append(new quadpatch[j - currentrow.length + 1]);
  currentrow[j] = p;
}

struct PatchInGrid {
  quadpatch p;
  PatchInGrid left = null;
  PatchInGrid right = null;
  PatchInGrid above = null;
  PatchInGrid below = null;
};

quadpatch operator cast(PatchInGrid pig) { return pig.p; }

PatchInGrid[] patches;

void addQuadPatch(quadpatch p) {
  assert(p.length == 4);

  PatchInGrid toAdd;
  toAdd.p = p;

  for (int i = patches.length - 1; i >= 0; --i) {
    PatchInGrid possibility = patches[i];
    if (possibility.above == null && toAdd.below == null && secondAbove(possibility, p)) {
      possibility.above = toAdd;
      toAdd.below = possibility;
    }
    if (possibility.below == null && toAdd.above == null && secondAbove(p, possibility)) {
      possibility.below = toAdd;
      toAdd.above = possibility;
    }
    if (possibility.left == null && toAdd.right == null && secondRight(p, possibility)) {
      possibility.left = toAdd;
      toAdd.right = possibility;
    }
    if (possibility.right == null && toAdd.left == null && secondRight(possibility, p)) {
      possibility.right = toAdd;
      toAdd.left = possibility;
    }
  }

  patches.push(toAdd);
}

triple[] temp;
for (real[] currentpoint : data) {
  temp.push((currentpoint[0], currentpoint[1], currentpoint[2]));
  if (temp.length == 4) {
    addQuadPatch(temp);
    temp = new triple[];
  }
}

/* Start at patches[0] and find the leftmost bottommost patch connected to it.
 */
bool leftrightcyclic = false;
bool updowncyclic = false;
PatchInGrid currentpatch = patches[0];
PatchInGrid firstpatch = currentpatch;
while (currentpatch.left != null) {
  currentpatch = currentpatch.left;
  if (currentpatch == firstpatch) {
    leftrightcyclic = true;
    break;
  }
}
firstpatch = currentpatch;
while (currentpatch.below != null) {
  currentpatch = currentpatch.below;
  if (currentpatch == firstpatch) {
    updowncyclic = true;
    break;
  }
}

firstpatch = currentpatch;
quadpatch[][] patchMatrix;
PatchInGrid currentbottompatch = currentpatch;
do {
  quadpatch[] currentStrip;
  currentpatch = currentbottompatch;
  PatchInGrid bottom = currentbottompatch;
  do {
    currentStrip.push(currentpatch);
    /*
      if (currentpatch.above == null) {
      currentData.push(topleft(currentpatch));
      break;
      }
      if (currentpatch.above == bottom) {
      currentData.cyclic = true;
      break;
      }
    */
    currentpatch = currentpatch.above;
  } while (currentpatch != null && currentpatch != bottom);

  patchMatrix.push(currentStrip);

  /*
    if (currentbottompatch.right == null) {
    currentData = new triple[];
    do {
    currentData.push(botright(currentpatch));
    if (currentpatch.above == null) {
    currentData.push(topright(currentpatch));
    break;
    }
    if (currentpatch.above == bottom) {
    currentData.cyclic = true;
    break;
    }
    currentpatch = currentpatch.above;
    } while (currentpatch != null && currentpatch != bottom);
    thepoints.push(currentData);
    break;
    }
  */

  if (currentbottompatch.right == firstpatch) {
    patchMatrix.cyclic = true;
    break;
  }

  currentbottompatch = currentbottompatch.right;
} while (currentbottompatch != null && currentbottompatch != firstpatch);

triple f(pair uv) {
  int u = floor(uv.x);
  int v = floor(uv.y);
  int du = 0, dv = 0;
  if (!patchMatrix.cyclic && u >= patchMatrix.length) {
    assert(u == patchMatrix.length);
    --u;
    du = 1;
  }
  if (!patchMatrix[0].cyclic && v >= patchMatrix[0].length) {
    assert(v == patchMatrix[0].length);
    --v;
    dv = 1;
  }

  quadpatch inquestion = patchMatrix[u][v];
  if (du == 0) {
    if (dv == 0) return botleft(inquestion);
    else return topleft(inquestion);
  } else {
    if (dv == 0) return botright(inquestion);
    else return topright(inquestion);
  }
}

int nu = patchMatrix.length;
int nv = patchMatrix[0].length;

surface tempEllipsoid = surface(f, (0,0), (nu, nv),
                nu=nu, nv=nv,
                usplinetype=Spline, vsplinetype=Spline);

triple g(pair uv) { return tempEllipsoid.point(uv.x, uv.y); }
surface ellipsoid = surface(g, (0,0), (nu,nv-.001), nu=25, nv=40,
                usplinetype=Spline, vsplinetype=Spline);


currentprojection = orthographic(camera=(14,14,10));




triple min = 1.1*min(tempEllipsoid);
triple max = 1.1*max(tempEllipsoid);

path3 xaxis = min.x*X -- max.x*X;
real[] xaxisIsectionTimes = transpose(intersections(xaxis, tempEllipsoid, fuzz=.01))[0];
path3 xaxisInFront = subpath(xaxis, 0, xaxisIsectionTimes[0]);
path3 xaxisBehind = subpath(xaxis, xaxisIsectionTimes[0], length(xaxis));

path3 yaxis = min.y*Y -- 1.5*max.y*Y;
real[] yaxisIsectionTimes = transpose(intersections(yaxis, tempEllipsoid, fuzz=.01))[0];
path3 yaxisInFront = subpath(yaxis, yaxisIsectionTimes[1], length(yaxis));
path3 yaxisBehind = subpath(yaxis, 0, yaxisIsectionTimes[1]);

path3 zaxis = scale3(2)*(min.z*Z -- max.z*Z);
real[] zaxisIsectionTimes = transpose(intersections(zaxis, tempEllipsoid, fuzz=.01))[0];
path3 zaxisInFront = subpath(zaxis, zaxisIsectionTimes[1], length(zaxis));
path3 zaxisBehind = subpath(zaxis, 0, zaxisIsectionTimes[1]);


draw(xaxisBehind, arrow=Arrow3, L=Label("$\sigma_{\parallel}$",position=EndPoint), p=linewidth(0.8pt));
dot(point(xaxis,xaxisIsectionTimes[1]), L=Label("$R_{\parallel z}$",align=NW));
draw(yaxisBehind, p=linewidth(0.8pt));
dot(point(yaxis,yaxisIsectionTimes[0]));
draw(zaxisBehind, p=linewidth(0.8pt));
dot(point(zaxis,zaxisIsectionTimes[0]));






surface newEllipsoid;
for (patch p : ellipsoid.s) {
  if (dot(p.normal(1/2,1/2), currentprojection.camera) <= 0) newEllipsoid.push(p);
}
ellipsoid = newEllipsoid;
draw(ellipsoid, surfacepen=lightgray+opacity(0.5), meshpen=gray(0.4)+linewidth(0.2pt));




draw(xaxisInFront);
dot(point(xaxis,xaxisIsectionTimes[0]), L=Label("$R_{\parallel d}$", align=SE));
draw(yaxisInFront, arrow=Arrow3, L=Label("$\sigma_{\perp}$",position=EndPoint));
dot(point(yaxis,yaxisIsectionTimes[1]));
draw(zaxisInFront, arrow=Arrow3, L=Label("$\tau_{\parallel \perp}$", position=EndPoint));
dot(point(zaxis, zaxisIsectionTimes[1]));

答案3

谢谢alfC，

在您的帮助下，经过进一步的研究，我找到了一个让我非常满意的解决方案：

    \documentclass{scrreprt}

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    \usepackage[utf8]{inputenc}
    \usepackage[T1]{fontenc}
    \usepackage[ngerman]{babel}
    \usepackage{pgfplots}
    \usepgfplotslibrary{patchplots}
    \usetikzlibrary{calc}

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    \begin{document}

    \begin{figure}[htb]
    \centering
    \pgfplotsset{failurecriterion/.style={%
    compat=1.10,
    colormap={whitered}{color(0cm)=(white); color(1cm)=(black!75!gray)},
    view = {150}{20},
    axis equal image,
    axis lines=center,
    xlabel={$\sigma_{\parallel}$},
    ylabel={$\sigma_{\perp}$},
    zlabel={$\tau_{\parallel\perp}$},
    every axis x label/.style={at={(axis cs:\pgfkeysvalueof{/pgfplots/xmax},0,0)},xshift=-1em},
    every axis y label/.style={at={(axis cs:0,\pgfkeysvalueof{/pgfplots/ymax},0)},xshift=2ex},
    every axis z label/.style={at={(axis cs:0,0,\pgfkeysvalueof{/pgfplots/zmax})},xshift=1em},
    xmin=-1250, xmax=1750,
    ymin=- 300, ymax= 550,
    zmin=- 200, zmax= 350,
    ticks=none,
    width=1.0\linewidth,
    clip mode=individual,
    }}
    \begin{tikzpicture}
    \begin{axis}[failurecriterion]
    % Festigkeiten
    \addplot3 [only marks, mark size=1pt] coordinates {(1500,0,0) (-1000,0,0) (0,-240,0) (0,0,-150)};
    \node [above  left                            ] at (axis cs: 1500,   0,   0) {$R_{\parallel z}$};
    \node [below right                            ] at (axis cs:-1000,   0,   0) {$R_{\parallel d}$};
    \node [below      , xshift=0.5em, yshift= -2ex] at (axis cs:    0, 180,   0) {$R_{\perp z}$};
    \node [above  left, xshift= -2em, yshift=1.0ex] at (axis cs:    0,-240,   0) {$R_{\perp d}$};
    \node [above  left                            ] at (axis cs:    0,   0, 150) {$R_{\parallel\perp}$};
    \node [below  left              , yshift= -1ex] at (axis cs:    0,   0,-150) {$R_{\parallel\perp}$};
    % Versagenskoerper
    \addplot3+[patch, mark=none, opacity=0.5, patch type=rectangle,z buffer=sort,patch refines=1,line width=0.25pt] file {criterion.txt};
    % Festigkeiten Vordergrund
    \addplot3 [only marks, mark size=1pt] coordinates {(0,180,0) (0,0,150)};
    \end{axis}
    \end{tikzpicture}
    \end{figure}

    \end{document}

结果如下：

在此处输入图片描述

唯一奇怪的是，图片不能按线宽缩放。

相关内容