src/parser/dsn6-parser.js

/**
 * @file Dsn6 Parser
 * @author Alexander Rose <alexander.rose@weirdbyte.de>
 * @private
 */
 
import { Matrix4 } from '../../lib/three.es6.js'
 
import { Debug, Log, ParserRegistry } from '../globals.js'
import { ensureBuffer } from '../utils.js'
import { degToRad } from '../math/math-utils.js'
import VolumeParser from './volume-parser.js'
 
class Dsn6Parser extends VolumeParser {
  get type () { return 'dsn6' }
  get isBinary () { return true }
 
  _parse () {
    // DSN6 http://www.uoxray.uoregon.edu/tnt/manual/node104.html
    // BRIX http://svn.cgl.ucsf.edu/svn/chimera/trunk/libs/VolumeData/dsn6/brix-1.html
 
    if (Debug) Log.time('Dsn6Parser._parse ' + this.name)
 
    const v = this.volume
    const header = {}
    let divisor, summand
 
    const bin = ensureBuffer(this.streamer.data)
    const intView = new Int16Array(bin)
    const byteView = new Uint8Array(bin)
    const brixStr = String.fromCharCode.apply(null, byteView.subarray(0, 512))
 
    if (brixStr.startsWith(':-)')) {
      header.xStart = parseInt(brixStr.substr(10, 5))  // NXSTART
      header.yStart = parseInt(brixStr.substr(15, 5))
      header.zStart = parseInt(brixStr.substr(20, 5))
 
      header.xExtent = parseInt(brixStr.substr(32, 5))  // NX
      header.yExtent = parseInt(brixStr.substr(38, 5))
      header.zExtent = parseInt(brixStr.substr(42, 5))
 
      header.xRate = parseInt(brixStr.substr(52, 5))  // MX
      header.yRate = parseInt(brixStr.substr(58, 5))
      header.zRate = parseInt(brixStr.substr(62, 5))
 
      header.xlen = parseFloat(brixStr.substr(73, 10)) * this.voxelSize
      header.ylen = parseFloat(brixStr.substr(83, 10)) * this.voxelSize
      header.zlen = parseFloat(brixStr.substr(93, 10)) * this.voxelSize
 
      header.alpha = parseFloat(brixStr.substr(103, 10))
      header.beta = parseFloat(brixStr.substr(113, 10))
      header.gamma = parseFloat(brixStr.substr(123, 10))
 
      divisor = parseFloat(brixStr.substr(138, 12)) / 100
      summand = parseInt(brixStr.substr(155, 8))
 
      header.sigma = parseFloat(brixStr.substr(170, 12)) * 100
    } else {
      // swap byte order when big endian
      if (intView[ 18 ] !== 100) {
        for (let i = 0, n = intView.length; i < n; ++i) {
          const val = intView[ i ]
          intView[ i ] = ((val & 0xff) << 8) | ((val >> 8) & 0xff)
        }
      }
 
      header.xStart = intView[ 0 ]  // NXSTART
      header.yStart = intView[ 1 ]
      header.zStart = intView[ 2 ]
 
      header.xExtent = intView[ 3 ]  // NX
      header.yExtent = intView[ 4 ]
      header.zExtent = intView[ 5 ]
 
      header.xRate = intView[ 6 ]  // MX
      header.yRate = intView[ 7 ]
      header.zRate = intView[ 8 ]
 
      const factor = 1 / intView[ 17 ]
      const scalingFactor = factor * this.voxelSize
 
      header.xlen = intView[ 9 ] * scalingFactor
      header.ylen = intView[ 10 ] * scalingFactor
      header.zlen = intView[ 11 ] * scalingFactor
 
      header.alpha = intView[ 12 ] * factor
      header.beta = intView[ 13 ] * factor
      header.gamma = intView[ 14 ] * factor
 
      divisor = intView[ 15 ] / 100
      summand = intView[ 16 ]
      header.gamma = intView[ 14 ] * factor
    }
 
    v.header = header
 
    if (Debug) Log.log(header, divisor, summand)
 
    const data = new Float32Array(
      header.xExtent * header.yExtent * header.zExtent
    )
 
    let offset = 512
    const xBlocks = Math.ceil(header.xExtent / 8)
    const yBlocks = Math.ceil(header.yExtent / 8)
    const zBlocks = Math.ceil(header.zExtent / 8)
 
    // loop over blocks
    for (var zz = 0; zz < zBlocks; ++zz) {
      for (var yy = 0; yy < yBlocks; ++yy) {
        for (var xx = 0; xx < xBlocks; ++xx) {
          // loop inside block
          for (var k = 0; k < 8; ++k) {
            var z = 8 * zz + k
            for (var j = 0; j < 8; ++j) {
              var y = 8 * yy + j
              for (var i = 0; i < 8; ++i) {
                var x = 8 * xx + i
 
                // check if remaining slice-part contains data
                if (x < header.xExtent && y < header.yExtent && z < header.zExtent) {
                  var idx = ((((x * header.yExtent) + y) * header.zExtent) + z)
                  data[ idx ] = (byteView[ offset ] - summand) / divisor
                  ++offset
                } else {
                  offset += 8 - i
                  break
                }
              }
            }
          }
        }
      }
    }
 
    v.setData(data, header.zExtent, header.yExtent, header.xExtent)
    if (header.sigma) {
      v.setStats(undefined, undefined, undefined, header.sigma)
    }
 
    if (Debug) Log.timeEnd('Dsn6Parser._parse ' + this.name)
  }
 
  getMatrix () {
    const h = this.volume.header
 
    const basisX = [
      h.xlen,
      0,
      0
    ]
 
    const basisY = [
      h.ylen * Math.cos(Math.PI / 180.0 * h.gamma),
      h.ylen * Math.sin(Math.PI / 180.0 * h.gamma),
      0
    ]
 
    const basisZ = [
      h.zlen * Math.cos(Math.PI / 180.0 * h.beta),
      h.zlen * (
        Math.cos(Math.PI / 180.0 * h.alpha) -
        Math.cos(Math.PI / 180.0 * h.gamma) *
        Math.cos(Math.PI / 180.0 * h.beta)
      ) / Math.sin(Math.PI / 180.0 * h.gamma),
      0
    ]
    basisZ[ 2 ] = Math.sqrt(
      h.zlen * h.zlen * Math.sin(Math.PI / 180.0 * h.beta) *
      Math.sin(Math.PI / 180.0 * h.beta) - basisZ[ 1 ] * basisZ[ 1 ]
    )
 
    const basis = [ 0, basisX, basisY, basisZ ]
    const nxyz = [ 0, h.xRate, h.yRate, h.zRate ]
    const mapcrs = [ 0, 1, 2, 3 ]
 
    const matrix = new Matrix4()
 
    matrix.set(
      basis[ mapcrs[1] ][0] / nxyz[ mapcrs[1] ],
      basis[ mapcrs[2] ][0] / nxyz[ mapcrs[2] ],
      basis[ mapcrs[3] ][0] / nxyz[ mapcrs[3] ],
      0,
      basis[ mapcrs[1] ][1] / nxyz[ mapcrs[1] ],
      basis[ mapcrs[2] ][1] / nxyz[ mapcrs[2] ],
      basis[ mapcrs[3] ][1] / nxyz[ mapcrs[3] ],
      0,
      basis[ mapcrs[1] ][2] / nxyz[ mapcrs[1] ],
      basis[ mapcrs[2] ][2] / nxyz[ mapcrs[2] ],
      basis[ mapcrs[3] ][2] / nxyz[ mapcrs[3] ],
      0,
      0, 0, 0, 1
    )
 
    matrix.multiply(
      new Matrix4().makeRotationY(degToRad(90))
    )
 
    matrix.multiply(new Matrix4().makeTranslation(
      -h.zStart, h.yStart, h.xStart
    ))
 
    matrix.multiply(new Matrix4().makeScale(
      -1, 1, 1
    ))
 
    return matrix
  }
}
 
ParserRegistry.add('dsn6', Dsn6Parser)
ParserRegistry.add('brix', Dsn6Parser)
 
export default Dsn6Parser