Surface Simplification with Quadric Error Metrics - for BlackMoon

gistfile1.txt

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