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