Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Here is cylinder data: file

We can mesh the surface, thanks to lalmei, Belisarius and Simon (link)

How to achieve uv mapping from a 3d straight and curved tube?

Working code:

  Needs["TetGenLink`"]
 file="http://pastebin.com/H9y9SqYy";
 dat = Import[file, "Table"];
 {pts, tetrahedra} = TetGenDelaunay[dat];
 csr[{aa_, bb_, cc_, dd_}] := 
 With[{a = aa - dd, b = bb - dd, c = cc - dd}, 
 Norm[a.a Cross[b, c] + b.b Cross[c, a] + 
 c.c Cross[a, b]]/(2 Norm[a.Cross[b, c]])]
 radii = csr[pts[[#]]] & /@ tetrahedra;
 alphashape[rmax_] := Pick[tetrahedra, radii, r_ /; r < rmax] 
faces[tetras_] := Flatten[tetras /. {a_, b_, c_, d_} :> {{a, b, c}, 
 {a, b, d}, {a, c,     d}, {b, c, d}}, 1]

 externalfaces[faces_] := Cases[Tally[Sort /@ faces], {face_, 1} :> face]
 polys = externalfaces@faces@alphashape[0.01];
 Graphics3D[GraphicsComplex[pts, Polygon@polys], Boxed -> False]

enter image description here

share|improve this question
    
The unwrapped 2D positions of the points are already the UV coordinates, aren't they? –  Rahul Jun 15 at 19:56
    
@RahulNarain yes but it is difficult to mesh those 2d points –  Being Human Jun 15 at 19:58
    
Is there a way to keep the points posted ? –  lalmei Jun 27 at 14:47
    
@lalmei is there any stack exchange repository, I can post it there. –  Being Human Jun 27 at 15:31
    
May be PasteBin? –  Alexey Popkov Jun 27 at 18:58

2 Answers 2

up vote 6 down vote accepted

Here's an unwrapping of the side of the cylinder, if that is what you are after. With the OP's polys and pts (CoordinateTransform taken from @Kuba's answer), we begin by selecting only the polygons on the side of the cylinder (sidepolys). (There's no noise in this data, at least as far as the end caps go, so the tolerance tol could be set to 0.)

pts2 = CoordinateTransform[
    "Cartesian" -> "Cylindrical", {#, #3, #2} & @@@ pts][[;; , {2, 3}]];

tol = 0.00001;
sidepolys = 
  Cases[polys, 
   p_ /; Max@pts2[[p]][[All, 2]] > Min[pts2[[All, 2]]] + tol && 
          Min@pts2[[p]][[All, 2]] < Max[pts2[[All, 2]]] - tol];

The x-coordinates of points in pts2 lie between -Pi and Pi.

Through[{Max, Min}[#]] & /@ Transpose@pts2
(*
  {{3.13762, -3.13777}, {0.1, 0.}}
*)

The boundaries x == -Pi and x == Pi correspond to cutting the cylinder in order to unfold it. Some polygons will cross the cut. We'll separate the ones that stay inside the unfolded cylinder from those that crossed the cut.

{polysinside, polyscrossing} =
  GatherBy[sidepolys, (Max[#] - Min[#] < Pi &@pts2[[#]][[All, 1]]) &];

Two of the vertices of a crossing polygon (triangle) will be near one edge (x == -Pi and x == Pi) and one will be near the other. We can tell which edge from the sign of the x-coordinate. We can use UnitStep to encode the signs and a simple rule for pick which of the three vertices needs to be moved by a distance of 2 Pi closer to its mates.

signs = UnitStep[pts2[[All, 1]]];

posRule = {{0, 0, 1} -> 3, {0, 1, 0} -> 2, {0, 1, 1} -> 
    1, {1, 0, 0} -> 1, {1, 0, 1} -> 2, {1, 1, 0} -> 3};

We can then calculate the new coordinates and the new indices for the fixed vertices. We first calculate which vertex in each polygon should be translated and store it in repl. From this we calculate the new points newpts2 to be added to the GraphicsComplex. Next, we use repl again to replace the old vertex index in the polygon with the index of the new vertex and store the fixed polygons in polysfixed.

repl = Replace[signs[[#]] & /@ polyscrossing, posRule, 1];
newpts2 = pts2[[#]] + Thread[{2 Pi - 4 Pi signs[[#]], 0}] &@
   MapThread[Part, {polyscrossing, repl}];

i = Length[pts2];   (* starting index *)
polysfixed = MapThread[ReplacePart[#1, #2 -> ++i] &, {polyscrossing, repl}];

We can add them to the OP's GraphicsComplex. Here I've colored the fixed polygons for purposes of illustration. One could just as well join them to the other ones.

Graphics[GraphicsComplex[
  Join[pts2, newpts2], {White, EdgeForm@Directive[Thin, Gray], 
   Polygon@polysinside, EdgeForm@Directive[Thin, Blue], Red, 
   Polygon@polysfixed}], AspectRatio -> 1/2, Frame -> True]

Mathematica graphics

I'm not sure what is sought in question 2. Perhaps this?

Graphics[GraphicsComplex[
  Join[pts2, newpts2], {Polygon@Join[polysinside, polysfixed]}, 
  VertexColors -> (ColorData["Rainbow"] /@ 
     Rescale[Join[pts2, newpts2][[All, 2]], {0, 0.1}])], 
 AspectRatio -> 1/2, Frame -> True]

Mathematica graphics

share|improve this answer

Here is the general answer for any shaped object of surface genus-0, though it can have holes as long as it's an outer boundary (maybe its more general and someone can correct me).

I will first describe the general UV mapping. This is usually done for a surface with a pre-chosen boundary, you need to choose which points are part of the boundary and give them values. In your case the boundary is the uncapped edges of the cylinder, and periodic line along the cylinder.

We have the 3-D positions $\vec{r_i}=\{v^1_i,v^2_i,v^3_i\}$ and we want the two-dimensional mapping $\tilde{r}_i=\{u_i,v_i\}$

This can be done in several ways in Mathematica. But let me describe what we are really trying to do. All the methods boil down to a type of optimization a problem.

For example:

Taking the mesh you created (without the caps), we could break a continous set of edges along a side of the cylinder (or along a line connecting both caps), and break it. Now we have a 2D mesh that we can simply apply a specific embedding of it using the GraphLayout option in Graph, you can then pick the coordinates from the graphics, and normalize it to you 2D-choice.

Taking this view, can we see this is all an optimization of positions of the mesh in 2D. Well we can go directly to the idea of SpringEmbemding from Graph. Here the main goal is assume the points have a "spring" that connects them along the edges of the mesh. This allows us to know how points are connected to each other, and doesn't let points bypass closer points while performing the UV map.

So taking the spring view, one embedding would be one that minimizes the energy of the mesh, given a set of boundary conditions. So we would have to minimize the following function:

$$ E= \frac{1}{2} \sum_{i=1}^{n+b} \sum_{j \in adj_i} T_{ij} || \tilde{r_i}-\tilde{r_j}|| $$

where $adj_{i}$ is the list of points adjacent to point $i$, this is given by the mesh, and $T_{ij}$ is a weight associated with the edge, generally it can simply be the adjacency matrix matrix value, (there is mathematical problem with it, I think it is called linear reproduction, but it tends to work) Sometimes you want to chose different weights depending on the properties you want to preserve during the mapping (relative areas of the polygons, angles between edges etc...)

So we want to find the $\tilde{r}$ that minimize this function. Here there is no information about the position of the points in 3-D, only which points are connected to each. This can be solved by

$$ \tilde{r}_i = \sum_{j\in adj_i} w_{ij} \tilde{r}_j $$ where $$w_{ij} = \frac{T_{ij}}{\sum_{k\in adj_i} T_{ik}}$$

This just a boundary problem, which we can solve using LinearSolve once you build the appropriate weight matrix.

If you have a pre-defined boundary, or say patch you want to parametrize then the solution is very straight forward. You can separate the boundary points $\tilde{r}_b=\{u_b,v_b\}$ from the interior points. You get two separate linear equations, one for each dimensions.

$$ \mathbf{A}\cdot \mathbf{u} = \mathbf{\bar{u}}\,\;,\; \mathbf{A}\cdot \mathbf{v} = \mathbf{\bar{v}}$$

where

$$A_{ij}= \delta_{ij}- w_{ij} $$ with $w_{ij} =0 $ for $i=j$ and if $ j \notin adj_i$ and $$ \bar{r}_i = \{\bar{u}_i,\bar{v}_i\}= \sum_{j\in adj_{i}, j\in bound} w_{ij} \tilde{r}_j $$

The main problem in your case, is that there is no unique way to pick the place where to place the periodic boundary. If your object already had cylindrical symmetry in 3D you can easy do this by picking a side as done another answer.

I will give a more concrete example with that twisted cylinder you posted before. Since the cylinder is done, and I am trying to give a more general answer and show you that it will work with mostly any shape. Here it is a bit non-trivial to pick a "side."

So we will keep the cylinder capped, and do a spherical mapping of the closed tube, we will get the spherical mapping plus one of the coordinates will give us a way to map a boundary along the tube. Maybe later I will rewrite it so that it gives you the cylindrical mapping, but that is for a later date.

First we obtain the mesh as you mentioned:

Needs["TetGenLink`"]
file="https://dl.dropboxusercontent.com/u/68983831/curved_pipe02.txt";
dat=Import[file,"Table"];
{pts,tetrahedra}=TetGenDelaunay[dat];
csr[{aa_,bb_,cc_,dd_}]:=With[{a=aa-dd,b=bb-dd,c=cc-dd},Norm[a.a Cross[b,c]+b.b Cross[c,a]+c.c Cross[a,b]]/(2 Norm[a.Cross[b,c]])]
radii=csr[pts[[#]]]&/@tetrahedra;
alphashape[rmax_]:=Pick[tetrahedra,radii,r_/;r<rmax]
faces[tetras_]:=Flatten[tetras/.{a_,b_,c_,d_}:>{{a,b,c},{a,b,d},{a,c,d},{b,c,d}},1]
 externalfaces[faces_]:=Cases[Tally[Sort/@faces],{face_,1}:>face]
 polys=externalfaces@faces@alphashape[.001];
 Graphics3D[GraphicsComplex[pts,Polygon@polys],Boxed->False]

enter image description here

Note that is the same as before, however we haven't meshed all the points from the original. Because how we created the non-convex surface, by trimming large polygons. There are some points inside the object that haven't been meshed properly. There are algorithms to add these points back to the mesh, but this is way over complicating it. I will assume you don't mind losing some resolution around the bends.

However this will require us to rewrite the point list and polygons:

 pointspresent = Sort@DeleteDuplicates@Flatten[polys];

 newpolys = 
  polys /. Rule @@@ (Thread[{#, Range[Length@#]}] &[Sort@DeleteDuplicates@Flatten[polys]])
 newpoints = 
   Table[pts[[pointspresent[[i]]]], {i, Range[Length@pointspresent]}];

and we will re-write the polygons into an adjancy list so that finding adjacent points is only a single call.

 newadjlist = 
     Table[{x, 
       DeleteCases[
         DeleteDuplicates[
            Flatten[Cases[newpolys, z_ /; MemberQ[z, x]], 1]], x]} , {x,Range@Length@newpoints}];

In order to create the periodic boundary condition without the mapping, we need one extra piece of information. The normal vectors pointing away from the surface (or in.) This will gives us local information to decide which points are to right or the left of boundary. To find the out direction we use the points in the tetrahedrals we threw away while meshing the surface.

We simply average the normal vectors of the all the adjunct faces. There may be more efficient ways to do this, but this works.

 newfaces[f_]:=Flatten[f/.{a_,b_,c_,d_}:>{{{a,b,c},d},{{a,b,d},c},{{a,c,d},b},{{b,c,d},a}},1];

 faceswithinpoints = 
   Cases[newfaces@alphashape[.001], z_ /; MemberQ[polys, Sort[z[[1]]]]];


 areavectors = 
    SortBy[Flatten[#, 1] &@
     Map[Function[{x},  
       Map[Function[{y}, 
             {y,If[#.(pts[[x[[2]]]] - pts[[y]]) >= 0, -#, #] &[ 
                Cross@@((pts[[#]] - pts[[y]]) & /@ 
                  DeleteCases[x[[1]], y])]}], x[[1]]]], faceswithinpoints], First];


 {thepoints, thenorms} = 
      Thread[({Mean[#[[1]]], Mean[#[[2]]]} & /@ (Thread /@ 
          GatherBy[areavectors, #[[1]] &]))];

   (* Super Quick refers to how fast i coded, not how efficient it is *)

 myquicknormals = #[[2]] & /@ 
   SortBy[Thread[{thepoints /.Rule @@@ (Thread[{#, Range[Length@#]}] &[
   Sort@DeleteDuplicates@Flatten[polys]]), thenorms}], First];

You can check these do indeed give you the normals

 dd = RandomChoice[Range@Length@newpoints];
 Show[Graphics3D[GraphicsComplex[newpoints, Polygon@newpolys], 
     Boxed -> False], 
     Graphics3D[{Red, Arrowheads[0.03], 
     Arrow[{newpoints[[dd]], (newpoints[[dd]] + (0.004*
     Normalize@(myquicknormals[[dd]])))}]}],
  Graphics3D[{Blue, Sphere[newpoints[[dd]], 0.001]}]]

enter image description here

Spherical Mapping

We compute the mapping by solving the linear equation of the points with specific boundary conditions of a sphere, by choosing a north and south pole. If we just wanted to map it to a sphere, we can pick any two points we want... any two. Here since the object is oriented we pick the max and min in that direction. Once we find the "latitudes" we follow the max gradient of the latitude to assign that point to the periodic boundary that connects north and south.

Then we "fix" the weight matrix along the periodic boundary so that points adjunct to the boundary on one side are separate from other by the appropriate amount, $2 \pi$ in our case, this where the normals come in.

Then is simply of using LinearSolve again for the longitudinal, with the poles having a defined longitude since it is ambiguous. One might be able to use LinearProgramming instead to make sure the values dont go beyond $-\pi$ to $\pi$ Here we simply use Mod.

Needs["Developer`"]
ClearAll[computeLatitudes];
computeLatitudes[points0_, adj0_, pointsnotpoles0_, north0_,south0_] := 
Module[{Amat, points = points0, adj = adj0, 
  pointsnotpoles = pointsnotpoles0, north = north0, south = south0, 
  b, res, StartT},
Amat =ToPackedArray[SparseArray[{}, {Length[points] - 2, Length[points] - 2}, 0.0]];

Scan[Function[{x},
 Amat[[Position[pointsnotpoles, x][[1, 1]], 
   Position[pointsnotpoles, x][[1, 1]]]] = N[Length[adj[[x]]]];
 Scan[Function[{y},
   Amat[[Position[pointsnotpoles, x][[1, 1]], 
      Position[pointsnotpoles, y][[1, 1]]]] = -1.; 
   ], Select[adj[[x]], (# != north && # != south) &]];
 ], pointsnotpoles];
 b = SparseArray[{}, Length[points] - 2, 0.0];
 Scan[b[[Position[pointsnotpoles, #][[1, 1]]]] +=3.141592653589793; &,adj[[south]]];
 res = LinearSolve[Amat, b, Method -> "Multifrontal"];
 Return[{Amat, res}];
 ];

 parametrize2[points0_, adj0_] := 
   Module[{points = points0, adj = adj0, north, south, pointsnotpoles,res, bl, Almat, previous, here, maximum, nextpos, prevpos, theta, phi, resl, StartT},
 StartT = SessionTime[];

 {south, north} = 
  Part[Sort[
   Range[Length[
     points]], ((points[[#1]][[3]] > 
       points[[#2]][[3]])) &], {1, -1}];
  pointsnotpoles =Select[Range[Length[points]], (# != north && # != south) &];
  {Almat, res} = computeLatitudes[points, adj, pointsnotpoles, north, south];   
  theta = SparseArray[Rule@@@Transpose[{pointsnotpoles, res}],Length[points]];
  theta[[north]] = 0.0;theta[[south]] = 3.141592653589793;

  Scan[Function[{x},
    Almat[[Position[pointsnotpoles, x][[1, 1]], 
    Position[pointsnotpoles, x][[1, 1]]]] -= 1.;
   ], Last@(Union[##] & @@@ ({adj[[#]] & /@ {north, south}}))];    
  bl = SparseArray[{}, {Length[points] - 2}, 0.0];
  previous = north;
  here = RandomChoice[adj[[north]]];
  maximum = 0.0;
 mydateline = {north};
 westofdateline = {};
 While[here != south,
   Scan[Function[{x},
   If[theta[[x]] > maximum,
    maximum = theta[[x]];
    nextpos = x;
   ];
  If[x == previous,
   prevpos = x;
   ];
  ], adj[[here]]];

 mydateline = Append[mydateline, nextpos];
 westofdateline = 
 Append[westofdateline, 
  If[Sign[Dot[points[[nextpos]] - points[[here]], 
      Cross[myquicknormals[[here]], 
       points[[prevpos]] - points[[here]]]]] > 0,
   Cases[adj[[here]], 
    x_ /; (Sign[
          Dot[points[[x]] - points[[here]], 
           Cross[myquicknormals[[here]], 
            points[[prevpos]] - points[[here]]]]] < 0 || 
        Sign[Dot[points[[x]] - points[[here]], 
           Cross[myquicknormals[[here]], 
            points[[nextpos]] - points[[here]]]]] >= 0) && 
      x != prevpos && x != nextpos && x != here],
   Cases[adj[[here]], 
    x_ /; (Sign[
          Dot[points[[x]] - points[[here]], 
           Cross[myquicknormals[[here]], 
            points[[prevpos]] - points[[here]]]]] < 0 && 
        Sign[Dot[points[[x]] - points[[here]], 
           Cross[myquicknormals[[here]], 
            points[[nextpos]] - points[[here]]]]] >= 0) && 
      x != prevpos && x != nextpos && x != here]
   ]];    
 Scan[Function[{tt},
  If[Length[Position[pointsnotpoles, tt]] != 0,
   bl[[Position[pointsnotpoles, tt][[1, 1]]]] += 2.0*3.1416;
   bl[[Position[pointsnotpoles, here][[1, 1]]]] += -2.0*3.1416;]
  ], If[
  Sign[Dot[points[[nextpos]] - points[[here]], 
     Cross[myquicknormals[[here]], 
      points[[prevpos]] - points[[here]]]]] > 0,
  Cases[adj[[here]], 
   x_ /; (Sign[
         Dot[points[[x]] - points[[here]], 
          Cross[myquicknormals[[here]], 
           points[[prevpos]] - points[[here]]]]] < 0 || 
       Sign[Dot[points[[x]] - points[[here]], 
          Cross[myquicknormals[[here]], 
           points[[nextpos]] - points[[here]]]]] >= 0) && 
     x != prevpos && x != nextpos && x != here],
  Cases[adj[[here]], 
   x_ /; (Sign[
         Dot[points[[x]] - points[[here]], 
          Cross[myquicknormals[[here]], 
           points[[prevpos]] - points[[here]]]]] < 0 && 
       Sign[Dot[points[[x]] - points[[here]], 
          Cross[myquicknormals[[here]], 
           points[[nextpos]] - points[[here]]]]] >= 0) && 
     x != prevpos && x != nextpos && x != here]
  ]];
 previous = here;
 here = nextpos;
 ];
resl = LinearSolve[Almat, bl, Method -> "Multifrontal"];
phi = SparseArray[Rule @@@ Transpose[{pointsnotpoles, resl}], 
 Length[adj]];
phi[[north]] = 0.0;
phi[[south]] = 0.0;
Return[Transpose[{phi, theta}]];
];

 ListPlot[
  Thread[{Mod[#[[1]], 2 Pi, -Pi], #[[2]]} &@ 
    Thread[parametrize2[newpoints, #[[2]] & /@ newadjlist]]], Frame -> True,Axes -> False];

enter image description here

To show you the solution is indeed spherical we map the points with the mesh on a sphere.

 spherepoints = {Cos[#[[1]]] Sin[#[[2]]], Sin[#[[1]]] Sin[#[[2]]], 
 Cos[#[[2]]]} & /@ parametrize2[newpoints, #[[2]] & /@ newadjlist];

 graphfad = 
  SparseArray[
   Flatten[MapIndexed[{Function[{x}, {First@#2, x} -> 1] /@ #1} &, 
     #[[2]] & /@ newadjlist]], Length[newpoints]];

 vertexrulz = MapIndexed[First@#2 -> #1 &, spherepoints];

 Show[Graphics3D[{White, Specularity[White], Sphere[{0, 0, 0}, 0.95]}],
    GraphPlot3D[graphfad, VertexCoordinateRules -> vertexrulz]]

enter image description here

We can also check where the periodic boundary is located ( mydateline in the code, which leaks out the module).

Show[Graphics3D[{Green, Sphere[newpoints[[ mydateline]], 0.0005]}], 
    Graphics3D[GraphicsComplex[newpoints, Polygon@newpolys]]]

enter image description here

One you have this spherical or cylindrical mapping, I suggest you then find a way to optimize/smooth out the parametrization, for example as it is done in this question, just dont re-mesh the points in 2D since your mesh keeps the correct 3D topology information

share|improve this answer
    
@BeingHuman fixed typo. Should be working hopefully. If there is any other problems let me know. Or if you find a way to make anything more efficient. The main code definitely needs cleaning up, and it can definitely become more functional. It can't be compiled atm since the adjancy list is ragged, and I haven't found a way to deal with that. So any solutions are very welcome. –  lalmei Jul 3 at 9:04
    
Impressive answer. I wish I could claim I understood all of it. :^) –  Mr.Wizard Jul 23 at 13:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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