Rosetta Code

Catmull–Clark subdivision surface: Difference between revisions

Page Discussion

Catmull–Clark subdivision surface (view source)

Revision as of 15:58, 23 May 2024

11,953 bytes added ,  28 days ago
m
Rewritten code to deal with holes in the surface.
m (syntax highlighting fixup automation)
m (Rewritten code to deal with holes in the surface.)
(4 intermediate revisions by 3 users not shown)
Line 1:
{{task|3D}}[[Category:Graphics algorithms]]{{requires|Graphics}}
[[Category:Geometry]]
{{task|3D}}{{requires|Graphics}}
Implement the Catmull-Clark surface subdivision ([[wp:Catmull–Clark_subdivision_surface|description on Wikipedia]]), which is an algorithm that maps from a surface (described as a set of points and a set of polygons with vertices at those points) to another more refined surface. The resulting surface will always consist of a mesh of quadrilaterals.
 
Line 42 ⟶ 44:
 
For edges and vertices not next to a hole, the standard algorithm from above is used.
 
=={{header|C}}==
[[file:catmull-C.png|center]]
Only the subdivision part. The [[Catmull–Clark_subdivision_surface/C|full source]] is way too long to be shown here. Lots of macros, you'll have to see the full code to know what's what.
<syntaxhighlight lang="c">vertex face_point(face f)
{
int i;
Line 139 ⟶ 140:
return nm;
}</syntaxhighlight>
 
=={{header|Go}}==
{{trans|Python}}
Line 146:
 
See the original version for full comments.
<syntaxhighlight lang="go">package main
 
