src/surface/surface-utils.js

/**
 * @file Surface Utils
 * @author Alexander Rose <alexander.rose@weirdbyte.de>
 * @private
 */

import { degToRad } from '../math/math-utils.js'
import {
    m4new, m4multiply, m4makeTranslation, m4makeScale, m4makeRotationY
} from '../math/matrix-utils.js'
import {
    v3addScalar, v3subScalar, v3divideScalar, v3multiplyScalar,
    v3floor, v3ceil, v3sub, v3negate,
    v3cross, v3fromArray, normalizeVector3array
} from '../math/vector-utils.js'

function laplacianSmooth (verts, faces, numiter, inflate) {
    // based on D. Xu, Y. Zhang (2009) Generating Triangulated Macromolecular
    // Surfaces by Euclidean Distance Transform. PLoS ONE 4(12): e8140.
    //
    // Permission to use, copy, modify, and distribute this program for
    // any purpose, with or without fee, is hereby granted, provided that
    // the notices on the head, the reference information, and this
    // copyright notice appear in all copies or substantial portions of
    // the Software. It is provided "as is" without express or implied
    // warranty.
    //
    // ported to JavaScript and adapted to NGL by Alexander Rose

  numiter = numiter || 1
  inflate = inflate || true

  var nv = verts.length / 3
  var nf = faces.length / 3
  var norms

  if (inflate) {
    norms = new Float32Array(nv * 3)
  }

  var tps = new Float32Array(nv * 3)

  var i
  var ndeg = 20
  var vertdeg = new Array(ndeg)

  for (i = 0; i < ndeg; ++i) {
    vertdeg[ i ] = new Uint32Array(nv)
  }

  for (i = 0; i < nv; ++i) {
    vertdeg[ 0 ][ i ] = 0
  }

  var j, jl
  var flagvert

    // for each face

  for (i = 0; i < nf; ++i) {
    var ao = i * 3
    var bo = i * 3 + 1
    var co = i * 3 + 2

        // vertex a

    flagvert = true
    for (j = 0, jl = vertdeg[ 0 ][ faces[ao] ]; j < jl; ++j) {
      if (faces[ bo ] === vertdeg[ j + 1 ][ faces[ ao ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ ao ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ ao ] ] ][ faces[ ao ] ] = faces[ bo ]
    }

    flagvert = true
    for (j = 0, jl = vertdeg[ 0 ][ faces[ ao ] ]; j < jl; ++j) {
      if (faces[ co ] === vertdeg[ j + 1 ][ faces[ ao ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ ao ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ ao ] ] ][ faces[ ao ] ] = faces[ co ]
    }

        // vertex b

    flagvert = true
    for (j = 0, jl = vertdeg[ 0 ][ faces[ bo ] ]; j < jl; ++j) {
      if (faces[ ao ] === vertdeg[ j + 1 ][ faces[ bo ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ bo ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ bo ] ] ][ faces[ bo ] ] = faces[ ao ]
    }

    flagvert = true
    for (j = 0, jl = vertdeg[ 0 ][ faces[ bo ] ]; j < jl; ++j) {
      if (faces[ co ] === vertdeg[ j + 1 ][ faces[ bo ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ bo ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ bo ] ] ][ faces[ bo ] ] = faces[ co ]
    }

        // vertex c

    flagvert = true
    for (j = 0; j < vertdeg[ 0 ][ faces[ co ] ]; ++j) {
      if (faces[ ao ] === vertdeg[ j + 1 ][ faces[ co ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ co ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ co ] ] ][ faces[ co ] ] = faces[ ao ]
    }

    flagvert = true
    for (j = 0, jl = vertdeg[ 0 ][ faces[ co ] ]; j < jl; ++j) {
      if (faces[ bo ] === vertdeg[ j + 1 ][ faces[ co ] ]) {
        flagvert = false
        break
      }
    }
    if (flagvert) {
      vertdeg[ 0 ][ faces[ co ] ]++
      vertdeg[ vertdeg[ 0 ][ faces[ co ] ] ][ faces[ co ] ] = faces[ bo ]
    }
  }

  var wt = 1.0
  var wt2 = 0.5
  var i3, vi3, vdi, wtvi, wt2vi
  var ssign = -1
  var scaleFactor = 1
  var outwt = 0.75 / (scaleFactor + 3.5)  // area-preserving

    // smoothing iterations

  for (var k = 0; k < numiter; ++k) {
        // for each vertex

    for (i = 0; i < nv; ++i) {
      i3 = i * 3
      vdi = vertdeg[ 0 ][ i ]

      if (vdi < 3) {
        tps[ i3 ] = verts[ i3 ]
        tps[ i3 + 1 ] = verts[ i3 + 1 ]
        tps[ i3 + 2 ] = verts[ i3 + 2 ]
      } else if (vdi === 3 || vdi === 4) {
        tps[ i3 ] = 0
        tps[ i3 + 1 ] = 0
        tps[ i3 + 2 ] = 0

        for (j = 0; j < vdi; ++j) {
          vi3 = vertdeg[ j + 1 ][ i ] * 3
          tps[ i3 ] += verts[ vi3 ]
          tps[ i3 + 1 ] += verts[ vi3 + 1 ]
          tps[ i3 + 2 ] += verts[ vi3 + 2 ]
        }

        tps[ i3 ] += wt2 * verts[ i3 ]
        tps[ i3 + 1 ] += wt2 * verts[ i3 + 1 ]
        tps[ i3 + 2 ] += wt2 * verts[ i3 + 2 ]

        wt2vi = wt2 + vdi
        tps[ i3 ] /= wt2vi
        tps[ i3 + 1 ] /= wt2vi
        tps[ i3 + 2 ] /= wt2vi
      } else {
        tps[ i3 ] = 0
        tps[ i3 + 1 ] = 0
        tps[ i3 + 2 ] = 0

        for (j = 0; j < vdi; ++j) {
          vi3 = vertdeg[ j + 1 ][ i ] * 3
          tps[ i3 ] += verts[ vi3 ]
          tps[ i3 + 1 ] += verts[ vi3 + 1 ]
          tps[ i3 + 2 ] += verts[ vi3 + 2 ]
        }

        tps[ i3 ] += wt * verts[ i3 ]
        tps[ i3 + 1 ] += wt * verts[ i3 + 1 ]
        tps[ i3 + 2 ] += wt * verts[ i3 + 2 ]

        wtvi = wt + vdi
        tps[ i3 ] /= wtvi
        tps[ i3 + 1 ] /= wtvi
        tps[ i3 + 2 ] /= wtvi
      }
    }

    verts.set(tps)  // copy smoothed positions

    if (inflate) {
      computeVertexNormals(verts, faces, norms)
      var nv3 = nv * 3

      for (i3 = 0; i3 < nv3; i3 += 3) {
                // if(verts[i].inout) ssign=1;
                // else ssign=-1;

        verts[ i3 ] += ssign * outwt * norms[ i3 ]
        verts[ i3 + 1 ] += ssign * outwt * norms[ i3 + 1 ]
        verts[ i3 + 2 ] += ssign * outwt * norms[ i3 + 2 ]
      }
    }
  }
}
laplacianSmooth.__deps = [ computeVertexNormals ]

