# Plot Surface from Curves and Shortest Distance

I have the following curves:

contour0 = {{0, 0, 0}, {75., 21.541481684759155,
0}, {76.68311055933344, 23.316889440666557,
0}, {76.78414599838258, 25., 0}, {87.7566066865461,
37.24339331345389, 0}, {88.40938986817645, 50.,
0}, {93.28219721271785, 56.71780278728215, 0}, {100, 100, 0}};
contour25 = {{0, 0, 25}, {50., 7.321705206770043,
25}, {63.07509587419139, 11.92490412580861,
25}, {74.45040471732287, 25.549595282677124, 25}, {75.,
26.63604811037152, 25}, {83.55789081275054, 41.44210918724946,
25}, {100, 100, 25}};
contour50 = {{0, 0, 50}, {50., 8.664514126213406,
50}, {59.961068635033, 15.038931364967,
50}, {69.90507899416409, 30.094921005835914, 50}, {100, 100,
50}};
contour75 = {{0, 0, 75}, {84.36885742716176, 90.63114257283824,
75}, {100, 100, 75}};
contour100 = {{0, 0, 100}, {22.758147078400903, 2.2418529215990968,
100}, {25., 2.869118873270885, 100}, {39.792724873953674,
10.207275126046326, 100}, {48.86643807518578, 25., 100}, {100,
100, 100}}


Now, if I plot the curves they look like this: How can I draw a surface from those curves?

Once the surface is drawn, is there then a way to plot the shortest path belonging to the surface and connecting points {0,0,0} and {100,100,100}?

• Do you want the surface smoothened? If not, move the y-coordinates and weave Polygons over them. – Nicholas G Nov 1 '16 at 10:36

You can get a pretty good surface representation using Interpolation.

First combining the data into a single list:

data = Join[contour0, contour25, contour50, contour75, contour100];


Then formatting it into something Interpolation needs for higher dimensions ({{x,y},z})

formatdata = Table[{{data[[i, 3]], data[[i, 1]]}, data[[i, 2]]}, {i, 1, Length@data}];


I should note here that given the form of your data a bit of rearranging is necessary to avoid issues related to non-uniqueness. Now constructing the Interpolation function:

fun = Interpolation[formatdata, InterpolationOrder -> 1];


Again accounting for the reordering of data required for Interpolation the resulting surface can be visualized with

Show[Plot3D[fun[z, y], {z, Min@data[[;; , 3]], Max@data[[;; , 3]]},
{y, Min@data[[;; , 2]], Max@data[[;; , 2]]}],
ListPointPlot3D[data[[;; , {3, 1, 2}]], PlotStyle -> PointSize[Large]]] Edit:

If you are really concerned with the shape of the surface of sparsely populated regions, like the contour @ 75, create more data points along each contour to drive the behavior.

contour0fun = Interpolation[contour0[[;; , {1, 2}]], InterpolationOrder -> 1];
contour25fun = Interpolation[contour25[[;; , {1, 2}]], InterpolationOrder -> 1];
contour50fun = Interpolation[contour50[[;; , {1, 2}]], InterpolationOrder -> 1];
contour75fun = Interpolation[contour75[[;; , {1, 2}]], InterpolationOrder -> 1];
contour100fun = Interpolation[contour100[[;; , {1, 2}]], InterpolationOrder -> 1];

contour0 = Table[{{0, x}, contour0fun[x]}, {x, 0, 100, .1}];
contour25 = Table[{{25, x}, contour25fun[x]}, {x, 0, 100, .1}];
contour50 = Table[{{50, x}, contour50fun[x]}, {x, 0, 100, .1}];
contour75 = Table[{{75, x}, contour75fun[x]}, {x, 0, 100, .1}];
contour100 = Table[{{100, x}, contour100fun[x]}, {x, 0, 100, .1}];

formatdata = Join[contour0, contour25, contour50, contour75, contour100];

fun = Interpolation[formatdata, InterpolationOrder -> 1];

Show[Plot3D[fun[z, y], {z, Min@data[[;; , 3]], Max@data[[;; , 3]]},
{y, Min@data[[;; , 2]], Max@data[[;; , 2]]}],
ListPointPlot3D[data[[;; , {3, 1, 2}]], PlotStyle -> PointSize[Large]]] EDIT 2:

# Finding the shortest path between two points on a surface

We have now gone from a visualization problem to one of differential geometry. I should say right off the bat that I am not super proficent in this area, but I gave it my best shot. I think you are looking for something called a Geodesic. Having known nothing about this topic prior to your question I went straight to the Google to do some digging. There has been quite a lot of work in the area, and great summary can be found here. To the best of my understanding this is still very much an open problem in mathematics, particularly for mesh based surfaces.

