Catmull–Clark subdivision surface: Difference between revisions

Content deleted Content added
Rdm (talk | contribs)
→‎{{header|Mathematica}}: Replaced code with more new more readable version.
Line 203:
 
=={{header|Mathematica}}==
This implementation supports tris, quads, and higher polys, as well as surfaces with holes.
The function relies on three externally defined general functionality functions:
 
<lang Mathematica>subSetQ[large_,small_] := MemberQ[large,small]
{{improve|Mathematica|Reformat the code to be more readable; to make its structure more visible.}}
subSetQ[large_,small_List] := And@@(MemberQ[large,#]&/@small)
 
containing[groupList_,item_]:= Flatten[Position[groupList,group_/;subSetQ[group,item]]]
This implementation supports tris, quads, and higher polys, as well as surfaces with holes.
 
ReplaceFace[face_]:=Transpose[Prepend[Transpose[{#[[1]],face,#[[2]]}&/@Transpose[Partition[face,2,1,1]//{#,RotateRight[#]}&]],face]]</lang>
subSetQ[small,large] is a boolean test for whether small is a subset of large. Note this is not a general purpose implimentation and only serves this purpose under the constrictions of the following program.
 
containing[{obj1,obj2,...},item] Will return a list of indices of the objects containing item, where objects are faces or edges and item is edges or vertexes.
faces containing a given vertex, faces containing a given edge, edges containing a given point.
It is used for each such purpose in the code called via infix notation, the specific usage is easily distinguised by variable names. For example faces~containing~edge would be a list of the indices for the faces containing the given edge.
 
ReplaceFace[face] replaces the face with a list of descriptions for the new faces. It will return a list containing mixed objects, vertexes, edges and faces where edges and faces referes to the new vertexes to be generated by the code. When the new vertexes have been appended to the updated old vertexes, these mixed objects will be recalcluated into correct indices into the new vertex list by the later defined function newIndex[].
 
 
<lang Mathematica>CatMullClark[{Points_,faces_}]:=Block[{avgFacePoints,avgEdgePoints,updatedPoints,newEdgePoints,newPoints,edges,newFaces,weights,pointUpdate,edgeUpdate,newIndex},
edges = DeleteDuplicates[Flatten[Partition[#,2,1,-1]&/@faces,1],Sort[#1]==Sort[#2]&];
fp avgFacePoints= Mean[vPoints[[#]]] & /@ ifaces;
avgEdgePoints=Mean[Points[[#]]] &/@ edges;
 
weights[vertex_]:= Count[faces,vertex,2]//{(#-3),1,2}/#&;
pointUpdate[vertex_]:=
If[Length[faces~containing~vertex]!=Length[edges~containing~vertex],
Mean[avgEdgePoints[[Select[edges~containing~vertex,holeQ[edges[[#]],faces]&]]]],
Total[weights[vertex]{ Points[[vertex]], Mean[avgFacePoints[[faces~containing~vertex]]], Mean[avgEdgePoints[[edges~containing~vertex]]]}]
];
 
edgeUpdate[edge_]:=
If[Length[faces~containing~edge]==1,
Mean[Points[[edge]]],
Mean[Points[[Flatten[{edge, faces[[faces~containing~edge]]}]]]]
];
 
updatedPoints = pointUpdate/@Range[1,Length[Points]];
newEdgePoints = edgeUpdate/@edges;
newPoints = Join[updatedPoints,avgFacePoints,newEdgePoints];
 
newIndex[edge_/;Length[edge]==2] := Length[Points]+Length[faces]+Position[Sort/@edges,Sort@edge][[1,1]]
newIndex[face_] := Length[Points]+Position[faces,face][[1,1]]
 
newFaces = Flatten[Map[newIndex[#,{Points,edges,faces}]&,ReplaceFace/@faces,{-2}],1];
{newPoints,newFaces}
]</lang>
 
The implimentation can be tested with polygons with and without holes by using the polydata
v<lang Mathematica>{points,faces}= PolyhedronData["Cube", {"VertexCoordinates","FaceIndices"}] // N;
 
Function[iteration,
<lang Mathematica>CatmullClark[{v_, i_}] := Block[{e, vc, fp, ep, vp},
Graphics3D[(Polygon[iteration[[1]][[#]]]&/@iteration[[2]])]
e = Function[a, {a, Select[Transpose[{i, Range@Length@i}],
]/@NestList[CatMullClark,{points,faces},3]//GraphicsRow</lang>
Length@Intersection[#[[1]], a] == 2 &][[All, 2]]}] /@
Union[Sort /@ Flatten[Partition[#, 2, 1, 1] & /@ i, 1]];
vc = Table[{n, Select[Transpose[{i, Range@Length@i}], MemberQ[#[[1]], n] &][[All, 2]],
Select[Transpose[{e[[All, 1]], Range@Length@e[[All, 1]]}], MemberQ[#[[1]], n] &][[All, 2]]}, {n, Length@v}];
fp = Mean[v[[#]]] & /@ i;
ep = If[Length[#[[2]]] == 1, Mean[v[[#[[1]]]]], Mean@Join[v[[#[[1]]]], fp[[#[[2]]]]]] & /@ e;
vp = If[Length[#[[2]]] != Length[#[[3]]], Mean@Join[{v[[#[[1]]]]}, ep[[Select[#[[3]],
Length[e[[#, 2]]] != 2 &]]]], ((Length@#[[2]] - 3) v[[#[[1]]]] + Mean@fp[[#[[2]]]] +
2 Mean@ep[[#[[3]]]])/Length@#[[2]]] & /@ vc;
{Join[vp, ep, fp], Flatten[Function[a, Function[
b, {a[[1]], #[[1]] + Length[vc], b + Length[vc] + Length[e], #[[2]] + Length[vc]} &@
Sort[Select[Transpose[{e, Range@Length@e}], MemberQ[#[[1, 1]], a[[1]]] && MemberQ[#[[1, 2]], b] &],
With[{f = i[[Intersection[#[[1, 2]], #2[[1, 2]]][[1]]]],
n = Intersection[#[[1, 1]], #2[[1, 1]]][[1]]},
Xor[Abs[#] == 1, # < 0] &@(Position[f, Complement[#[[1, 1]], {n}][[1]]] -
Position[f, n])[[1, 1]]] &][[All, 2]]] /@ a[[2]]] /@ vc, 1]}]
v = PolyhedronData["Cube", "VertexCoordinates"] // N
i = PolyhedronData["Cube", "FaceIndices"]
NestList[CatmullClark, {v, i}, 4];
Graphics3D[{FaceForm[{Opacity[0.3]}, {Opacity[0.1]}], GraphicsComplex[#[[1]], Polygon[#[[2]]]]}] & /@ %
Graphics3D[{EdgeForm[], FaceForm[White, Black],
GraphicsComplex[#[[1]], Polygon[#[[2]]], VertexNormals -> #[[1]]]}, Boxed -> False] & /@ %%
</lang>
 
This code was written in Mathematica 8.
The last few lines, after the function definition, do a test by using the built-in polyhedron data to generate the vertices and face indices. Then it repeatedly applies the method and graphs the results. Note that this was written in Mathematica 7, although it should be easy enough to port to maybe v5.2.
 
=={{header|OCaml}}==