There are several possible interpretations of your question. The first is that you want to draw a surface representing the "outer shell" (convex hull? interpolation by a smooth surface?) of some data points given in no particular order. Mathematical algorithms for producing such surfaces virtually always produce triangulated surfaces. Since you insist on a rectangular mesh, I am not at all surprised that you were unable to find anything that does what you want.
However, the points in the data file you supply are not ordered randomly. Instead, as realized by alfC, the first consecutive four give the corners of a quadrilateral; the next consecutive four give the corners of a second quadrilateral; and so on. If you interpret the data file as giving a collection of quadrilaterals in no particular order, and the goal is to draw the surface that is their union, then this task is much easier. User alfC has already shown how it may be done using pgfplots; here is an Asymptote version, with comments to explain what is being done:
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));
Here's the result:
I also produced an alternative, that is designed to add several additional features:
- The output is a vector graphic rather than a rasterized one.
- The surface shown is a smooth surface.
- The density of the mesh can be adjusted, and need not be based on the number of points actually given.
The second criterion, in particular, requires a lot of extra programming, since I am assuming that the quadrilaterals are given in no particular order. Essentially, I have to reconstruct that order, and then tell Asymptote to use a Spline interpolation to get a (mostly) smooth surface out of it.
Here's the code:
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]));
http://asymptote.sourceforge.net/instead. – pluton Mar 7 '14 at 14:47