Now, on to this this particular problem. In this case we have a $\left( u, v \right)$ parameterized surface defined as $\left<u,v,f\left(u,v\right)\right>$. With this type of form, the geodestics of the surface may be found through minimizing the arc length functional. For our particular surface definition, with $z=f\left(u,v\right)$, there exists a system of differential equations that can be solve to determine a geodestic (the full form of one of these equations can be found here).

I started by looking an example problem: What is the shortest distance between $\left(0,0\right)$ and $\left(10,10\right)$ on $\left<u,v,u\left(v-1\right)^2\right>$?

We can start by looking at the surface:

Plot3D[u*(v - 1)^2, {u, 0, 10}, {v, 0, 10}, ImageSize -> 400, AxesLabel -> Automatic] Now we need to setup up our geodesic calucation. Practically, this means setting up a system of coupled differential equations. For this particular example our system looks

sol2=NDSolve[{(x^\[Prime]\[Prime])[t]==-(((-1+y[t])^2 (4 (-1+y[t])
(x^\[Prime])[t] (y^\[Prime])[t]+2 x[t] (y^\[Prime])[t]^2))/(1+4 x[t]^2
(-1+y[t])^2+(-1+y[t])^4)),(y^\[Prime]\[Prime])[t]==-((2 x[t] (-1+y[t]) (4
(-1+y[t]) (x^\[Prime])[t] (y^\[Prime])[t]+2 x[t] (y^\[Prime])[t]^2))/
(1+4 x[t]^2 (-1+y[t])^2+(-1+y[t])^4)),x==0,y==0,x==10,y==10},
{x[t],y[t]},t]


Note that in the system we are looking for parametric equations x[t] and y[t] with $t \in \left[0,1\right]$. Solving this system gives the shortest path between our two points. We can extract an approximate representation of the path with

points2=Table[{sol2[[1,1,2]],sol2[[1,2,2]],sol2[[1,1,2]]*(sol2[[1,2,2]]-1)^2},
{t,0,1,.001}];


and visualize the path on the surface with