function computeVertexNormals (position, index, normal) {
  var i, il

  if (normal === undefined) {
    normal = new Float32Array(position.length)
  } else {
        // reset existing normals to zero
    for (i = 0, il = normal.length; i < il; i++) {
      normal[ i ] = 0
    }
  }

  var a = new Float32Array(3)
  var b = new Float32Array(3)
  var c = new Float32Array(3)
  var cb = new Float32Array(3)
  var ab = new Float32Array(3)

  if (index) {
        // indexed elements
    for (i = 0, il = index.length; i < il; i += 3) {
      var ai = index[ i ] * 3
      var bi = index[ i + 1 ] * 3
      var ci = index[ i + 2 ] * 3

      v3fromArray(a, position, ai)
      v3fromArray(b, position, bi)
      v3fromArray(c, position, ci)

      v3sub(cb, c, b)
      v3sub(ab, a, b)
      v3cross(cb, cb, ab)

      normal[ ai ] += cb[ 0 ]
      normal[ ai + 1 ] += cb[ 1 ]
      normal[ ai + 2 ] += cb[ 2 ]

      normal[ bi ] += cb[ 0 ]
      normal[ bi + 1 ] += cb[ 1 ]
      normal[ bi + 2 ] += cb[ 2 ]

      normal[ ci ] += cb[ 0 ]
      normal[ ci + 1 ] += cb[ 1 ]
      normal[ ci + 2 ] += cb[ 2 ]
    }
  } else {
        // non-indexed elements (unconnected triangle soup)
    for (i = 0, il = position.length; i < il; i += 9) {
      v3fromArray(a, position, i)
      v3fromArray(b, position, i + 3)
      v3fromArray(c, position, i + 6)

      v3sub(cb, c, b)
      v3sub(ab, a, b)
      v3cross(cb, cb, ab)

      normal[ i ] = cb[ 0 ]
      normal[ i + 1 ] = cb[ 1 ]
      normal[ i + 2 ] = cb[ 2 ]

      normal[ i + 3 ] = cb[ 0 ]
      normal[ i + 4 ] = cb[ 1 ]
      normal[ i + 5 ] = cb[ 2 ]

      normal[ i + 6 ] = cb[ 0 ]
      normal[ i + 7 ] = cb[ 1 ]
      normal[ i + 8 ] = cb[ 2 ]
    }
  }

  normalizeVector3array(normal)

  return normal
}
computeVertexNormals.__deps = [
  v3sub, v3cross, v3fromArray, normalizeVector3array
]

