# Plotting contours plot for f(x,y,z)=c

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!

-
Would you please include a sample of your data, or provide a function that generates data that can serve as an example? –  Mr.Wizard Aug 10 '12 at 21:47
Also, do you want a 2D contour plot, or a 3D curve through space? –  Mr.Wizard Aug 10 '12 at 22:00
The answers to this question might be useful to you. –  Ｊ. Ｍ. Aug 11 '12 at 3:53
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) –  Nikita Aug 13 '12 at 20:50
After merging your accounts you should now be able to comment. –  Sjoerd C. de Vries Aug 13 '12 at 22:45
show 1 more comment

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]]}]


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}]


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.

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]}]}]


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

-
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) –  Nikita Aug 13 '12 at 20:50
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) –  Nikita Aug 13 '12 at 22:51
@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! –  halirutan Aug 14 '12 at 1:14
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! –  Nikita Aug 14 '12 at 14:31
@Nikita Here we go. See my update. –  halirutan Aug 14 '12 at 15:56
show 1 more comment

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}]


-
You removed my glorious infix! +1 for the method however. (It's nice to see you posting again.) –  Mr.Wizard Aug 11 '12 at 7:00
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? –  Nikita Aug 13 '12 at 22:58

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]


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 &)]


-
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. –  withparadox2 Aug 11 '12 at 7:12
@paradox2 hence: "I'm not claiming this is a good method ..." –  Mr.Wizard Aug 11 '12 at 7:13