Modifying a Graphics3D object generated by ParametricPlot3D

Question

Here is a set of structured 3D points. Now we can form a BSpline using these points as knots.

dat=Import["3DFoil.mat", "Data"]
fu=BSplineFunction[dat]

Here we can do a ParametricPlot3D with these points.

pic=ParametricPlot3D[fu[u,v],{u, 0, 1}, {v, 0, 1}, Mesh -> All, AspectRatio -> 
Automatic,PlotPoints->10,Boxed-> False,Axes-> False]

Question

If we carefully look at the 3D geometry coming out of the spline we can see that it is a hollow structure. This hole appears in both side of the symmetric profile. How can we perfectly (not visually!) fill up this hole and create a unified Graphics3D object where holes in both sides are patched.

What I am able to get so far is the following. Holes are not fully patched.

I am asking too many questions recently and I am sorry for that. But if any of you get interested I hope you will help.

Update

Here is the problem with belisarius method. It generates triangles with almost negligible areas.

dat = Import[NotebookDirectory[] <> "/3DFoil.mat", "Data"];
(*With your points in "dat"*)
fd = First@Dimensions@dat;
check = ParametricPlot3D[{BSplineFunction[dat][u, v], 
BSplineFunction[{dat[[1]], Reverse@dat[[1]]}][u, v], 
BSplineFunction[{dat[[fd]], Reverse@dat[[fd]]}][u, v]}, {u, 0, 
1}, {v, 0, 1}, Mesh -> All, AspectRatio -> Automatic, 
PlotPoints -> 10, Boxed -> False, Axes -> False]

output is here

Export[NotebookDirectory[]<>"myres.obj",check];
cd=Import[NotebookDirectory[]<>"myres.obj"];
middle=
check[[1]][[2]][[1]][[1(* Here are the numbers of different Graphics group*)]][[2,1,1,1]];
sidePatch1=check[[1]][[2]][[1]][[2]][[2,1,1,1]];
sidePatch2=check[[1]][[2]][[1]][[3]][[2,1,1,1]];

There are three Graphics groups rest are empty. Now lets see the area of the triangles in those groups.