import (
Line 487:
[ 4 23 13 20]
</pre>
 
=={{header|Haskell}}==
<syntaxhighlight lang="haskell">{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
 
Line 853 ⟶ 852:
(22,[24,15,5,8])
(23,[22,8,5,7])</pre>
 
=={{header|J}}==
 
<syntaxhighlight lang="j">avg=: +/ % #
 
havePoints=: e."1/~ i.@#
Line 882 ⟶ 880:
Example use:
 
<syntaxhighlight lang="j">NB.cube
points=: _1+2*#:i.8
mesh=: 1 A."1 I.(,1-|.)8&$@#&0 1">4 2 1
Line 915 ⟶ 913:
│ │ 0.555556 0.555556 0.555556│
└──────────┴─────────────────────────────┘</syntaxhighlight>
 
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Objects;
import java.util.function.Function;
import java.util.stream.Collectors;
 
public final class CatmullClarkSurfaceSubdivision {
 
public static void main(String[] args) {
List<Face> faces = List.of( new Face(List.of( new Point(-1.0, 1.0, 1.0),
new Point(-1.0, -1.0, 1.0),
new Point( 1.0, -1.0, 1.0),
new Point( 1.0, 1.0, 1.0) )),
new Face(List.of( new Point( 1.0, 1.0, 1.0),
new Point( 1.0, -1.0, 1.0),
new Point( 1.0, -1.0, -1.0),
new Point( 1.0, 1.0, -1.0) )),
new Face(List.of( new Point( 1.0, 1.0, -1.0),
new Point( 1.0, -1.0, -1.0),
new Point(-1.0, -1.0, -1.0),
new Point(-1.0, 1.0, -1.0) )),
new Face(List.of( new Point(-1.0, 1.0, -1.0),
new Point(-1.0, 1.0, 1.0),
new Point( 1.0, 1.0, 1.0),
new Point( 1.0, 1.0, -1.0) )),
new Face(List.of( new Point(-1.0, 1.0, -1.0),
new Point(-1.0, -1.0, -1.0),
new Point(-1.0, -1.0, 1.0),
new Point(-1.0, 1.0, 1.0) )),
new Face(List.of( new Point(-1.0, -1.0, -1.0),
new Point(-1.0, -1.0, 1.0),
new Point( 1.0, -1.0, 1.0),
new Point( 1.0, -1.0, -1.0) )) );
displaySurface(faces);
final int iterations = 1;
for ( int i = 0; i < iterations; i++ ) {
faces = catmullClarkSurfaceSubdivision(faces);
}
displaySurface(faces);
}
// The Catmull-Clarke surface subdivision algorithm.
private static List<Face> catmullClarkSurfaceSubdivision(List<Face> faces) {
// Determine, for each edge, whether or not it is an edge of a hole, and set its edge point accordingly.
List<Edge> edges = faces.stream().map( face -> face.edges.stream() ).flatMap(Function.identity()).toList();
for ( Edge edge : edges ) {
List<Point> facePointsForEdge =
faces.stream().filter( face -> face.contains(edge) ).map( face -> face.facePoint ).toList();
if ( facePointsForEdge.size() == 2 ) {
edge.holeEdge = false;
edge.edgePoint = centroid(List.of( edge.midEdge, centroid(facePointsForEdge) ));
} else {
edge.holeEdge = true;
edge.edgePoint = edge.midEdge;
}
}
Map<Point, Point> nextVertices = nextVertices(edges, faces);
List<Face> nextFaces = new ArrayList<Face>();
for ( Face face : faces ) { // The face may contain any number of points
if ( face.vertices.size() >= 3 ) { // A face with 2 or fewer points does not contribute to the surface
Point facePoint = face.facePoint;
for ( int i = 0; i < face.vertices.size(); i++ ) {
nextFaces.addLast( new Face(List.of(
nextVertices.get(face.vertices.get(i)),
face.edges.get(i).edgePoint,
facePoint,
face.edges.get(Math.floorMod(i - 1, face.vertices.size())).edgePoint )) );
}
}
}
return nextFaces;
}
// Return a map containing, for each vertex,
// the new vertex created by the current iteration of the Catmull-Clark surface subdivision algorithm.
private static Map<Point, Point> nextVertices(List<Edge> edges, List<Face> faces) {
Map<Point, Point> nextVertices = new HashMap<Point, Point>();
List<Point> vertices =
faces.stream().map( face -> face.vertices.stream() ).flatMap(Function.identity()).distinct().toList();
for ( Point vertex : vertices ) {
List<Face> facesForVertex = faces.stream().filter( face -> face.contains(vertex) ).toList();
List<Edge> edgesForVertex = edges.stream().filter( edge -> edge.contains(vertex) ).distinct().toList();
if ( facesForVertex.size() != edgesForVertex.size() ) {
List<Point> midEdgeOfHoleEdges = edgesForVertex.stream().filter( edge -> edge.holeEdge )
.map( edge -> edge.midEdge ).collect(Collectors.toList());
midEdgeOfHoleEdges.add(vertex);
nextVertices.put(vertex, centroid(midEdgeOfHoleEdges));
} else {
final int faceCount = facesForVertex.size();
final double multipleOne = (double) ( faceCount - 3 ) / faceCount;
final double multipleTwo = 1.0 / faceCount;
final double multipleThree = 2.0 / faceCount;
Point nextVertexOne = vertex.multiply(multipleOne);
List<Point> facePoints = facesForVertex.stream().map( face -> face.facePoint ).toList();
Point nextVertexTwo = centroid(facePoints).multiply(multipleTwo);
List<Point> midEdges = edgesForVertex.stream().map( edge -> edge.midEdge ).toList();
Point nextVertexThree = centroid(midEdges).multiply(multipleThree);
Point nextVertexFour = nextVertexOne.add(nextVertexTwo);
nextVertices.put(vertex, nextVertexFour.add(nextVertexThree));
}
}
return nextVertices;
}
// Return the centroid point of the given list of points.
private static Point centroid(List<Point> points) {
return points.stream().reduce(Point.ZERO, (left, right) -> left.add(right) ).divide(points.size());
}
// Display the current Catmull-Clark surface on the console.
private static void displaySurface(List<Face> faces) {
System.out.println("Surface {");
faces.stream().forEach(System.out::println);
System.out.println("}" + System.lineSeparator());
}
// A point of the current Catmull-Clark surface.
private static final class Point implements Comparable<Point> {
public Point(double aX, double aY, double aZ) {
x = aX; y = aY; z = aZ;
}
@Override
public int compareTo(Point other) {
if ( x == other.x ) {
if ( y == other.y ) {
return Double.compare(z, other.z);
}
return Double.compare(y, other.y);
}
return Double.compare(x, other.x);
}
@Override
public boolean equals(Object other) {
return switch ( other ) {
case Point point -> x == point.x && y == point.y && z == point.z;
default -> false;
};
}
@Override
public int hashCode() {
return Objects.hash(x, y, z);
}
public Point add(Point other) {
return new Point(x + other.x, y + other.y, z + other.z);
}
public Point multiply(double factor) {
return new Point(x * factor, y * factor, z * factor);
}
public Point divide(double factor) {
return multiply(1.0 / factor);
}
public String toString() {
return "(" + format(x) + ", " + format(y) + ", " + format(z) + ")";
}
public static Point ZERO = new Point(0.0, 0.0, 0.0);
private String format(double value) {
return ( value >= 0 ) ? String.format(" %.3f", value) : String.format("%.3f", value);
}
private final double x, y, z;
}
// An edge of the current Catmull-Clark surface.
private static final class Edge {
public Edge(Point aBegin, Point aEnd) {
if ( aBegin.compareTo(aEnd) <= 0 ) {
begin = aBegin; end = aEnd;
} else {
begin = aEnd; end = aBegin;
}
midEdge = centroid(List.of( begin, end ));
}
 
@Override
public boolean equals(Object other) {
return switch ( other ) {
case Edge edge -> begin.equals(edge.begin) && end.equals(edge.end);
default -> false;
};
}
@Override
public int hashCode() {
return Objects.hash(begin, end);
}
public boolean contains(Point point) {
return point.equals(begin) || point.equals(end);
}
public boolean holeEdge;
public Point edgePoint;
public final Point begin, end, midEdge;
}
// A face of the current Catmull-Clark surface.
private static final class Face {
public Face(List<Point> aVertices) {
vertices = new ArrayList<Point>(aVertices);
facePoint = centroid(vertices);
edges = new ArrayList<Edge>();
for ( int i = 0; i < vertices.size() - 1; i++ ) {
edges.addLast( new Edge(vertices.get(i), vertices.get(i + 1)) );
}
edges.addLast( new Edge(vertices.getLast(), vertices.getFirst()) );;
}
public boolean contains(Point vertex) {
return vertices.contains(vertex);
}
public boolean contains(Edge edge) {
return contains(edge.begin) && contains(edge.end);
}
public String toString() {
return "Face: " + vertices.stream().map( point -> point.toString() ).collect(Collectors.joining("; "));
}
public final List<Point> vertices;
public final Point facePoint;
public final List<Edge> edges;
}
 
}
</syntaxhighlight>
{{ out }}
<pre>
Surface {
Face: (-1.000, 1.000, 1.000); (-1.000, -1.000, 1.000); ( 1.000, -1.000, 1.000); ( 1.000, 1.000, 1.000)
Face: ( 1.000, 1.000, 1.000); ( 1.000, -1.000, 1.000); ( 1.000, -1.000, -1.000); ( 1.000, 1.000, -1.000)
Face: ( 1.000, 1.000, -1.000); ( 1.000, -1.000, -1.000); (-1.000, -1.000, -1.000); (-1.000, 1.000, -1.000)
Face: (-1.000, 1.000, -1.000); (-1.000, 1.000, 1.000); ( 1.000, 1.000, 1.000); ( 1.000, 1.000, -1.000)
Face: (-1.000, 1.000, -1.000); (-1.000, -1.000, -1.000); (-1.000, -1.000, 1.000); (-1.000, 1.000, 1.000)
Face: (-1.000, -1.000, -1.000); (-1.000, -1.000, 1.000); ( 1.000, -1.000, 1.000); ( 1.000, -1.000, -1.000)
}
 
Surface {
Face: (-0.556, 0.556, 0.556); (-0.750, 0.000, 0.750); ( 0.000, 0.000, 1.000); ( 0.000, 0.750, 0.750)
Face: (-0.556, -0.556, 0.556); ( 0.000, -0.750, 0.750); ( 0.000, 0.000, 1.000); (-0.750, 0.000, 0.750)
Face: ( 0.556, -0.556, 0.556); ( 0.750, 0.000, 0.750); ( 0.000, 0.000, 1.000); ( 0.000, -0.750, 0.750)
Face: ( 0.556, 0.556, 0.556); ( 0.000, 0.750, 0.750); ( 0.000, 0.000, 1.000); ( 0.750, 0.000, 0.750)
Face: ( 0.556, 0.556, 0.556); ( 0.750, 0.000, 0.750); ( 1.000, 0.000, 0.000); ( 0.750, 0.750, 0.000)
Face: ( 0.556, -0.556, 0.556); ( 0.750, -0.750, 0.000); ( 1.000, 0.000, 0.000); ( 0.750, 0.000, 0.750)
Face: ( 0.556, -0.556, -0.556); ( 0.750, 0.000, -0.750); ( 1.000, 0.000, 0.000); ( 0.750, -0.750, 0.000)
Face: ( 0.556, 0.556, -0.556); ( 0.750, 0.750, 0.000); ( 1.000, 0.000, 0.000); ( 0.750, 0.000, -0.750)
Face: ( 0.556, 0.556, -0.556); ( 0.750, 0.000, -0.750); ( 0.000, 0.000, -1.000); ( 0.000, 0.750, -0.750)
Face: ( 0.556, -0.556, -0.556); ( 0.000, -0.750, -0.750); ( 0.000, 0.000, -1.000); ( 0.750, 0.000, -0.750)
Face: (-0.556, -0.556, -0.556); (-0.750, 0.000, -0.750); ( 0.000, 0.000, -1.000); ( 0.000, -0.750, -0.750)
Face: (-0.556, 0.556, -0.556); ( 0.000, 0.750, -0.750); ( 0.000, 0.000, -1.000); (-0.750, 0.000, -0.750)
Face: (-0.556, 0.556, -0.556); (-0.750, 0.750, 0.000); ( 0.000, 1.000, 0.000); ( 0.000, 0.750, -0.750)
Face: (-0.556, 0.556, 0.556); ( 0.000, 0.750, 0.750); ( 0.000, 1.000, 0.000); (-0.750, 0.750, 0.000)
Face: ( 0.556, 0.556, 0.556); ( 0.750, 0.750, 0.000); ( 0.000, 1.000, 0.000); ( 0.000, 0.750, 0.750)
Face: ( 0.556, 0.556, -0.556); ( 0.000, 0.750, -0.750); ( 0.000, 1.000, 0.000); ( 0.750, 0.750, 0.000)
Face: (-0.556, 0.556, -0.556); (-0.750, 0.000, -0.750); (-1.000, 0.000, 0.000); (-0.750, 0.750, 0.000)
Face: (-0.556, -0.556, -0.556); (-0.750, -0.750, 0.000); (-1.000, 0.000, 0.000); (-0.750, 0.000, -0.750)
Face: (-0.556, -0.556, 0.556); (-0.750, 0.000, 0.750); (-1.000, 0.000, 0.000); (-0.750, -0.750, 0.000)
Face: (-0.556, 0.556, 0.556); (-0.750, 0.750, 0.000); (-1.000, 0.000, 0.000); (-0.750, 0.000, 0.750)
Face: (-0.556, -0.556, -0.556); (-0.750, -0.750, 0.000); ( 0.000, -1.000, 0.000); ( 0.000, -0.750, -0.750)
Face: (-0.556, -0.556, 0.556); ( 0.000, -0.750, 0.750); ( 0.000, -1.000, 0.000); (-0.750, -0.750, 0.000)
Face: ( 0.556, -0.556, 0.556); ( 0.750, -0.750, 0.000); ( 0.000, -1.000, 0.000); ( 0.000, -0.750, 0.750)
Face: ( 0.556, -0.556, -0.556); ( 0.000, -0.750, -0.750); ( 0.000, -1.000, 0.000); ( 0.750, -0.750, 0.000)
}
</pre>
 
=={{header|Julia}}==
<syntaxhighlight lang="julia">using Makie, Statistics
 
# Point3f0 is a 3-tuple of 32-bit floats for 3-dimensional space, and all Points are 3D.
Line 1,063 ⟶ 1,360:
println("Press Enter to continue", readline())
</syntaxhighlight>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
This implementation supports tris, quads, and higher polys, as well as surfaces with holes.
The function relies on three externally defined general functionality functions:
 
<syntaxhighlight lang=Mathematica"mathematica">subSetQ[large_,small_] := MemberQ[large,small]
subSetQ[large_,small_List] := And@@(MemberQ[large,#]&/@small)
 
Line 1,083 ⟶ 1,379:
 
 
<syntaxhighlight lang=Mathematica"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]&];
avgFacePoints=Mean[Points[[#]]] &/@ faces;
Line 1,113 ⟶ 1,409:
 
The implimentation can be tested with polygons with and without holes by using the polydata
<syntaxhighlight lang=Mathematica"mathematica">{points,faces}=PolyhedronData["Cube",{"VertexCoordinates","FaceIndices"}];
 
Function[iteration,
Line 1,121 ⟶ 1,417:
 
For a surface with holes the resulting iterative subdivision will be:
<syntaxhighlight lang=Mathematica"mathematica">faces = Delete[faces, 6];
Function[iteration, Graphics3D[
(Polygon[iteration[[1]][[#]]] & /@ iteration[[2]])
Line 1,128 ⟶ 1,424:
 
This code was written in Mathematica 8.
 
=={{header|Nim}}==
{{trans|Go}}
{{trans|Python}}
<syntaxhighlight lang=Nim"nim">import algorithm
import tables
 
Line 1,542 ⟶ 1,837:
[ 2, 20, 13, 17]
[ 4, 23, 13, 20]</pre>
 
=={{header|OCaml}}==
{{incorrect|OCaml|wrong output data}}
Line 1,552 ⟶ 1,846:
In the [[Catmull–Clark subdivision surface/OCaml|sub-page]] there is also a program in OCaml+OpenGL which displays a cube subdivided 2 times with this algorithm.
 
<syntaxhighlight lang="ocaml">open Dynar
 
let add3 (x1, y1, z1) (x2, y2, z2) (x3, y3, z3) =
Line 1,691 ⟶ 1,985:
Another implementation which should work with holes, but has only been tested on a cube
{{works with|OCaml|4.02+}}
<syntaxhighlight lang=OCaml"ocaml">type point = { x: float; y : float; z : float }
let zero = { x = 0.0; y = 0.0; z = 0.0 }
let add a b = { x = a.x+.b.x; y = a.y+.b.y; z = a.z+.b.z }
Line 1,800 ⟶ 2,094:
Face: (0.7500, 0.0000, 0.7500) (0.5556, -0.5556, 0.5556) (0.0000, -0.7500, 0.7500) (0.0000, 0.0000, 1.0000)
}</pre>
 
=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang=Phix"phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">-- demo\rosetta\Catmull_Clark_subdivision_surface.exw</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
Line 2,074 ⟶ 2,367:
{ 4,23,13,20}}
</pre>
 
=={{header|Python}}==
<syntaxhighlight lang="python">
"""
 
Line 2,634 ⟶ 2,926:
graph_output(output_points, output_faces)
</syntaxhighlight>
 
=={{header|Rust}}==
{{trans|Python}}
<syntaxhighlight lang=Rust"rust">
pub struct Vector3 {pub x: f64, pub y: f64, pub z: f64, pub w: f64}
 
Line 2,922 ⟶ 3,213:
}
</syntaxhighlight>
 
=={{header|Tcl}}==
This code handles both holes and arbitrary polygons in the input data.
<syntaxhighlight lang="tcl">package require Tcl 8.5
 
# Use math functions and operators as commands (Lisp-like).
Line 3,051 ⟶ 3,341:
 
<center>[[File:Tcl-Catmull.png]]</center>
 
=={{header|Wren}}==
{{trans|Go}}
Line 3,058 ⟶ 3,347:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang=ecmascript"wren">import "./dynamic" for Tuple, Struct
import "./sort" for Sort
import "./math" for Int
import "./fmt" for Fmt
 
var Point = Tuple.create("Point", ["x", "y", "z"])
Line 3,383 ⟶ 3,672:
[ 4 23 13 20]
</pre>
 
[[Category:Geometry]]
PSNOW123
908

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/327193...364663"