Catmull–Clark subdivision surface: Difference between revisions
m
Minor edit.
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
m (Minor edit.) |
||
| (3 intermediate revisions by 3 users not shown) | |||
Line 1:
[[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.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.stream.Collectors;
import java.util.stream.Stream;
public final class CatmullClarkSurfaceSubdivision {
public static void main(String[] args) {
List<Point> points = 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 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) );
List<List<Integer>> faces = List.of( List.of( 0, 1, 2, 3 ),
List.of( 3, 2, 4, 5 ),
List.of( 5, 4, 6, 7 ),
List.of( 7, 0, 3, 5 ),
List.of( 7, 6, 1, 0 ),
List.of( 6, 1, 2, 4 ) );
ListPair result = new ListPair(points, faces);
final int iterations = 1;
for ( int i = 0; i < iterations; i++ ) {
result = catmullClarkSurfaceSubdivision(result.points, result.faces);
}
for ( Point point : result.points ) {
System.out.println(point);
}
System.out.println();
for ( List<Integer> face : result.faces ) {
System.out.println(face);
}
}
private static ListPair catmullClarkSurfaceSubdivision(List<Point> points, List<List<Integer>> faces) {
List<Point> facePoints = facePoints(points, faces);
List<Edge> edges = edges(points, faces);
List<Point> edgePoints = edgePoints(points, facePoints, edges);
List<Point> averageFacePoints = averageFacePoints(points, faces, facePoints);
List<Point> averageEdgeCentres = averageEdgeCentres(points, edges);
List<Integer> facesCount = facesCount(points, faces);
List<Point> newPoints = newPoints(points, facesCount, averageFacePoints, averageEdgeCentres);
List<Integer> facePointIndexes = new ArrayList<Integer>();
int facePointIndex = newPoints.size();
for ( Point facePoint : facePoints ) {
newPoints.add(facePoint);
facePointIndexes.addLast(facePointIndex);
facePointIndex += 1;
}
Map<Long, Integer> edgePointIndexes = new HashMap<Long, Integer>();
for ( int edgeIndex = 0; edgeIndex < edges.size(); edgeIndex++ ) {
final int pointOne = edges.get(edgeIndex).pointOne;
final int pointTwo = edges.get(edgeIndex).pointTwo;
Point edgePoint = edgePoints.get(edgeIndex);
newPoints.add(edgePoint);
edgePointIndexes.put(szudzikHash( new IndexPair(pointOne, pointTwo) ), facePointIndex);
facePointIndex += 1;
}
List<List<Integer>> newFaces = new ArrayList<List<Integer>>();
for ( List<Integer> face : faces ) { // The face can contain any number of points
if ( face.size() >= 3 ) { // A face with 2 or less points does not contribute to the surface
final int fullFace = facePointIndexes.get(faces.indexOf(face));
IndexPair pair = order( new IndexPair(face.getLast(), face.getFirst()) );
int edgePoint = edgePointIndexes.get(szudzikHash(pair));
for ( int i = 0; i < face.size(); i++ ) {
IndexPair previousPair = order( new IndexPair(face.get(i), face.get(( i + 1 ) % face.size())) );
final int previousEdgePoint = edgePointIndexes.get(szudzikHash(previousPair));
newFaces.addLast(List.of( face.get(i), previousEdgePoint, fullFace, edgePoint ));
edgePoint = previousEdgePoint;
}
}
}
return new ListPair(newPoints, newFaces);
}
// Return the new points created by the current iteration of the Catmull-Clark surface subdivision algorithm.
private static List<Point> newPoints(List<Point> points, List<Integer> facesCount,
List<Point> averageFacePoints, List<Point> averageEdgeCentres) {
List<Point> newPoints = new ArrayList<Point>();
for ( int pointIndex = 0; pointIndex < points.size(); pointIndex++ ) {
final int faceCount = facesCount.get(pointIndex);
final double multipleOne = (double) ( faceCount - 3 ) / faceCount;
final double multipleTwo = 1.0 / faceCount;
final double multipleThree = 2.0 / faceCount;
Point currentPoint = points.get(pointIndex);
Point pointOne = currentPoint.scalarMultiply(multipleOne);
Point averageFacePoint = averageFacePoints.get(pointIndex);
Point pointTwo = averageFacePoint.scalarMultiply(multipleTwo);
Point averageEdgeCentre = averageEdgeCentres.get(pointIndex);
Point pointThree = averageEdgeCentre.scalarMultiply(multipleThree);
Point pointFour = pointOne.add(pointTwo);
newPoints.addLast(pointFour.add(pointThree));
}
return newPoints;
}
// Return a list containing a value for each point, which is the number of faces containing that point.
private static List<Integer> facesCount(List<Point> points, List<List<Integer>> faces) {
List<Integer> pointsFaces = Stream.generate( () -> 0 ).limit(points.size()).collect(Collectors.toList());
for ( List<Integer> face : faces ) {
for ( int pointIndex : face ) {
pointsFaces.set(pointIndex, pointsFaces.get(pointIndex) + 1);
}
}
return pointsFaces;
}
// Return a list containing a value for each point,
// which is the average of the centres of the edges containing that point.
private static List<Point> averageEdgeCentres(List<Point> points, List<Edge> edges) {
List<PointWithEdgeCount> pointExtras = Stream.generate( () -> new PointWithEdgeCount(Point.ZERO, 0) )
.limit(points.size()).collect(Collectors.toList());
for ( Edge edge : edges ) {
Point centrePoint = edge.centrePoint;
for ( int pointIndex : List.of( edge.pointOne, edge.pointTwo ) ) {
Point point = pointExtras.get(pointIndex).point;
pointExtras.get(pointIndex).point = point.add(centrePoint);
pointExtras.get(pointIndex).edgeCount += 1;
}
}
List<Point> averageEdgeCentres = new ArrayList<Point>();
for ( PointWithEdgeCount pointExtra : pointExtras ) {
averageEdgeCentres.addLast(pointExtra.point.scalarDivide(pointExtra.edgeCount));
};
return averageEdgeCentres;
}
// Return a list containing a value for each point,
// which is the average of the face points for the faces containing that point.
private static List<Point> averageFacePoints(List<Point> points, List<List<Integer>> faces,
List<Point> facePoints) {
List<PointWithEdgeCount> pointExtras = Stream.generate( () -> new PointWithEdgeCount(Point.ZERO, 0) )
.limit(points.size()).collect(Collectors.toList());
for ( List<Integer> face : faces ) {
Point facePoint = facePoints.get(faces.indexOf(face));
for ( int pointIndex : face ) {
Point pointExtra = pointExtras.get(pointIndex).point;
pointExtras.get(pointIndex).point = pointExtra.add(facePoint);
pointExtras.get(pointIndex).edgeCount += 1;
}
}
List<Point> averageFacePoints = new ArrayList<Point>();
for ( PointWithEdgeCount pointExtra : pointExtras ) {
averageFacePoints.addLast(pointExtra.point.scalarDivide(pointExtra.edgeCount));
}
return averageFacePoints;
}
// Return a list of edge points, where an edge point is the average of the centre of an edge
// and the centre of the line segment joining the face points of its two adjacent faces.
private static List<Point> edgePoints(List<Point> points, List<Point> facePoints, List<Edge> edges) {
List<Point> edgePoints = new ArrayList<Point>();
for ( Edge edge : edges ) {
Point edgeCentre = edge.centrePoint;
Point facePoint1 = facePoints.get(edge.faceOne);
Point facePoint2 = ( edge.faceTwo == FACE_NOT_SET ) ? facePoint1 : facePoints.get(edge.faceTwo);
Point centreFacePoint = facePoint1.centrePoint(facePoint2);
edgePoints.addLast(edgeCentre.centrePoint(centreFacePoint));
}
return edgePoints;
}
// Return a list of edges.
private static List<Edge> edges(List<Point> points, List<List<Integer>> faces) {
List<Edge> partialEdges = new ArrayList<Edge>();
for ( List<Integer> face : faces ) {
for ( int pointIndex = 0; pointIndex < face.size(); pointIndex++ ) {
final int pointIndexOne = face.get(pointIndex);
final int pointIndexTwo = face.get((pointIndex + 1) % face.size());
IndexPair pair = order( new IndexPair(pointIndexOne, pointIndexTwo) );
partialEdges.addLast( new Edge(pair.first, pair.second, faces.indexOf(face), 0, Point.ZERO) );
}
}
Collections.sort(partialEdges);
List<Edge> mergedEdges = new ArrayList<Edge>();
int edgeIndex = 0;
while ( edgeIndex < partialEdges.size() ) {
Edge edgeOne = partialEdges.get(edgeIndex);
if ( edgeIndex < partialEdges.size() - 1 ) {
Edge edgeTwo = partialEdges.get(edgeIndex + 1);
if ( edgeOne.pointOne == edgeTwo.pointOne && edgeOne.pointTwo == edgeTwo.pointTwo ) {
mergedEdges.addLast( new Edge(
edgeOne.pointOne, edgeOne.pointTwo, edgeOne.faceOne, edgeTwo.faceOne, Point.ZERO) );
edgeIndex += 2;
} else {
mergedEdges.addLast( new Edge(
edgeOne.pointOne, edgeOne.pointTwo, edgeOne.faceOne, FACE_NOT_SET, Point.ZERO) );
edgeIndex += 1;
}
} else {
mergedEdges.add( new Edge(
edgeOne.pointOne, edgeOne.pointTwo, edgeOne.faceOne, FACE_NOT_SET, Point.ZERO) );
edgeIndex += 1;
}
}
List<Edge> edges = new ArrayList<Edge>();
for ( Edge mergedEdge : mergedEdges ) {
Point pointOne = points.get(mergedEdge.pointOne);
Point pointTwo = points.get(mergedEdge.pointTwo);
Point centrePoint = pointOne.centrePoint(pointTwo);
edges.addLast( new Edge(
mergedEdge.pointOne, mergedEdge.pointTwo, mergedEdge.faceOne, mergedEdge.faceTwo, centrePoint) );
}
return edges;
}
// Return a list of face points, where a face point is the average of all the points of a face.
private static List<Point> facePoints(List<Point> points, List<List<Integer>> faces) {
List<Point> facePoints = new ArrayList<Point>();
for ( List<Integer> face : faces ) {
Point facePoint = Point.ZERO;
for ( int pointIndex : face ) {
Point point = points.get(pointIndex);
facePoint = facePoint.add(point);
}
facePoints.addLast(facePoint.scalarDivide(face.size()));
}
return facePoints;
}
// Return an IndexPair with its values sorted into ascending order.
private static IndexPair order(IndexPair pair) {
if ( pair.first < pair.second ) {
return pair;
}
return new IndexPair(pair.second, pair.first);
}
// Return a hashed value of the two numbers in the given IndexPair.
private static long szudzikHash(IndexPair pair) {
return ( pair.first >= pair.second ) ?
( pair.first * pair.first ) + pair.first + pair.second : ( pair.second * pair.second ) + pair.first;
}
// A three dimensional point.
private static class Point {
public Point(double aX, double aY, double aZ) {
x = aX; y = aY; z = aZ;
}
public Point add(Point other) {
return new Point(x + other.x, y + other.y, z + other.z);
}
public Point scalarMultiply(double factor) {
return new Point(x * factor, y * factor, z * factor);
}
public Point scalarDivide(double factor) {
return scalarMultiply(1.0 / factor);
}
public Point centrePoint(Point other) {
return add(other).scalarDivide(2.0);
}
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(" %.5f", value) : String.format("%.5f", value);
}
private double x, y, z;
}
// An edge consists of two end points, together with it's one or two adjacent faces
// and the centre point of that edge.
private static class Edge implements Comparable<Edge> {
public Edge(int aPointOne, int aPointTwo, int aFaceOne, int aFaceTwo, Point aCentrePoint) {
pointOne = aPointOne; pointTwo = aPointTwo; faceOne = aFaceOne; faceTwo = aFaceTwo;
centrePoint = aCentrePoint;
}
@Override
public int compareTo(Edge other) {
if ( pointOne == other.pointOne ) {
if ( pointTwo == other.pointTwo ) {
return Integer.compare(faceOne, other.faceOne);
}
return Integer.compare(pointTwo, other.pointTwo);
}
return Integer.compare(pointOne, other.pointOne);
}
private int pointOne, pointTwo, faceOne, faceTwo;
private Point centrePoint;
}
// A point together with the number of edges containing that point.
private static class PointWithEdgeCount {
public PointWithEdgeCount(Point aPoint, int aEdgeCount) {
point = aPoint; edgeCount = aEdgeCount;
}
private Point point;
private int edgeCount;
}
private static record IndexPair(int first, int second) {}
private static record ListPair(List<Point> points, List<List<Integer>> faces) {}
private static final int FACE_NOT_SET = -1;
}
</syntaxhighlight>
{{ out }}
<pre>
(-0.55556, 0.55556, 0.55556)
(-0.55556, -0.55556, 0.55556)
( 0.55556, -0.55556, 0.55556)
( 0.55556, 0.55556, 0.55556)
( 0.55556, -0.55556, -0.55556)
( 0.55556, 0.55556, -0.55556)
(-0.55556, -0.55556, -0.55556)
(-0.55556, 0.55556, -0.55556)
( 0.00000, 0.00000, 1.00000)
( 1.00000, 0.00000, 0.00000)
( 0.00000, 0.00000, -1.00000)
( 0.00000, 1.00000, 0.00000)
(-1.00000, 0.00000, 0.00000)
( 0.00000, -1.00000, 0.00000)
(-0.75000, 0.00000, 0.75000)
( 0.00000, 0.75000, 0.75000)
(-0.75000, 0.75000, 0.00000)
( 0.00000, -0.75000, 0.75000)
(-0.75000, -0.75000, 0.00000)
( 0.75000, 0.00000, 0.75000)
( 0.75000, -0.75000, 0.00000)
( 0.75000, 0.75000, 0.00000)
( 0.75000, 0.00000, -0.75000)
( 0.00000, -0.75000, -0.75000)
( 0.00000, 0.75000, -0.75000)
(-0.75000, 0.00000, -0.75000)
[0, 14, 8, 15]
[1, 17, 8, 14]
[2, 19, 8, 17]
[3, 15, 8, 19]
[3, 19, 9, 21]
[2, 20, 9, 19]
[4, 22, 9, 20]
[5, 21, 9, 22]
[5, 22, 10, 24]
[4, 23, 10, 22]
[6, 25, 10, 23]
[7, 24, 10, 25]
[7, 16, 11, 24]
[0, 15, 11, 16]
[3, 21, 11, 15]
[5, 24, 11, 21]
[7, 25, 12, 16]
[6, 18, 12, 25]
[1, 14, 12, 18]
[0, 16, 12, 14]
[6, 18, 13, 23]
[1, 17, 13, 18]
[2, 20, 13, 17]
[4, 23, 13, 20]
</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,464:
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=
subSetQ[large_,small_List] := And@@(MemberQ[large,#]&/@small)
Line 1,083 ⟶ 1,483:
<syntaxhighlight lang=
edges = DeleteDuplicates[Flatten[Partition[#,2,1,-1]&/@faces,1],Sort[#1]==Sort[#2]&];
avgFacePoints=Mean[Points[[#]]] &/@ faces;
Line 1,113 ⟶ 1,513:
The implimentation can be tested with polygons with and without holes by using the polydata
<syntaxhighlight lang=
Function[iteration,
Line 1,121 ⟶ 1,521:
For a surface with holes the resulting iterative subdivision will be:
<syntaxhighlight lang=
Function[iteration, Graphics3D[
(Polygon[iteration[[1]][[#]]] & /@ iteration[[2]])
Line 1,128 ⟶ 1,528:
This code was written in Mathematica 8.
=={{header|Nim}}==
{{trans|Go}}
{{trans|Python}}
<syntaxhighlight lang=
import tables
Line 1,542 ⟶ 1,941:
[ 2, 20, 13, 17]
[ 4, 23, 13, 20]</pre>
=={{header|OCaml}}==
{{incorrect|OCaml|wrong output data}}
Line 1,552 ⟶ 1,950:
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 ⟶ 2,089:
Another implementation which should work with holes, but has only been tested on a cube
{{works with|OCaml|4.02+}}
<syntaxhighlight lang=
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,198:
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=
<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,471:
{ 4,23,13,20}}
</pre>
=={{header|Python}}==
<syntaxhighlight lang="python">
"""
Line 2,634 ⟶ 3,030:
graph_output(output_points, output_faces)
</syntaxhighlight>
=={{header|Rust}}==
{{trans|Python}}
<syntaxhighlight lang=
pub struct Vector3 {pub x: f64, pub y: f64, pub z: f64, pub w: f64}
Line 2,922 ⟶ 3,317:
}
</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,445:
<center>[[File:Tcl-Catmull.png]]</center>
=={{header|Wren}}==
{{trans|Go}}
Line 3,058 ⟶ 3,451:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang=
import "./sort" for Sort
import "./math" for Int
import "./fmt" for Fmt
var Point = Tuple.create("Point", ["x", "y", "z"])
Line 3,383 ⟶ 3,776:
[ 4 23 13 20]
</pre>
| |||