polygonArea[pts_List?
(Length[#]==3&)]:=Norm[Cross[pts[[2]]-pts[[1]],pts[[3]]-pts[[1]]]]/2
TriangleMaker[{a_,b_,c_}]:={vertices[[a]],vertices[[b]],vertices[[c]]}
tring=Map[polygonArea[TriangleMaker[#]]&,middle];
tring//Min

For the middle large group output is

0.000228007

This is therefore a permissible triangulation. But for the side patches we get zero areas.

Map[polygonArea[TriangleMaker[#]] &, sidePatch1] // Min
Map[polygonArea[TriangleMaker[#]] &, sidePatch2] // Min

Any way out here belisarius ?

My partial solution

First download the package for simplifying complex polygon from Wolfram archive.

fu = BSplineFunction[dat];
pic =(*ParametricPlot3D[fu[u,v],{u,0,1},{v,0,1},Mesh->None,
AspectRatio->Automatic,PlotPoints->25,Boxed->False,Axes->False,
BoundaryStyle->Red]*)
ParametricPlot3D[fu[u, v], {u, 0, 1}, {v, 0, 1}, Mesh -> None, 
AspectRatio -> Automatic, PlotPoints -> 10, Boxed -> False, 
Axes -> False, BoundaryStyle -> Black];
bound = First@Cases[Normal[pic], Line[pts_] :> pts, Infinity];
corners = Flatten[Table[fu[u, v], {u, 0, 1}, {v, 0, 1}], 1];
nf = Nearest[bound -> Automatic]; {a1, a2} = 
Union@Flatten@(nf /@ corners);
sets = {bound[[2 ;; a1]], bound[[a1 ;; a2]],bound[[a2 ;; a2 + a1]]};
CorrectOneNodeNumber = Polygon[sets[[{1, 3}]]][[1]][[1]] // Length;
CorrectOneNodes1 = 
Polygon[sets[[{1, 3}]]][[1]][[1]]; CorrectOneNodes2 = 
Take[Polygon[sets[[{1, 3}]]][[1]][[2]], CorrectOneNodeNumber];
<< PolygonTriangulation`SimplePolygonTriangulation`
ver1 = CorrectOneNodes1;
ver2 = CorrectOneNodes2;
triang1 = SimplePolygonTriangulation3D[ver1];
triang2 = SimplePolygonTriangulation3D[ver2];
Show[Graphics3D[{PointSize[Large], Point[CorrectOneNodes1]},Boxed -> False,
BoxRatios -> 1], Graphics3D[{PointSize[Large], Point[CorrectOneNodes2]},
Boxed -> False, BoxRatios -> 1],
Graphics3D[GraphicsComplex[ver1, Polygon[triang1]], Boxed -> False,
BoxRatios -> 1],
Graphics3D[GraphicsComplex[ver2, Polygon[triang2]], Boxed -> False,
BoxRatios -> 1]]

We get nice triangles here.

picfin=ParametricPlot3D[fu[u,v],{u,0,1},  {v,0,1},Mesh->All,AspectRatio->Automatic,PlotPoints->10,Boxed->False,Axes->False,BoundaryStyle->None];pic3D=Show[Graphics3D[GraphicsComplex[ver1,Polygon[triang1]]],picfin,Graphics3D[GraphicsComplex[ver2,Polygon[triang2]]],Boxed->False,Axes->False]

Now this has just one problem. Here irrespective of the PlotPoints there are four triangles always appearing that just shares only one edge with any other neighboring triangle. But we expect all of the triangles to share at least two edges with other trangles. That happens if we use belisarius method. But it creates too small triangles that my panel solver rejects as tingles with zero area.

One can check here the problem of my method. Here we will use the method from the solution by Sjoerd.

Export[NotebookDirectory[]<>"myres.obj",pic3D];
cd=Import[NotebookDirectory[]<>"myres.obj"];
polygons=(cd[[1]][[2]]/.GraphicsComplex-> List)[[2]][[1]][[1,1]];
pt=(cd[[1]][[2]]/.GraphicsComplex-> List)[[1]];
vertices=pt;
(*Split every triangle in 3 edges,with nodes in each edge sorted*)
triangleEdges=(Sort/@Subsets[#,{2}])&/@polygons;
(*Generate a list of edges*)
singleEdges=Union[Flatten[triangleEdges,1]];
(*Define a function which,given an edge (node number list),returns the bordering*)
(*triangle numbers.It's done by working through each of the triangles' edges*)
ClearAll[edgesNeighbors]
edgesNeighbors[_]={};
MapIndexed[(edgesNeighbors[#1[[1]]]=Flatten[{edgesNeighbors[#1[[1]]],#2[[1]]}];
edgesNeighbors[#1[[2]]]=Flatten[{edgesNeighbors[#1[[2]]],#2[[1]]}];
edgesNeighbors[#1[[3]]]=Flatten[{edgesNeighbors[#1[[3]]],#2[[1]]}];)&,triangleEdges];

(*Build a triangle relation table.Each'1' indicates a triangle relation*)
relations=ConstantArray[0,{triangleEdges//Length,triangleEdges//Length}];
Scan[(n=edgesNeighbors[##];
If[Length[n]==2,{n1,n2}=n;
relations[[n1,n2]]=1;relations[[n2,n1]]=1];)&,singleEdges]
(*Build a neighborhood list*)
triangleNeigbours=Table[Flatten[Position[relations[[i]],1]],{i,triangleEdges//Length}];
trires=Table[Flatten[{polygons[[i]],triangleNeigbours[[i]]}],{i,1,Length@polygons}];
Cases[Cases[trires,x_:>Length[x]],4]

Output shows always there are four triangles that shares only one edges with others.

{4,4,4,4}

In case of belisarius method we don't see this happening but there we get triangles with numerically zero areas.

BR

Brutus

Take a look at reference.wolfram.com/mathematica/TetGenLink/tutorial/…. It is designed for this kind of things – Shamus

@belisarius I have used TetGenLink but it is not meant for this. It makes tetrahedral solid mesh. I need a surface mesh. I actually want to use the Mathematica generated surface mesh. But need to get a solid body defined by Graphics3D or GraphicsComplex. I have written code to do the rest. It works well. But here I cant patch the hole in the first place. Finally again TetGen is not a solution. – Brutus

How could you obtain a solution where all of the triangles to share at least two edges with other triangles if you have (for example) only four vertices? – Shamus

Answer 1

Your data set looks like this:

Graphics3D[Point@Flatten[dat, 1]]

It consists of 22 sections of 50 points.

Adding a mid-line in each end section (which is actually the end section flattened):

dat2 = Append[Prepend[dat, 
                      Table[(dat[[1, i]] + dat[[1, -i]])/2, {i, Length[dat[[1]]]}]
              ], 
              Table[(dat[[-1, i]] + dat[[-1, -i]])/2, {i, Length[dat[[-1]]]}]
       ];

Graphics3D[{Point@Flatten[dat, 1], Red, Point@dat2[[1]], Green, Point@dat2[[-1]]}]

Now add some weights to the wingtip rim:

sw = Table[1, {24}, {50}];
sw[[2]] = 1000 sw[[1]];
sw[[-2]] = 1000 sw[[1]];
fu = BSplineFunction[dat2, SplineWeights -> sw];

Show[
  ParametricPlot3D[fu[u, v], {u, 0, 1}, {v, 0, 1}, Mesh -> All, 
                      AspectRatio -> Automatic, PlotPoints -> 20, Boxed -> False, 
                      Axes -> False, Lighting -> "Neutral"
  ], 
  Graphics3D[{PointSize -> 0.025, Green, Point@dat2[[-1]], Red,Point@dat2[[-2]]}]
]

Note that I increased the PlotPoints value to 20.

Befriend

Answer 2

Import the data and construct the BSpline function as before:

dat = Import["Downloads/3DFoil.mat", "Data"];

fu = BSplineFunction[dat]

Generate the surface, making sure to include (only) the boundary line, which will follow the edge of the surface. Make sure to set Mesh to either All or None.

pic = ParametricPlot3D[fu[u, v], {u, 0, 1}, {v, 0, 1}, Mesh -> None, 
  AspectRatio -> Automatic, PlotPoints -> 10, Boxed -> False, 
  Axes -> False, BoundaryStyle -> Red]

Extract the points from the boundary line:

bound = First@Cases[Normal[pic], Line[pts_] :> pts, Infinity]

Find the "corners", based on your parameter space:

corners = Flatten[Table[fu[u, v], {u, 0, 1}, {v, 0, 1}], 1]

Find the edge points best corresponding to the corners, keeping in mind that ParametricPlot3D doesn't use the limits exactly, so we can't just use Position:

nf = Nearest[bound -> Automatic];
nf /@ corners

Figure our which range of points on the boundary correspond to the areas you need to fill up. This step involved some manual inspection.

sets = {bound[[2 ;; 22]], bound[[22 ;; 52]], bound[[52 ;; 72]], 
  bound[[72 ;;]]}

Construct new polygons corresponding to the holes:

Graphics3D[Polygon[sets[[{1, 3}]]], Boxed -> False, BoxRatios -> 1]

Show[pic, Graphics3D[Polygon[sets[[{1, 3}]]]]]

Note that there is probably still a hole that can't be seen where the edge runs between the holes you mentioned, and I haven't tried to fill it in, but you should have enough information to do that if needed.

Croteau

nf /@ corners returns {{22},{22},{52},{52}} and with what logic you can find sets = {bound[[2 ;; 22]], bound[[22 ;; 52]], bound[[52 ;; 72]], bound[[72 ;;]]} from that? I can not get this. Anyways it was very helpful still I need to get rid of the "manual inspection". I will try.. – Brutus

@plato in this case the ends all have roughly the same z-value, which might generalize for you. In this case I'd initially used ;;22, 22;;52 and 52;;, but it included too much. I threw out the first point, and made the third segment the same length based on the idea that it's probably symmetric (enough.) – Croteau

Answer 3

(*With your points in "dat"*)
fu = BSplineFunction[dat[[1 ;; 2]]];
Show[{ParametricPlot3D[fu[u, v], {u, 0, 1}, {v, 0, 1}, 
                      Mesh -> All, AspectRatio -> Automatic, PlotPoints -> 30], 
      ListPlot3D[dat[[1]]]}]

And with

InputForm[%]

you get the "unified" graphics object.

Edit

Another way, probably better:

(*With your points in "dat"*)
fu = BSplineFunction[dat];
Show[

{ ParametricPlot3D[fu[u, v], {u, 0, 1}, {v, 0, 1}, 
                       Mesh -> All, AspectRatio -> Automatic, 
                       PlotPoints -> 10, Boxed -> False, Axes -> False], 
  ParametricPlot3D[
   BSplineFunction[{First@dat, Reverse@First@dat}][u, v], {u, 0, 1}, {v, 0, 1},
                    Mesh -> None, PlotStyle -> Yellow], 
  ParametricPlot3D[
   BSplineFunction[{dat[[First@Dimensions@dat]],
                    Reverse@dat[[First@Dimensions@dat]]}]
                    [u, v], {u, 0, 1}, {v, 0, 1}]}]

In just one structure:

(*With your points in "dat"*)
fd = First@Dimensions@dat;
ParametricPlot3D[
 {BSplineFunction[dat][u, v],
  BSplineFunction[{dat[[1]],  Reverse@dat[[1]]}] [u, v],
  BSplineFunction[{dat[[fd]], Reverse@dat[[fd]]}][u, v]},
 {u, 0, 1}, {v, 0, 1},
 Mesh -> All, AspectRatio -> Automatic,
 PlotPoints -> 10, Boxed -> False, Axes -> False]

Edit

You can check that there are small triangles, but they are triangles indeed and not zero area polygons:

fu = BSplineFunction[dat];
check = ParametricPlot3D[{BSplineFunction[{First@dat, Reverse@dat[[1]]}][u, v]}, 
                         {u, 0, 1}, {v, 0, 1}, Mesh -> All, 
                         PlotStyle -> Yellow, Mesh -> All, AspectRatio -> Automatic, 
                         PlotPoints -> 10, Boxed -> False, Axes -> False];
pts = check /. Graphics3D[GraphicsComplex[a_, b__], ___] -> a;
m = check[[1]][[2]][[1]][[1]] /. {___, GraphicsGroup[{Polygon[a_]}]} -> a;
t = Replace[m, {a_, b_, c_} -> {pts[[a]], pts[[b]], pts[[c]]}, {1}];
polygonArea[pts_List?(Length[#] == 3 &)] := 
                                 Norm[Cross[pts[[2]] - pts[[1]], pts[[3]] - pts[[1]]]]/2;

t[[Position[Ordering[polygonArea /@ t], 1][[1]]]]

(*
->{{{-4.93236, 0.0989696, -2.91748}, 
    {-4.92674, 0.0990546, -2.91748}, 
    {-4.93456, 0.100181, -2.91748}}}
*)

Shamus

I updated my question. Have a look in case you have some idea with those triangles with zero areas. – Brutus