Created
January 17, 2014 00:37
-
-
Save Aatch/8466307 to your computer and use it in GitHub Desktop.
Surface Simplification with Quadric Error Metrics - for BlackMoon
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
; This is written pseudo-code, so hopefully should be fairly self-evident. | |
; | |
; A few notes about the general maths which I think will be needed: | |
; | |
; Matrices use what is called a "row-major" convention for sizes/indexes, | |
; this just means that a 2x3 matrix has 2 rows and 3 columns, it also | |
; means that for a matrix, M, M[1][2] is the number in the second column | |
; in the first row (I will also use base-1 indexing, since it better | |
; matches mathematical convention) | |
; | |
; Vectors can be either row vectors or column vectors, this just means | |
; that a n-length vector can represent either a 1xn matrix (row vector) | |
; or a nx1 matrix (column vector). I'll write all the vectors as row | |
; vectors, however the 'transpose' method will convert them between row | |
; and column forms, which is necessary at times. | |
; | |
; Since the cross product is only defined for 3-dimensional vectors, | |
; multiplication will be standard matrix multiplication. The cross | |
; product, 'A x B' will be done as A.cross(B). | |
; Simplifies the given mesh, performing n_contractions contractions | |
; | |
; threshold is used for controlling whether | |
; two vertices that are *not* connected by an edge are considered for | |
; contraction. It is the maximum distance between the two vertices, a value | |
; of 0.0 will make it only consider edge-connected vertices. | |
fun simplify_mesh(mesh, n_contractions, threshold) { | |
errors = []; | |
; Build up the error quadrics | |
for each vertex in mesh.vertices { | |
errors.append(calculate_quadric(vertex)) | |
} | |
to_contract = [] | |
; Go through and check all the pairs of vertices, in reality | |
; there would be a pre-calculated list of valid pairs already | |
for each (v1, v2) in mesh.vertices.pairs { | |
Q = errors.for(v1) + errors.for(v2) | |
; So there is a way to find an optimum position to contract to | |
; by treating it as a system of linear equations to be solved, | |
; I'll outline it here, however it may not be worthwhile (or | |
; possible) to do in realtime. | |
; The optimum position for the target vertex, vt, is the minimisation | |
; of the error, this satisfies the matrix equation: | |
; | |
; ( Q11 Q12 Q13 Q14 ) ( 0 ) | |
; ( Q21 Q22 Q23 Q24 ) * vt = ( 0 ) | |
; ( Q13 Q14 Q33 Q34 ) ( 0 ) | |
; ( 0 0 0 1 ) ( 1 ) | |
; | |
; Which can be rearranged to: | |
; | |
; ( Q11 Q12 Q13 Q14 ) -1 ( 0 ) | |
; vt = ( Q21 Q22 Q23 Q24 ) * ( 0 ) | |
; ( Q13 Q14 Q33 Q34 ) ( 0 ) | |
; ( 0 0 0 1 ) ( 1 ) | |
; | |
; The '-1' up there means "inversion", the reason the bottom row is | |
; empty is because vt is in homogenous coordinates and therefore has | |
; a 'w' component of '1' | |
; | |
; However, not all matrices are invertible, so we have to fallback to | |
; an approximation when that happens. | |
if Q.is_invertible() { | |
Qt = Q.copy() | |
Qt.set_row(4, [0, 0, 0, 1]) | |
vt = Qt.invert() * vector(0, 0, 0, 1).transpose() | |
; vt is a column vector here, so this is row_vector*matrix*column_vector, | |
; which produces a 1x1 matrix, which is just a number. | |
error = vt.transpose() * Q * vt | |
} else { | |
; If we can't find an optimum, then use the minimum error for the | |
; three options of v1, v2 and (v1 - v2)/2 | |
vavg = (v1 - v2)/2 | |
v1_error = v1 * Q * v1.transpose() | |
v2_error = v2 * Q * v2.transpose() | |
vavg_error = vavg * Q * vavg.transpose() | |
if v1_error is smallest { | |
error = v1_error | |
vt = v1 | |
} else if v2_error is smallest { | |
error = v2_error | |
vt = v2 | |
} else { | |
error = vavg_error | |
vt = vavg | |
} | |
} | |
; Add this pair to the list of contractions, dropping other pairs | |
; out if they have a higher error | |
if to_contract.length < n_contractions { | |
to_contract.append([vt, error, v1, v2]) | |
} else if error is smaller than largest error in to_contract { | |
to_contract.remove_largest_error(); | |
to_contract.append([vt, error, v1, v2]) | |
} | |
} | |
} | |
; Now go and contract each pair to the target position | |
for each (target, error, v1, v2) in to_contract { | |
contract(v1, v2, target) | |
} | |
} | |
; Calculates the quadric error metric for the given vertex | |
funcalculate_quadric(vertex) { | |
; Q is the final error metric, it is a symmetric matrix, meaning | |
; that M[i][j] = M[j][i] | |
Q = new SymmetrixMatrix(4,4) | |
; Loop through each of the faces that share this vertex | |
; and get the two other vertices for that face. | |
for each (v1, v2) in vertex.faces.other_vertices { | |
; Kp is the "fundamental error quadric", this is just a matrix | |
; that can be used to find the square of the distance between and | |
; a point. | |
Kp = new SymmetricMatrix(4,4) | |
plane = plane_vector(vertex, v1, v2) | |
; So the next bit is what the following operation is *actually* doing, | |
; note that this isn't the dot product as the left hand side is a column | |
; vector, meaning you get a (in this case) 4x4 matrix. | |
; Kp = plane.transpose()*plane | |
Kp.set_row(1, [plane[1]*plane[1], plane[1]*plane[2], plane[1]*plane[3], plane[1]*plane[4]]) | |
Kp.set_row(2, [plane[2]*plane[1], plane[2]*plane[2], plane[2]*plane[3], plane[2]*plane[4]]) | |
Kp.set_row(3, [plane[3]*plane[1], plane[3]*plane[2], plane[3]*plane[3], plane[3]*plane[4]]) | |
Kp.set_row(4, [plane[4]*plane[1], plane[4]*plane[2], plane[4]*plane[3], plane[4]*plane[4]]) | |
; In reality, the above only needs to do 10 of those multiplication, since the other 6 are | |
; actually identical numbers | |
; Turns out that the sum of the fundamental error quadrics produces a quadric | |
; that can be used to find the sum of the squared distances to all the planes | |
; at once. I.E. if Kp1 can be used to find the squared distance between p1 and v, | |
; then Kp1+Kp2 can be used to find the squared distance between p1 and v plus the | |
; squared distance between p2 and v. | |
Q += Kp | |
} | |
return Q | |
} | |
; Returns a vector (a, b, c, d) for the plane defined by the | |
; equation 'ax + by + cz + d = 0', where a^2 + b^2 + c^2 = 1; | |
; | |
; The definition of a plane above is: the point (x, y, z) is on | |
; the plane, p, if the equation `ax + by + cz + d = 0` holds. This | |
; means that the coefficients a, b and c can be interpreted as | |
; rotations, such that the vector (a, b, c) is the normal of the | |
; plane. 'd' here is the offset of the plane from the origin in the | |
; direction of the normal (since otherwise all planes would contain | |
; the origin) | |
fun plane_vector(A, B, C) { | |
edge_ab = (B - A) ; Edge from A to B | |
edge_ac = (C - A) ; Edge from A to C | |
; The cross product of the two edges is normal of the plane | |
; they share, this gets the a, b and c components | |
normal = edge_ab.cross(edge_ac).normalize() | |
; since A is a row vector, we need to transpose the normal into a | |
; column vector to get the equivalent to the dot product. | |
; This gets us the offset from the origin, i.e. 'd' | |
d = -A*normal.transpose() | |
return vector(normal[1], normal[2], normal[3], d) | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment