diff --git a/M2.png b/M2.png index e795093..3c22d97 100644 Binary files a/M2.png and b/M2.png differ diff --git a/README.md b/README.md index 99d1315..d1b1776 100644 --- a/README.md +++ b/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. diff --git a/bench.mjs b/bench.mjs index 36e4084..4c51cee 100644 --- a/bench.mjs +++ b/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 diff --git a/graphM2.py b/graphM2.py index aae754d..1860266 100755 --- a/graphM2.py +++ b/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) diff --git a/streamlineIO.mjs b/streamlineIO.mjs index cfe1413..894311d 100644 --- a/streamlineIO.mjs +++ b/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,