9
$\begingroup$

I have the following question: I have a file that has structure:

x1 y1 z1 f1
x2 y2 z2 f2
...
xn yn zn fn

I can easily visualize it with Mathematica using ListContourPlot3D. But could you please tell me how I can plot contour plot for this surface? I mean with these data I have a set of surfaces corresponding to different isovalues (f) and I want to plot intersection between all these surfaces and some certain plane. I tried to Google but didn't get any results. Any help and suggestions are really appreciated. Thanks in advance!

$\endgroup$
6
  • $\begingroup$ Would you please include a sample of your data, or provide a function that generates data that can serve as an example? $\endgroup$
    – Mr.Wizard
    Commented Aug 10, 2012 at 21:47
  • $\begingroup$ Also, do you want a 2D contour plot, or a 3D curve through space? $\endgroup$
    – Mr.Wizard
    Commented Aug 10, 2012 at 22:00
  • 2
    $\begingroup$ The answers to this question might be useful to you. $\endgroup$ Commented Aug 11, 2012 at 3:53
  • $\begingroup$ Thanks for all your replies! Sorry for the delay, I had problems with an internet acces. But the problem is when I'm trying to use suggestions by halirutan I get empty box (and I need exactly this plot). I'd love to provide the sample of my data but it's big, it has 2744 lines and I don't know how to upload this file here (any suggestions?). Thanks! PS Original questions was asked by me. I do not know how to comment other answers (I don't see button add commment) $\endgroup$
    – Nikita
    Commented Aug 13, 2012 at 20:50
  • $\begingroup$ After merging your accounts you should now be able to comment. $\endgroup$ Commented Aug 13, 2012 at 22:45

3 Answers 3

23
$\begingroup$

Ok, lets give this a try. @Mr.Wizard already showed you, how you can use Interpolate to make a function from your discrete data and since you didn't provide some test-data, I just assume we are speaking of an isosurface of a function $f(x,y,z)=c$ which is defined in some box in $\mathbb{R}^3$.

For testing we use $$f(x,y,z) = x^3+y^2-z^2\;\;\mathrm{and}\;\; -2\leq x,y,z \leq 2$$ which accidently happens to be the first example of ContourPlot3D.

The idea behind the following is pretty easy: As you may know from school, there is a simple representation of a plane in 3d which uses a point vector $p_0$ and two direction vectors $v_1$ and $v_2$. Every point on this plane can be reached through the $(s,t)$ parametrization

$$p(s,t)=p_0+s\cdot v_1+t\cdot v_2$$

Please note that $p_0, p, v_1, v_2$ are vectors in 3d and $s,t$ are scalars. The other form which we will use only for illustration is called the normal form of a plane. It is give by

$$n\cdot (p-p_0)=0$$

where $n$ is the vector normal to the plane, which can easily be calculated with the cross-product by $v_1\times v_2$. Lets start by looking at our example. To draw the plane inside ContourPlot3D we use the normal form which is plane2:

f[{x_, y_, z_}] := x^3 + y^2 - z^2;
v1 = {1, 1, 0};
v2 = {0, 0, 1};
p0 = {0, 0, 0};
plane1 = p0 + s*v1 + t*v2;
plane2 = Cross[v1, v2].({x, y, z} - p0);
gr3d = ContourPlot3D[{f[{x, y, z}], plane2}, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
  Contours -> {0}, 
  ContourStyle -> {ColorData[22, 3], Directive[Opacity[0.5], ColorData[22, 4]]}]

Mathematica graphics

What we do now is, that we try to find the contour value (which is 0 here) of $f(x,y,z)$ for all points, that lie on our plane. This is like doing a normal ContourPlot because our plane is 2d (although placed in 3d space). Therefore, we use the 2d to 3d mapping of plane1

gr2d = ContourPlot[f[plane1], {s, -2, 2}, {t, -2, 2}, Contours -> {0}, 
 ContourShading -> None, ContourStyle -> {ColorData[22, 1], Thick}]

Mathematica graphics

Look at the intersection. It is exactly the loop we would have expected from the 3d illustration. Now you could object something like "ew.. but I really like a curve in 3d..". Again, the mapping from this 2d curve to 3d is given in the plane equation. You can simply extract the Line[..] directives from the above plot and transfer it back to 3d:

Show[{gr3d, 
  Graphics3D[{Red, Cases[Normal[gr2d], Line[__], Infinity] /. 
     Line[pts_] :> Tube[p0 + #1*v1 + #2*v2 & @@@ pts, .05]}]
}]

I extract the Lines with Cases and then use the exact same definition of plane1 as pure function to transform the pts.

Mathematica graphics

When I'm not completely wasted at 5:41 in the morning, than this approach should work for your interpolated data too.

Apply method on test-data

I uploaded your test-data to our Git-repository and therefore, the code below should work without downloading anything. The approach is the same as above but some small things have changed, since we work on interpolated data now. I'll explain only the differences.

First we import the data and since we have a long list of {x,y,z,f} values, we transform them to {{x,y,z},f} as required by the Interpolation function. I'm not using the interpolation-function directly. I wrap a kind of protection around it which tests whether a given {x,y,z} is numeric and lies inside the interpolation box. Otherwise I just return 0.

data = {Most[#], Last[#]} & /@ 
   Import["https://raw.github.com/stackmma/Attachments/master/data_9304_187.m"];
ip = Interpolation[data];
fip[{x_?NumericQ, y_?NumericQ, z_?NumericQ}] := 
 If[Apply[And, #2[[1]] < #1 < #2[[2]] & @@@ 
    Transpose[{{x, y, z}, First[ip]}]],
  ip[x, y, z], 0.0]

The code below is almost the same. I only adapted the plane that it goes through your interpolation box. Furthermore, if you inspect your data you see that the value run from 0 to 1.2. Therefore I'm plotting the 0.5 contour by subtracting 0.5 from the function value and using Contours -> {0}. Remember that when I would simply plot the 0.5 contour, it would draw me a different plane, since we use one combined ContourPlot3D call.

Furthermore, notice that I normalized the direction vectors of the plane. This makes it easier to adjust the 2d plot of the contour. The rest should be the same.

v1 = Normalize[{30, 30, 0}];
v2 = Normalize[{0, 0, 21}];
p0 = {26, 26, 17};
plane1 = p0 + s*v1 + t*v2;
plane2 = Cross[v1, v2].({x, y, z} - p0);
gr3d = ContourPlot3D[{fip[{x, y, z}] - 0.5, plane2}, {x, 27, 30}, {y, 
   27, 30}, {z, 17.3, 21}, Contours -> {0}, 
  ContourStyle -> {Directive[Opacity[.5], ColorData[22, 3]], 
    Directive[Opacity[.8], ColorData[22, 5]]}]

gr2d = ContourPlot[fip[plane1] - 0.5, {s, 2, 5}, {t, 1, 4}, 
   Contours -> {0}, ContourShading -> None, 
   ContourStyle -> {ColorData[22, 1], Thick}];
Show[{gr3d, 
  Graphics3D[{Red, 
    Cases[Normal[gr2d], Line[__], Infinity] /. 
     Line[pts_] :> Tube[p0 + #1*v1 + #2*v2 & @@@ pts, .05]}]}]

Mathematica graphics

As you can see, your sphere has a whole inside.

Mathematica graphics

$\endgroup$
6
  • $\begingroup$ Thanks for all your replies! Sorry for the delay, I had problems with an internet acces. But the problem is when I'm trying to use suggestions by halirutan I get empty box (and I need exactly this plot). I'd love to provide the sample of my data but it's big, it has 2744 lines and I don't know how to upload this file here (any suggestions?). Thanks! PS Original questions was asked by me. I do not know how to comment other answers (I don't see button add commment) $\endgroup$
    – Nikita
    Commented Aug 13, 2012 at 20:50
  • $\begingroup$ thanks for solving account issue! So my question if it is possible to load file on this forum or how I can upload my data which are quite big (app. 3000 lines) $\endgroup$
    – Nikita
    Commented Aug 13, 2012 at 22:51
  • $\begingroup$ @Nikita I asked in the chat and there seems to be no standard-upload place. Oleksandr suggested ge.tt. Before uploading a very large file, please read the FAQ and remember that the questions on this site should be helpful for everyone and should not be too localized! $\endgroup$
    – halirutan
    Commented Aug 14, 2012 at 1:14
  • $\begingroup$ here is a link (depositfiles.com/files/r0pckyaiu) to zip file that contains 2 absolutely identical data (one is in format: x y z f, and the other is in the same format but can be directly load to mathematica). Yes, I understand that questions shouldn't be too localized, but the data I'm working on are cube files and it can be useful for many people. Thanks! $\endgroup$
    – Nikita
    Commented Aug 14, 2012 at 14:31
  • $\begingroup$ @Nikita Here we go. See my update. $\endgroup$
    – halirutan
    Commented Aug 14, 2012 at 15:56
7
$\begingroup$

You can use the options MeshFunctions in combination with Mesh for this.

I'm borrowing Mr.Wizard's data here for a moment:

data = Flatten[Table[{x, y, z, x^2 + y^2 - z^2}, {x, -2, 2, 0.2}, 
        {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];

Suppose you want to plot the intersection of the contours of data with the plane x - y == 0, then you could do something like

ListContourPlot3D[data, Contours -> {0.5, 2}, 
 ContourStyle -> Opacity[0.3],
 BoundaryStyle -> Opacity[0.3],
 MeshFunctions -> {(#1 - #2) &}, Mesh -> {{0}},
 MeshStyle -> {Thick, Orange}]

Mathematica graphics

$\endgroup$
2
  • 1
    $\begingroup$ You removed my glorious infix! +1 for the method however. (It's nice to see you posting again.) $\endgroup$
    – Mr.Wizard
    Commented Aug 11, 2012 at 7:00
  • $\begingroup$ thanks for your reply! These tricks with meshfunctions are brilliant, but when I use ListContourPlot3D on my data (which are exactly in the same format) I have empty box. But when I apply ListContourPlot3D on data that have only f values (without x,y,z) Mathematica plots correct figures but I cannot use MeshFunction in this case (because there is no information about x,y,z values). Do you know what can solve this problem? $\endgroup$
    – Nikita
    Commented Aug 13, 2012 at 22:58
5
$\begingroup$

I'm not claiming this is a good method, I'm just getting some ink the page:

data = Table[{x, y, z, x^2 + y^2 - z^2}, {x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}] ~Flatten~ 2;

ListContourPlot3D[data, Contours -> {0.5, 2}, Mesh -> None]

Mathematica graphics

int = Interpolation[data];

ContourPlot3D[int[x, y, z], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
 Contours -> {0.5, 2}, RegionFunction -> (-0.02 < #2 - # < 0.02 &)]

Mathematica graphics

$\endgroup$
3
  • $\begingroup$ The second piece of codes using RegionFunction doesn't creat a pretty plot,which contains two outer curves with small intervals, and neither PlotPoints nor MaxRecursion could help to improve it. $\endgroup$ Commented Aug 11, 2012 at 7:12
  • $\begingroup$ @paradox2 hence: "I'm not claiming this is a good method ..." $\endgroup$
    – Mr.Wizard
    Commented Aug 11, 2012 at 7:13
  • $\begingroup$ Your answer is tricky.Smart guy^-^ $\endgroup$ Commented Aug 11, 2012 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.