surfacenets.js

"use strict"

module.exports = surfaceNets

var generateContourExtractor = require("ndarray-extract-contour")
var triangulateCube = require("triangulate-hypercube")
var zeroCrossings = require("zero-crossings")

function buildSurfaceNets(order, dtype) {
  var dimension = order.length
  var code = ["'use strict';"]
  var funcName = "surfaceNets" + order.join("_") + "d" + dtype

  //Contour extraction function
  code.push(
    "var contour=genContour({",
      "order:[", order.join(), "],",
      "scalarArguments: 3,",
      "phase:function phaseFunc(p,a,b,c) { return (p > c)|0 },")
  if(dtype === "generic") {
    code.push("getters:[0],")
  }

  //Generate vertex function
  var cubeArgs = []
  var extraArgs = []
  for(var i=0; i<dimension; ++i) {
    cubeArgs.push("d" + i)
    extraArgs.push("d" + i)
  }
  for(var i=0; i<(1<<dimension); ++i) {
    cubeArgs.push("v" + i)
    extraArgs.push("v" + i)
  }
  for(var i=0; i<(1<<dimension); ++i) {
    cubeArgs.push("p" + i)
    extraArgs.push("p" + i)
  }
  cubeArgs.push("a", "b", "c")
  extraArgs.push("a", "c")
  code.push("vertex:function vertexFunc(", cubeArgs.join(), "){")
  //Mask args together
  var maskStr = []
  for(var i=0; i<(1<<dimension); ++i) {
    maskStr.push("(p" + i + "<<" + i + ")")
  }
  //Generate variables and giganto switch statement
  code.push("var m=(", maskStr.join("+"), ")|0;if(m===0||m===", (1<<(1<<dimension))-1, "){return}")
  var extraFuncs = []
  var currentFunc = []
  if(1<<(1<<dimension) <= 128) {
    code.push("switch(m){")
    currentFunc = code
  } else {
    code.push("switch(m>>>7){")
  }
  for(var i=0; i<1<<(1<<dimension); ++i) {
    if(1<<(1<<dimension) > 128) {
      if((i%128)===0) {
        if(extraFuncs.length > 0) {
          currentFunc.push("}}")
        }
        var efName = "vExtra" + extraFuncs.length
        code.push("case ", (i>>>7), ":", efName, "(m&0x7f,", extraArgs.join(), ");break;")
        currentFunc = [
          "function ", efName, "(m,", extraArgs.join(), "){switch(m){"
        ]
        extraFuncs.push(currentFunc)
      }  
    }
    currentFunc.push("case ", (i&0x7f), ":")
    var crossings = new Array(dimension)
    var denoms = new Array(dimension)
    var crossingCount = new Array(dimension)
    var bias = new Array(dimension)
    var totalCrossings = 0
    for(var j=0; j<dimension; ++j) {
      crossings[j] = []
      denoms[j] = []
      crossingCount[j] = 0
      bias[j] = 0
    }
    for(var j=0; j<(1<<dimension); ++j) {
      for(var k=0; k<dimension; ++k) {
        var u = j ^ (1<<k)
        if(u > j) {
          continue
        }
        if(!(i&(1<<u)) !== !(i&(1<<j))) {
          var sign = 1
          if(i&(1<<u)) {
            denoms[k].push("v" + u + "-v" + j)
          } else {
            denoms[k].push("v" + j + "-v" + u)
            sign = -sign
          }
          if(sign < 0) {
            crossings[k].push("-v" + j + "-v" + u)
            crossingCount[k] += 2
          } else {
            crossings[k].push("v" + j + "+v" + u)
            crossingCount[k] -= 2            
          }
          totalCrossings += 1
          for(var l=0; l<dimension; ++l) {
            if(l === k) {
              continue
            }
            if(u&(1<<l)) {
              bias[l] += 1
            } else {
              bias[l] -= 1
            }
          }
        }
      }
    }
    var vertexStr = []
    for(var k=0; k<dimension; ++k) {
      if(crossings[k].length === 0) {
        vertexStr.push("d" + k + "-0.5")
      } else {
        var cStr = ""
        if(crossingCount[k] < 0) {
          cStr = crossingCount[k] + "*c"
        } else if(crossingCount[k] > 0) {
          cStr = "+" + crossingCount[k] + "*c"
        }
        var weight = 0.5 * (crossings[k].length / totalCrossings)
        var shift = 0.5 + 0.5 * (bias[k] / totalCrossings)
        vertexStr.push("d" + k + "-" + shift + "-" + weight + "*(" + crossings[k].join("+") + cStr + ")/(" + denoms[k].join("+") + ")")
        
      }
    }
    currentFunc.push("a.push([", vertexStr.join(), "]);",
      "break;")
  }
  code.push("}},")
  if(extraFuncs.length > 0) {
    currentFunc.push("}}")
  }

  //Create face function
  var faceArgs = []
  for(var i=0; i<(1<<(dimension-1)); ++i) {
    faceArgs.push("v" + i)
  }
  faceArgs.push("c0", "c1", "p0", "p1", "a", "b", "c")
  code.push("cell:function cellFunc(", faceArgs.join(), "){")

  var facets = triangulateCube(dimension-1)
  code.push("if(p0){b.push(",
    facets.map(function(f) {
      return "[" + f.map(function(v) {
        return "v" + v
      }) + "]"
    }).join(), ")}else{b.push(",
    facets.map(function(f) {
      var e = f.slice()
      e.reverse()
      return "[" + e.map(function(v) {
        return "v" + v
      }) + "]"
    }).join(),
    ")}}});function ", funcName, "(array,level){var verts=[],cells=[];contour(array,verts,cells,level);return {positions:verts,cells:cells};} return ", funcName, ";")

  for(var i=0; i<extraFuncs.length; ++i) {
    code.push(extraFuncs[i].join(""))
  }

  //Compile and link
  var proc = new Function("genContour", code.join(""))
  return proc(generateContourExtractor)
}

//1D case: Need to handle specially
function mesh1D(array, level) {
  var zc = zeroCrossings(array, level)
  var n = zc.length
  var npos = new Array(n)
  var ncel = new Array(n)
  for(var i=0; i<n; ++i) {
    npos[i] = [ zc[i] ]
    ncel[i] = [ i ]
  }
  return {
    positions: npos,
    cells: ncel
  }
}

var CACHE = {}

function surfaceNets(array,level) {
  if(array.dimension <= 0) {
    return { positions: [], cells: [] }
  } else if(array.dimension === 1) {
    return mesh1D(array, level)
  }
  var typesig = array.order.join() + "-" + array.dtype
  var proc = CACHE[typesig]
  var level = (+level) || 0.0
  if(!proc) {
    proc = CACHE[typesig] = buildSurfaceNets(array.order, array.dtype)
  }
  return proc(array,level)
}