Merge pull request #4 from neurolabusc/main

Add tt support, accelerate TRK and TCK reading

README.md

@@ -68,9 +68,11 @@ There are several important considerations regarding supporting the TRX format w
 
 The included JavaScript `bench` provides a method to evaluate performance. This benchmark is likely specific to JavaScript and so caution should be excercised in evaluating relative performance. The script will report the time to load a TRK, TCK, VTK or TRX file 10 times (it loads the tracts 11 times, and ignores the first run).
 
-The graph below shows the time load the [left inferior fronto-occipital fasciculus (IFOF) ](https://brain.labsolver.org/hcp_trk_atlas.html) from the HCP1065 Population-Averaged Tractography Atlas (Yeh, 2022). This has with 21209 streamlines and 4915166 vertices. The different formats are generated with the [tff_convert_tractogram.py](https://github.com/tee-ar-ex/trx-python) script.The benchmark was run on a MacBook laptop with a M2 Pro CPU. The ideal format would be both fast to load (to the left on the horizontal axis) and have a small file size (toward the bottom in the right axis). However, compression typically trades load time for file size. Here all data is loaded from a local solid state drive, whereas smaller files would benefit if data was loaded using a slow internet connection. The following file formats are illustrated (except where noted, both positions and indices are stored with 32-bit precision):
+The graph below shows the time load the [left inferior fronto-occipital fasciculus (IFOF) ](https://brain.labsolver.org/hcp_trk_atlas.html) from the HCP1065 Population-Averaged Tractography Atlas (Yeh, 2022). This has with 21209 streamlines and 4915166 vertices. The different formats are generated with the [tff_convert_tractogram.py](https://github.com/tee-ar-ex/trx-python) script. The benchmark was run on a MacBook laptop with a M2 Pro CPU. The ideal format would be both fast to load (to the left on the horizontal axis) and have a small file size (toward the bottom in the right axis). However, compression typically trades load time for file size. Here all data is loaded from a local solid state drive, whereas smaller files would benefit if data was loaded using a slow internet connection. The following file formats are illustrated (except where noted, both positions and indices are stored with 32-bit precision):
 
- - vtk: streamlines saved by an [undocumented](https://discourse.vtk.org/t/upcoming-changes-to-vtkcellarray/2066) extension to the [VTK legacy file format](https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf).
+ - tt: The [DSI Studio format](https://dsi-studio.labsolver.org/doc/cli_data.html) is very compact. The vertex position is stored as 1/32nd of a voxel, so it may be slightly lossy for some data.
+ - tt.gz: Gzip compressed DSI-Studio.
+ - vtk: streamlines saved by an [undocumented](https://discourse.vtk.org/t/upcoming-changes-to-vtkcellarray/2066) extension to the [VTK legacy file format](https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf) (converted with `--offsets_dtype uint32`).
  - tck: [MRtrix format](https://mrtrix.readthedocs.io/en/latest/getting_started/image_data.html#tracks-file-format-tck).
  - trk: The popular [TrackVis format](http://trackvis.org/docs/?subsect=fileformat).
  - trk.gz: Gzip compressed TrackVis.

bench.mjs

@@ -22,15 +22,21 @@ async function main() {
     let re = /(?:\.([^.]+))?$/;
     let ext = re.exec(fnm)[1];
     ext = ext.toUpperCase();
+    if (ext === 'GZ') {
+      ext = re.exec(fnm.slice(0, -3))[1]; // img.trk.gz -> img.trk
+      ext = ext.toUpperCase();
+    }
     let obj = [];
     let d = Date.now();
     let nrepeats = 11; //11 iterations, ignore first
     for (let i = 0; i < nrepeats; i++) {
         if (i == 1) d = Date.now(); //ignore first run for interpretting/disk
-        if (ext === "FIB" || ext === "VTK" || ext === "TCK" || ext === "TRK" || ext === "GZ" || ext === "ZSTD" || ext === "ZST") {
+        if (ext === "FIB" || ext === "TT" ||ext === "VTK" || ext === "TCK" || ext === "TRK" || ext === "ZSTD" || ext === "ZST") {
             const buf = fs.readFileSync(fnm);
             if (ext === "TCK")
                 obj = streamline.readTCK(new Uint8Array(buf).buffer);
+            else if (ext === "TT")
+                obj = streamline.readTT(new Uint8Array(buf).buffer);
             else if (ext === "FIB" || ext === "VTK")
                 obj = streamline.readVTK(new Uint8Array(buf).buffer);
             else

graphM2.py

@@ -4,9 +4,9 @@
 
 # Create a dataframe
 data = {
-    'format': ['vtk', 'tck', 'trk', 'trk.gz', 'trk.zst', 'trx', 'z.trx', '16.trx', '16z.trx'],
-    'time': [459, 866, 924, 3881, 4012, 159, 2763, 373, 2708],
-    'bytes': [113218694, 59236654, 59067828, 28746432, 28291562, 59152204, 28502799, 29661208, 24245538]
+    'format': ['tt', 'tt.gz', 'vtk', 'tck', 'trk', 'trk.gz', 'trk.zst', 'trx', 'z.trx', '16.trx', '16z.trx'],
+    'time': [2068, 3894, 416, 395, 573, 3483, 3827, 159, 2763, 373, 2708],
+    'bytes': [15022315, 7249360, 93473172, 59236654, 59067828, 28746432, 28291562, 59152204, 28502799, 29661208, 24245538]
 }
 
 df = pd.DataFrame(data)

streamlineIO.mjs

@@ -1,7 +1,7 @@
 //Install dependencies
 // npm install gl-matrix fflate fzstd
 
-export { readTRK, readTCK, readVTK, readTRX };
+export { readTRK, readTCK, readVTK, readTRX, readTT };
 import { mat3, mat4, vec3, vec4 } from "gl-matrix"; //for trk
 import * as fs from "fs";
 import * as fflate from "fflate";
@@ -12,6 +12,163 @@ function alert(str) { //for node.js which does not have a GUI alert
   process.exit()
 }
 
+//Read a Matlab V4 file, n.b. does not support modern versions
+//https://www.mathworks.com/help/pdf_doc/matlab/matfile_format.pdf
+function readMatV4(buffer) {
+  let len = buffer.byteLength
+  if (len < 40)
+    throw new Error("File too small to be MAT v4: bytes = " + buffer.byteLength)
+  let reader = new DataView(buffer)
+  let magic = reader.getUint16(0, true)
+  let _buffer = buffer
+  if (magic === 35615 || magic === 8075) {
+    // gzip signature 0x1F8B in little and big endian
+    const raw = fflate.decompressSync(new Uint8Array(buffer))
+    reader = new DataView(raw.buffer)
+    magic = reader.getUint16(0, true)
+    _buffer = raw.buffer
+    len = _buffer.byteLength
+  }
+  const textDecoder = new TextDecoder('utf-8')
+  let bytes = new Uint8Array(_buffer)
+  let pos = 0
+  let mat = []
+  function getTensDigit(v) {
+    return (Math.floor(v/10) % 10)
+  }
+  function readArray(tagDataType, tagBytesStart, tagBytesEnd) {
+    const byteArray = new Uint8Array(bytes.subarray(tagBytesStart, tagBytesEnd))
+    if (tagDataType === 1)
+      return new Float32Array(byteArray.buffer)
+    if (tagDataType === 2)
+      return new Int32Array(byteArray.buffer)
+    if (tagDataType === 3)
+      return new Int16Array(byteArray.buffer)
+    if (tagDataType === 4)
+      return new Uint16Array(byteArray.buffer)
+    if (tagDataType === 5)
+      return new Uint8Array(byteArray.buffer)
+    return new Float64Array(byteArray.buffer)
+  }
+  function readTag() {
+    let mtype = reader.getUint32(pos, true)
+    let mrows = reader.getUint32(pos+4, true)
+    let ncols = reader.getUint32(pos+8, true)
+    let imagf = reader.getUint32(pos+12, true)
+    let namlen = reader.getUint32(pos+16, true)
+    pos+= 20; //skip header
+    if (imagf !== 0)
+      throw new Error("Matlab V4 reader does not support imaginary numbers")
+    let tagArrayItems = mrows * ncols
+    if (tagArrayItems < 1)
+      throw new Error("mrows * ncols must be greater than one")
+    const byteArray = new Uint8Array(bytes.subarray(pos, pos+namlen))
+    let tagName = textDecoder.decode(byteArray).trim().replaceAll('\x00','')
+    let tagDataType = getTensDigit(mtype)
+    //0 double-precision (64-bit) floating-point numbers
+    //1 single-precision (32-bit) floating-point numbers
+    //2 32-bit signed integers
+    //3 16-bit signed integers
+    //4 16-bit unsigned integers
+    //5 8-bit unsigned integers
+    let tagBytesPerItem = 8
+    if ((tagDataType >= 1) && (tagDataType <= 2))
+      tagBytesPerItem = 4
+    else if ((tagDataType >= 3) && (tagDataType <= 4))
+      tagBytesPerItem = 2
+    else if (tagDataType === 5)
+      tagBytesPerItem = 1
+    else if (tagDataType !== 0)
+      throw new Error("impossible Matlab v4 datatype")
+    pos+= namlen; //skip name
+    if (mtype > 50)
+      throw new Error("Does not appear to be little-endian V4 Matlab file")
+    let posEnd = pos + (tagArrayItems * tagBytesPerItem)
+    mat[tagName] = readArray(tagDataType, pos, posEnd)
+    pos = posEnd
+  }
+  while ((pos + 20) < len)
+    readTag()
+  return mat
+} // readMatV4()
+
+// https://dsi-studio.labsolver.org/doc/cli_data.html
+// https://brain.labsolver.org/hcp_trk_atlas.html
+function readTT(buffer) {
+  let offsetPt0 = []
+  let pts = []
+  const mat = readMatV4(buffer);
+  if (!('trans_to_mni' in mat))
+    throw new Error("TT format file must have 'trans_to_mni'")
+  if (!('voxel_size' in mat))
+    throw new Error("TT format file must have 'voxel_size'")
+  if (!('track' in mat))
+    throw new Error("TT format file must have 'track'")
+  let trans_to_mni = mat4.create()
+  let m = mat.trans_to_mni
+  trans_to_mni = mat4.fromValues(m[0],m[1],m[2],m[3],  m[4],m[5],m[6],m[7],  m[8],m[9],m[10],m[11],  m[12],m[13],m[14],m[15])
+  mat4.transpose(trans_to_mni, trans_to_mni)
+  let zoomMat = mat4.create()
+  zoomMat = mat4.fromValues(1 / mat.voxel_size[0],0,0,-0.5,
+        0, 1 / mat.voxel_size[1], 0, -0.5,
+        0, 0, 1 / mat.voxel_size[2], -0.5,
+        0, 0, 0, 1)
+  mat4.transpose(zoomMat, zoomMat)
+  function parse_tt(track) {
+    let dv = new DataView(track.buffer)
+    let pos = []
+    let nvert3 = 0
+    let i = 0
+    while(i < track.length) {
+      pos.push(i)
+      let newpts = dv.getUint32(i, true)
+      i = i + newpts+13
+      nvert3 += newpts
+    }
+    offsetPt0 = new Uint32Array(pos.length+1)
+    pts = new Float32Array(nvert3)
+    let npt = 0
+    for (let i = 0; i < pos.length; i++) {
+      offsetPt0[i] = npt / 3
+      let p = pos[i]
+      let sz = dv.getUint32(p, true)/3
+      let x = dv.getInt32(p+4, true)
+      let y = dv.getInt32(p+8, true)
+      let z = dv.getInt32(p+12, true)
+      p += 16
+      pts[npt++] = x
+      pts[npt++] = y
+      pts[npt++] = z
+      for (let j = 2; j <= sz; j++) {
+          x = x + dv.getInt8(p++)
+          y = y + dv.getInt8(p++)
+          z = z + dv.getInt8(p++)
+          pts[npt++] = x
+          pts[npt++] = y
+          pts[npt++] = z
+      }
+    } //for each streamline
+    for (let i = 0; i < npt; i++)
+      pts[i] = pts[i]/32.0
+    let vox2mmMat = mat4.create()
+    mat4.mul(vox2mmMat, zoomMat, trans_to_mni)
+    let v = 0
+    for (let i = 0; i < npt / 3; i++) {
+      const pos = vec4.fromValues(pts[v], pts[v+1], pts[v+2], 1)
+      vec4.transformMat4(pos, pos, vox2mmMat)
+      pts[v++] = pos[0]
+      pts[v++] = pos[1]
+      pts[v++] = pos[2]
+    }
+    offsetPt0[pos.length] = npt / 3; //solve fence post problem, offset for final streamline
+  } // parse_tt()
+  parse_tt(mat.track)
+  return {
+    pts,
+    offsetPt0,
+  }
+} // readTT()
+
 function readTRK(buffer) {
   // http://trackvis.org/docs/?subsect=fileformat
   // http://www.tractometer.org/fiberweb/
@@ -89,7 +246,7 @@ function readTRK(buffer) {
     mat4.identity(mat);
   }
   let vox2mmMat = mat4.create();
-  mat4.mul(vox2mmMat, mat, zoomMat);
+  mat4.mul(vox2mmMat, zoomMat, mat);
   let i32 = null;
   let f32 = null;
   i32 = new Int32Array(buffer.slice(hdr_sz));
@@ -99,35 +256,36 @@ function readTRK(buffer) {
   //read and transform vertex positions
   let i = 0;
   let npt = 0;
-  let offsetPt0 = [];
-  let pts = [];
+  //over-provision offset array to store number of segments
+  let offsetPt0 = new Uint32Array(i32.length);
+  let noffset = 0;;
+  //over-provision points array to store vertex positions
+  let npt3 = 0;
+  let pts = new Float32Array(i32.length);
   while (i < ntracks) {
     let n_pts = i32[i];
     i = i + 1; // read 1 32-bit integer for number of points in this streamline
-    offsetPt0.push(npt); //index of first vertex in this streamline
+    offsetPt0[noffset++] = npt; //index of first vertex in this streamline
     for (let j = 0; j < n_pts; j++) {
       let ptx = f32[i + 0];
       let pty = f32[i + 1];
       let ptz = f32[i + 2];
       i += 3; //read 3 32-bit floats for XYZ position
-      pts.push(
+      pts[npt3++] =
         ptx * vox2mmMat[0] +
           pty * vox2mmMat[1] +
           ptz * vox2mmMat[2] +
-          vox2mmMat[3]
-      );
-      pts.push(
+          vox2mmMat[3];
+      pts[npt3++] =
         ptx * vox2mmMat[4] +
           pty * vox2mmMat[5] +
           ptz * vox2mmMat[6] +
-          vox2mmMat[7]
-      );
-      pts.push(
+          vox2mmMat[7];
+      pts[npt3++] =
         ptx * vox2mmMat[8] +
           pty * vox2mmMat[9] +
           ptz * vox2mmMat[10] +
-          vox2mmMat[11]
-      );
+          vox2mmMat[11];
       if (n_scalars > 0) {
         for (let s = 0; s < n_scalars; s++) {
           dpv[s].vals.push(f32[i]);
@@ -143,7 +301,11 @@ function readTRK(buffer) {
       }
     }
   } //for each streamline: while i < n_count
-  offsetPt0.push(npt); //add 'first index' as if one more line was added (fence post problem)
+  //add 'first index' as if one more line was added (fence post problem)
+  offsetPt0[noffset++] = npt;
+  //resize offset/vertex arrays that were initially over-provisioned
+  pts = pts.slice(0, npt3);
+  offsetPt0 = offsetPt0.slice(0, noffset); 
   return {
     pts,
     offsetPt0,
@@ -176,9 +338,13 @@ function readTCK(buffer) {
   let reader = new DataView(buffer);
   //read and transform vertex positions
   let npt = 0;
-  let offsetPt0 = [];
-  offsetPt0.push(npt); //1st streamline starts at 0
-  let pts = [];
+  //over-provision offset array to store number of segments
+  let offsetPt0 = new Uint32Array(len / 4);
+  let noffset = 0;
+  //over-provision points array to store vertex positions
+  let npt3 = 0;
+  let pts = new Float32Array(len / 4);
+  offsetPt0[0] = 0; //1st streamline starts at 0
   while (pos + 12 < len) {
     let ptx = reader.getFloat32(pos, true);
     pos += 4;
@@ -188,17 +354,20 @@ function readTCK(buffer) {
     pos += 4;
     if (!isFinite(ptx)) {
       //both NaN and Inifinity are not finite
-      offsetPt0.push(npt);
+      offsetPt0[noffset++] = npt;
       if (!isNaN(ptx))
         //terminate if infinity
         break;
     } else {
-      pts.push(ptx);
-      pts.push(pty);
-      pts.push(ptz);
+      pts[npt3++] = ptx;
+      pts[npt3++] = pty;
+      pts[npt3++] = ptz;
       npt++;
     }
   }
+  //resize offset/vertex arrays that were initially over-provisioned
+  pts = pts.slice(0, npt3);
+  offsetPt0 = offsetPt0.slice(0, noffset); 
   return {
     pts,
     offsetPt0,