Show[Plot3D[x*(y-1)^2,{x,0,10},{y,0,10},ImageSize->400],
ListPointPlot3D[points2]/.Point[a___]:>{Thick,Line[a]}] You will notice, in general, the shortest path between two points on a surface is not necessarily a straight line, like in this example. Now, moving onto the question of the shortest distance between $\left\{0,0,0\right\}$ and $\left\{100,100,100\right\}$ on the surface defined by the 5 contours we need to reformulate our geodesic equations such that the surface descriptions are consistent. I ran into some issues in NDSolve with respect to extrapolation on the surface (there may be a work around with an application of WhenEvent, but I didn't look too closely into it), so I expanded the domain. In particular, we have

contour0fun=Interpolation[contour0[[;;,{1,2}]],InterpolationOrder->1];
contour25fun=Interpolation[contour25[[;;,{1,2}]],InterpolationOrder->1];
contour50fun=Interpolation[contour50[[;;,{1,2}]],InterpolationOrder->1];
contour75fun=Interpolation[contour75[[;;,{1,2}]],InterpolationOrder->1];
contour100fun=Interpolation[contour100[[;;,{1,2}]],InterpolationOrder->1];

contour10=Table[{{-10,x},contour0fun[x]},{x,0,100,.01}];
contour0=Table[{{0,x},contour0fun[x]},{x,0,100,.01}];
contour25=Table[{{25,x},contour25fun[x]},{x,0,100,.01}];
contour50=Table[{{50,x},contour50fun[x]},{x,0,100,.01}];
contour75=Table[{{75,x},contour75fun[x]},{x,0,100,.01}];
contour100=Table[{{100,x},contour100fun[x]},{x,0,100,.01}];
contour110=Table[{{110,x},contour100fun[x]},{x,0,100,.01}];

bound1=Table[{{i,-10},0},{i,-10.1,110,.0033}];
bound2=Table[{{i,110},100},{i,-10.1,110.1,.0033}];
bound3=Table[{{i,0},0},{i,-10.1,110.1,.0033}];
bound4=Table[{{i,100},100},{i,-10.1,110.1,.0033}];


Combining the data and creating a new expanded surface

formatdata2 = Join[contour10, contour0, contour25, contour50, contour75,
contour100, contour110, bound1, bound2, bound3, bound4];

fun2 = Interpolation[formatdata2, InterpolationOrder -> 1];


Because our surface is an interpolation object we have to take care in constructing the derivatives in the geodesic equations found in NDSolve. Making each derivative its own function

dx[x1_, y1_] := Evaluate[D[fun2[x, y], x]] /. x -> x1 /. y -> y1
dy[x1_, y1_] := Evaluate[D[fun2[x, y], y]] /. x -> x1 /. y -> y1
dxdx[x1_, y1_] := Evaluate[D[fun2[x, y], {x, 2}]] /. x -> x1 /. y -> y1
dydy[x1_, y1_] := Evaluate[D[fun2[x, y], {y, 2}]] /. x -> x1 /. y -> y1
dxdy[x1_, y1_] := Evaluate[D[D[fun2[x, y], x], y]] /. x -> x1 /. y -> y1


Now we can solve our system of equations directly using numeric derivatives of our surface

shortpath=NDSolve[{x''[t]==-(dx[x[t],y[t]]*(y'[t]^2*dydy[x[t],y[t]]
+2*x'[t]*y'[t]*dxdy[x[t],y[t]]+x'[t]^2*dxdx[x[t],y[t]]))
/(1+dx[x[t],y[t]]^2+dy[x[t],y[t]]^2),
y''[t]==-(dy[x[t],y[t]]*(x'[t]^2*dxdx[x[t],y[t]]
+2*y'[t]*x'[t]*dxdy[x[t],y[t]]+y'[t]^2*dydy[x[t],y[t]]))
/(1+dy[x[t],y[t]]^2+dx[x[t],y[t]]^2),x==0,y==0,x==100,y==100},
{x[t],y[t]},{t,0,1}];


We reconstuct our path with

shortpathpoints=Table[{shortpath[[1,1,2]],shortpath[[1,2,2]],
fun2[shortpath[[1,1,2]],shortpath[[1,2,2]]]},{t,0,1,.0001}];


It is possible to calculate the distance of the path by summing the individual line lengths.

Total@(EuclideanDistance @@@Partition[shortpathpoints,2,1])
(*192.702*)


Note that this is exactly the same distance as the straight line path between the two points

Total@(EuclideanDistance @@@Partition[Table[{y,y,fun2[y,y]},{y,0,100,.01}],2,1])
(*192.702*)


Visualizing

Show[Plot3D[fun2[z,y],{z,0,100},{y,0,100}],
ListPointPlot3D[shortpathpoints]/.Point[a___]:>{Thick,Line[a]},
ListPointPlot3D[Table[{y,y,fun2[y,y]},{y,0,100,.01}]]/.Point[a___]:>{Thick,Red,Line[a]}] This shortest path is entirely dependent on the description of the surface though. If the shapes of the contours or the interpolation between then change, then the shortest path is also likely to change.

• I do not know if the syntax with the unmatched square brackets after Show works in version 11 or if it is a typo, but in version 10 you'd Show[{ Plot3D[fun[z, y], {z, Min@data[[;; , 3]], Max@data[[;; , 3]]}, {y, Min@data[[;; , 2]], Max@data[[;; , 2]]}], ListPointPlot3D[data[[;; , {3, 1, 2}]], PlotStyle -> PointSize[Large]] }] – Nicholas G Nov 1 '16 at 15:16
• Definitely a typo. I think it is fixed now – Marchi Nov 1 '16 at 16:19
• Exactly what I was looking for, thank you! I have added a part to the question you might also be able to answer... – GEF Nov 1 '16 at 19:53
• Great stuff on shortest path. I vaguely recalled path of steepest ascent and planned to look something up tonight but your work above have saved me. +1 – Edmund Nov 2 '16 at 23:36
• @GEF To the best of my recollection, the expanded surface definition was to account for NDSolve guessing points outside the original domain. There might be a more elegant solution within NDSolve, but supplying the solver with more points was the quickest, workable solution I thought of at the time. – Marchi Jul 31 '17 at 15:34

You may use ListPlot3D and Union.

From your image it appears that the points are defined as $\{y,z,x\}$. Therefore I rearrange them to $\{x,y,z\}$ by Mapping (/@) #[[{3, 1, 2}]] &. InterpolationOrder can be adjusted to your liking.

ListPlot3D[#[[{3, 1, 2}]] & /@
Union @@ {contour0, contour25, contour50, contour75, contour100},
InterpolationOrder -> 2, AspectRatio -> 1] You can display your basis contours by setting mesh options; MeshFunctions, Mesh, and MeshStyle.

ListPlot3D[#[[{3, 1, 2}]] & /@
Union @@ {contour0, contour25, contour50, contour75, contour100},
InterpolationOrder -> 2, AspectRatio -> 1,
MeshFunctions -> {#1 &},
Mesh -> {Range[0, 100, 25]},
MeshStyle -> Directive[{Thick, Blue}]] Hope this helps.

• Thank you @Edmund ! Is there then a way to find the shortest distance on the surface connecting points {0,0,0} and {100,100,100}? – GEF Nov 2 '16 at 11:04