function getRadiusDict (radiusList) {
  var radiusDict = {}
  for (var i = 0, il = radiusList.length; i < il; ++i) {
    radiusDict[ radiusList[ i ] ] = true
  }
  return radiusDict
}

function getSurfaceGrid (min, max, maxRadius, scaleFactor, extraMargin) {
    // need margin to avoid boundary/round off effects
  var margin = (1 / scaleFactor) * 3
  margin += maxRadius

  v3subScalar(min, min, extraMargin + margin)
  v3addScalar(max, max, extraMargin + margin)

  v3multiplyScalar(min, min, scaleFactor)
  v3floor(min, min)
  v3divideScalar(min, min, scaleFactor)

  v3multiplyScalar(max, max, scaleFactor)
  v3ceil(max, max)
  v3divideScalar(max, max, scaleFactor)

  var dim = new Float32Array(3)
  v3sub(dim, max, min)
  v3multiplyScalar(dim, dim, scaleFactor)
  v3ceil(dim, dim)
  v3addScalar(dim, dim, 1)

  var maxSize = Math.pow(10, 6) * 256
  var tmpSize = dim[ 0 ] * dim[ 1 ] * dim[ 2 ] * 3

  if (maxSize <= tmpSize) {
    scaleFactor *= Math.pow(maxSize / tmpSize, 1 / 3)

    v3multiplyScalar(min, min, scaleFactor)
    v3floor(min, min)
    v3divideScalar(min, min, scaleFactor)

    v3multiplyScalar(max, max, scaleFactor)
    v3ceil(max, max)
    v3divideScalar(max, max, scaleFactor)

    v3sub(dim, max, min)
    v3multiplyScalar(dim, dim, scaleFactor)
    v3ceil(dim, dim)
    v3addScalar(dim, dim, 1)
  }

  var tran = new Float32Array(min)
  v3negate(tran, tran)

    // coordinate transformation matrix
  var matrix = m4new()
  var mroty = m4new()
  m4makeRotationY(mroty, degToRad(90))
  m4multiply(matrix, matrix, mroty)

  var mscale = m4new()
  m4makeScale(
    mscale,
    -1 / scaleFactor,
    1 / scaleFactor,
    1 / scaleFactor
  )
  m4multiply(matrix, matrix, mscale)

  var mtrans = m4new()
  m4makeTranslation(
    mtrans,
    -scaleFactor * tran[2],
    -scaleFactor * tran[1],
    -scaleFactor * tran[0]
  )
  m4multiply(matrix, matrix, mtrans)

  return {
    dim: dim,
    tran: tran,
    matrix: matrix,
    scaleFactor: scaleFactor
  }
}
getSurfaceGrid.__deps = [
  degToRad,
  v3subScalar, v3addScalar, v3divideScalar, v3multiplyScalar,
  v3floor, v3ceil, v3sub, v3negate,
  m4new, m4multiply, m4makeTranslation, m4makeScale, m4makeRotationY
]

export {
    laplacianSmooth,
    computeVertexNormals,
    getRadiusDict,
    getSurfaceGrid
}