From 0a399342037f9862066965ef0d5d99d425f8552a Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 17 Jun 2026 14:10:14 -0400 Subject: [PATCH 1/4] Modification for larges files and improve support across dtype --- .gitignore | 5 + streamlineIO.mjs | 243 +++++++++++++++++++++++++++++++++++++---------- 2 files changed, 199 insertions(+), 49 deletions(-) diff --git a/.gitignore b/.gitignore index a1c2a23..2474c5c 100644 --- a/.gitignore +++ b/.gitignore @@ -21,3 +21,8 @@ # virtual machine crash logs, see http://www.java.com/en/download/help/error_hotspot.xml hs_err_pid* + +# JS modules +node_modules/ +package-lock.json +package.json diff --git a/streamlineIO.mjs b/streamlineIO.mjs index 894311d..744f8e9 100644 --- a/streamlineIO.mjs +++ b/streamlineIO.mjs @@ -247,30 +247,36 @@ function readTRK(buffer) { } let vox2mmMat = mat4.create(); mat4.mul(vox2mmMat, zoomMat, mat); - let i32 = null; - let f32 = null; - i32 = new Int32Array(buffer.slice(hdr_sz)); - f32 = new Float32Array(i32.buffer); + let dataView = new DataView(buffer, hdr_sz); + let totalWords = dataView.byteLength / 4; - let ntracks = i32.length; - //read and transform vertex positions - let i = 0; + let num_streamlines = 0; + let num_points = 0; + let w = 0; + while (w < totalWords) { + let n_pts = dataView.getInt32(w * 4, true); + num_streamlines++; + num_points += n_pts; + w += 1 + n_pts * (3 + n_scalars) + n_properties; + } + + let offsetPt0 = new Uint32Array(num_streamlines + 1); + let pts = new Float32Array(num_points * 3); + + w = 0; let npt = 0; - //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 noffset = 0; 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 + + while (w < totalWords) { + let n_pts = dataView.getInt32(w * 4, true); + w = w + 1; // read 1 32-bit integer for number of points 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 + let ptx = dataView.getFloat32(w * 4, true); + let pty = dataView.getFloat32((w + 1) * 4, true); + let ptz = dataView.getFloat32((w + 2) * 4, true); + w += 3; //read 3 32-bit floats for XYZ position pts[npt3++] = ptx * vox2mmMat[0] + pty * vox2mmMat[1] + @@ -288,24 +294,21 @@ function readTRK(buffer) { vox2mmMat[11]; if (n_scalars > 0) { for (let s = 0; s < n_scalars; s++) { - dpv[s].vals.push(f32[i]); - i++; + dpv[s].vals.push(dataView.getFloat32(w * 4, true)); + w++; } } npt++; } // for j: each point in streamline if (n_properties > 0) { for (let j = 0; j < n_properties; j++) { - dps[j].vals.push(f32[i]); - i++; + dps[j].vals.push(dataView.getFloat32(w * 4, true)); + w++; } } - } //for each streamline: while i < n_count + } //for each streamline: while w < totalWords //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, @@ -345,6 +348,7 @@ function readTCK(buffer) { let npt3 = 0; let pts = new Float32Array(len / 4); offsetPt0[0] = 0; //1st streamline starts at 0 + noffset = 1; while (pos + 12 < len) { let ptx = reader.getFloat32(pos, true); pos += 4; @@ -666,6 +670,87 @@ function readVTK (buffer) { }; }; // readVTK() +function getZip64OriginalSizes(zipData) { + // Parse the ZIP central directory to get correct uncompressed sizes. + // fflate ZIP64 bug: file.originalSize in the filter callback returns 0xFFFFFFFF + // (the 32-bit sentinel) instead of reading the real size from the ZIP64 + // extended information extra field (tag 0x0001) in the central directory. + // + // Supports ZIP32 (entries < 4 GB, CD < 4 GB) and ZIP64 (entries or CD + // offset up to 2^53 bytes, the JavaScript safe-integer limit). + let u8; + if (zipData instanceof ArrayBuffer) { + u8 = new Uint8Array(zipData); + } else { + u8 = new Uint8Array(zipData.buffer, zipData.byteOffset, zipData.byteLength); + } + const dv = new DataView(u8.buffer, u8.byteOffset, u8.byteLength); + const n = u8.byteLength; + const sizes = {}; + + // Read a 64-bit little-endian uint as a JS number (safe up to 2^53) + function getUint64(offset) { + const lo = dv.getUint32(offset, true); + const hi = dv.getUint32(offset + 4, true); + return hi * 0x100000000 + lo; + } + + // Find End of Central Directory record (PK\x05\x06) scanning backwards. + // Stop 22 bytes from end (minimum EOCD size); allow up to 64 KB ZIP comment. + let eocd = -1; + for (let i = n - 22; i >= Math.max(0, n - 65558); i--) { + if (dv.getUint32(i, true) === 0x06054b50) { eocd = i; break; } + } + if (eocd === -1) return sizes; + + let cdCount = dv.getUint16(eocd + 8, true); + let cdOffset = dv.getUint32(eocd + 16, true); + + // Check for ZIP64 EOCD locator (PK\x06\x07), which must sit exactly 20 + // bytes before the EOCD when no ZIP comment is present, but the ZIP spec + // only guarantees it precedes the EOCD — scan the last 64 KB for it. + for (let i = eocd - 20; i >= Math.max(0, eocd - 65558); i--) { + if (dv.getUint32(i, true) === 0x07064b50) { // PK\x06\x07 + const eocd64off = getUint64(i + 8); // offset of ZIP64 EOCD + if (eocd64off + 56 <= n && dv.getUint32(eocd64off, true) === 0x06064b50) { + cdCount = getUint64(eocd64off + 32); // 8-byte total entry count + cdOffset = getUint64(eocd64off + 48); // 8-byte CD start offset + } + break; + } + } + + // Parse central directory records + let pos = cdOffset; + for (let i = 0; i < cdCount; i++) { + if (pos + 46 > n || dv.getUint32(pos, true) !== 0x02014b50) break; // PK\x01\x02 + let origSize = dv.getUint32(pos + 24, true); + const fnLen = dv.getUint16(pos + 28, true); + const exLen = dv.getUint16(pos + 30, true); + const cmLen = dv.getUint16(pos + 32, true); + const fname = new TextDecoder().decode(u8.subarray(pos + 46, pos + 46 + fnLen)); + // If 0xFFFFFFFF, read real size from ZIP64 extended info extra field (tag 0x0001). + // The ZIP64 extra field encodes: origSize (8), compSize (8), localOffset (8), + // diskStart (4) — but only the fields that were 0xFFFFFFFF in the CD are present. + if (origSize === 0xFFFFFFFF) { + let ep = pos + 46 + fnLen; + const epEnd = ep + exLen; + while (ep + 4 <= epEnd) { + const tag = dv.getUint16(ep, true); + const sz = dv.getUint16(ep + 2, true); + if (tag === 0x0001 && sz >= 8) { + origSize = getUint64(ep + 4); // first 8-byte field is original size + break; + } + ep += 4 + sz; + } + } + sizes[fname] = origSize; + pos += 46 + fnLen + exLen + cmLen; + } + return sizes; +} // getZip64OriginalSizes() + async function readTRX(url, urlIsLocalFile = false) { //Javascript does not support float16, so we convert to float32 //https://stackoverflow.com/questions/5678432/decompressing-half-precision-floats-in-javascript @@ -690,19 +775,50 @@ async function readTRX(url, urlIsLocalFile = false) { let npt = 0; let pts = []; let offsetPt0 = []; - let dpg = []; - let dps = []; let dpv = []; + let dps = []; + let dpg = []; + let groups = []; let header = []; let isOverflowUint64 = false; + let positions_dtype = "float32"; let data = []; + function getAlignedArray(constructor, dataArray) { + const bytes = constructor.BYTES_PER_ELEMENT; + if (dataArray.byteOffset % bytes === 0) { + return new constructor(dataArray.buffer, dataArray.byteOffset, dataArray.byteLength / bytes); + } else { + return new constructor(dataArray.slice().buffer); + } + } if (urlIsLocalFile) { - data = fs.readFileSync(url); + const stats = fs.statSync(url); + if (stats.size >= 2 * 1024 * 1024 * 1024) { + const size = stats.size; + const arrayBuffer = new ArrayBuffer(size); + const uint8Array = new Uint8Array(arrayBuffer); + const fd = fs.openSync(url, 'r'); + const chunkSize = 512 * 1024 * 1024; + let offset = 0; + while (offset < size) { + const bytesToRead = Math.min(chunkSize, size - offset); + const bytesRead = fs.readSync(fd, uint8Array, offset, bytesToRead, offset); + if (bytesRead === 0) break; + offset += bytesRead; + } + fs.closeSync(fd); + data = Buffer.from(arrayBuffer); + } else { + data = fs.readFileSync(url); + } } else { let response = await fetch(url); if (!response.ok) throw Error(response.statusText); data = await response.arrayBuffer(); } + // Parse the ZIP central directory ourselves: fflate's file.originalSize + // returns 0xFFFFFFFF for ZIP64 entries instead of reading the ZIP64 extra field. + const sizeMap = getZip64OriginalSizes(data); const decompressed = fflate.unzipSync(data, { filter(file) { return file.originalSize > 0; @@ -717,22 +833,29 @@ async function readTRX(url, urlIsLocalFile = false) { let tag = fname.split(".")[0]; // "positions.3.float16 -> "positions" //todo: should tags be censored for invalide characters: https://stackoverflow.com/questions/8676011/which-characters-are-valid-invalid-in-a-json-key-name let data = decompressed[keys[i]]; + // Truncate to the size recorded in the ZIP central directory to work + // around a fflate ZIP64 bug where extra bytes from subsequent entries + // are appended to the decompressed output. + const correctByteLength = sizeMap[keys[i]]; + if (correctByteLength !== undefined && data.byteLength > correctByteLength) { + data = data.slice(0, correctByteLength); + } if (fname.includes("header.json")) { let jsonString = new TextDecoder().decode(data); - header = JSON.parse(jsonString); + // Strip UTF-8 BOM if present (some tools prepend \uFEFF) + if (jsonString.charCodeAt(0) === 0xFEFF) jsonString = jsonString.slice(1); + header = JSON.parse(jsonString.trim()); continue; } + + //next read arrays for all possible datatypes: int8/16/32/64 uint8/16/32/64 float16/32/64 let nval = 0; let vals = []; if (fname.endsWith(".uint64") || fname.endsWith(".int64")) { - //javascript does not have 64-bit integers! read lower 32-bits - //note for signed int64 we only read unsigned bytes - //for both signed and unsigned, generate an error if any value is out of bounds - //one alternative might be to convert to 64-bit double that has a flintmax of 2^53. nval = data.length / 8; //8 bytes per 64bit input vals = new Uint32Array(nval); - var u32 = new Uint32Array(data.buffer); + const u32 = getAlignedArray(Uint32Array, data); let j = 0; for (let i = 0; i < nval; i++) { vals[i] = u32[j]; @@ -740,35 +863,44 @@ async function readTRX(url, urlIsLocalFile = false) { j += 2; } } else if (fname.endsWith(".uint32")) { - vals = new Uint32Array(data.buffer); + vals = getAlignedArray(Uint32Array, data); } else if (fname.endsWith(".uint16")) { - vals = new Uint16Array(data.buffer); + vals = getAlignedArray(Uint16Array, data); } else if (fname.endsWith(".uint8")) { - vals = new Uint8Array(data.buffer); + vals = new Uint8Array(data.buffer, data.byteOffset, data.byteLength); } else if (fname.endsWith(".int32")) { - vals = new Int32Array(data.buffer); + vals = getAlignedArray(Int32Array, data); } else if (fname.endsWith(".int16")) { - vals = new Int16Array(data.buffer); + vals = getAlignedArray(Int16Array, data); } else if (fname.endsWith(".int8")) { - vals = new Int8Array(data.buffer); + vals = new Int8Array(data.buffer, data.byteOffset, data.byteLength); } else if (fname.endsWith(".float64")) { - vals = new Float64Array(data.buffer); + vals = getAlignedArray(Float64Array, data); } else if (fname.endsWith(".float32")) { - vals = new Float32Array(data.buffer); + vals = getAlignedArray(Float32Array, data); } else if (fname.endsWith(".float16")) { - //javascript does not have 16-bit floats! Convert to 32-bits nval = data.length / 2; //2 bytes per 16bit input vals = new Float32Array(nval); - var u16 = new Uint16Array(data.buffer); + const u16 = getAlignedArray(Uint16Array, data); const lut = new Float32Array(65536) for (let i = 0; i < 65536; i++) lut[i] = decodeFloat16(i) for (let i = 0; i < nval; i++) vals[i] = lut[u16[i]] } else continue; //not a data array nval = vals.length; //next: read data_per_group - if (pname.includes("dpg")) { + if (parts[0] === "dpg") { dpg.push({ + id: parts[1] + "/" + tag, // e.g. "AF_R/volume" + fname: parts.slice(1).join("/"), // e.g. "AF_R/volume.uint32" + vals: vals.slice(), + }); + continue; + } + //next: read groups + if (pname === "groups") { + groups.push({ id: tag, + fname: fname, vals: vals.slice(), }); continue; @@ -777,6 +909,7 @@ async function readTRX(url, urlIsLocalFile = false) { if (pname.includes("dpv")) { dpv.push({ id: tag, + fname: fname, vals: vals.slice(), }); continue; @@ -785,6 +918,7 @@ async function readTRX(url, urlIsLocalFile = false) { if (pname.includes("dps")) { dps.push({ id: tag, + fname: fname, vals: vals.slice(), }); continue; @@ -800,18 +934,29 @@ async function readTRX(url, urlIsLocalFile = false) { if (fname.startsWith("positions.3.")) { npt = nval; //4 bytes per 32bit input pts = vals.slice(); + if (fname.endsWith(".float64")) positions_dtype = "float64"; + else if (fname.endsWith(".float16")) positions_dtype = "float16"; + else positions_dtype = "float32"; } } if (noff === 0 || npt === 0) alert("Failure reading TRX format"); if (isOverflowUint64) alert("Too many vertices: JavaScript does not support 64 bit integers"); - offsetPt0[noff] = npt / 3; //solve fence post problem, offset for final streamline + + if (offsetPt0[noff - 1] === npt / 3) { + offsetPt0 = offsetPt0.slice(0, noff); + } else { + offsetPt0[noff] = npt / 3; //solve fence post problem, offset for final streamline + } + return { pts, offsetPt0, - dpg, dps, dpv, + dpg, + groups, header, + positions_dtype, }; }; // readTRX() From fa0fd77aaaed69a0e830fd74a32883908b03e8e4 Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 22 Jun 2026 10:10:39 -0400 Subject: [PATCH 2/4] resolve TRK coordinate mismatch by applying inverse affine translation to vertices during export --- streamlineIO.mjs | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/streamlineIO.mjs b/streamlineIO.mjs index 744f8e9..9a0816f 100644 --- a/streamlineIO.mjs +++ b/streamlineIO.mjs @@ -309,13 +309,25 @@ function readTRK(buffer) { } //for each streamline: while w < totalWords //add 'first index' as if one more line was added (fence post problem) offsetPt0[noffset++] = npt; + let header = { + DIMENSIONS: [reader.getInt16(6, true), reader.getInt16(8, true), reader.getInt16(10, true)], + VOXEL_SIZES: [voxel_sizeX, voxel_sizeY, voxel_sizeZ], + VOXEL_TO_RASMM: [ + [mat[0], mat[1], mat[2], mat[3]], + [mat[4], mat[5], mat[6], mat[7]], + [mat[8], mat[9], mat[10], mat[11]], + [mat[12], mat[13], mat[14], mat[15]] + ] + }; + return { pts, offsetPt0, dps, dpv, + header }; -}; //readTRK() +} // readTRK() function readTCK(buffer) { //https://mrtrix.readthedocs.io/en/latest/getting_started/image_data.html#tracks-file-format-tck From c4d3858925865cc55ebdaab2f4338f2132f46b3a Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 22 Jun 2026 10:32:17 -0400 Subject: [PATCH 3/4] Centralize legacy streamline IO functions and implement inverse affine corrections for valid TRK export --- streamlineIO.mjs | 387 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 386 insertions(+), 1 deletion(-) diff --git a/streamlineIO.mjs b/streamlineIO.mjs index 9a0816f..71c5525 100644 --- a/streamlineIO.mjs +++ b/streamlineIO.mjs @@ -1,7 +1,7 @@ //Install dependencies // npm install gl-matrix fflate fzstd -export { readTRK, readTCK, readVTK, readTRX, readTT }; +export { readTRK, readTCK, readVTK, readTRX, readTT, saveTCK, saveTRK, saveVTK, saveTRX }; import { mat3, mat4, vec3, vec4 } from "gl-matrix"; //for trk import * as fs from "fs"; import * as fflate from "fflate"; @@ -972,3 +972,388 @@ async function readTRX(url, urlIsLocalFile = false) { positions_dtype, }; }; // readTRX() +// Fast float32 to float16 conversion +function encodeFloat16(val) { + const floatView = new Float32Array(1); + const int32View = new Int32Array(floatView.buffer); + floatView[0] = val; + const f = int32View[0]; + + const sign = (f >> 16) & 0x8000; + let exponent = ((f >> 23) & 0xff) - 127; + let mantissa = f & 0x007fffff; + + if (exponent <= -15) { + if (exponent < -24) { + return sign; // underflow + } + mantissa = (mantissa | 0x00800000) >> (-14 - exponent); + return sign | (mantissa >> 13); + } else if (exponent >= 16) { + return sign | 0x7c00; // overflow to infinity + } + + return sign | ((exponent + 15) << 10) | (mantissa >> 13); +} + +function float32ToFloat16(float32Array) { + const out = new Uint16Array(float32Array.length); + for (let i = 0; i < float32Array.length; i++) { + out[i] = encodeFloat16(float32Array[i]); + } + return out; +} + +function buildTckHeader(numStreamlines) { + let offset = 80; + while (true) { + let h = `mrtrix tracks\ncount: ${String(numStreamlines).padStart(10, '0')}\ndatatype: Float32LE\nfile: . ${offset}\nEND\n`; + if (h.length <= offset) { + return h.padEnd(offset, ' '); + } + offset = h.length; + } +} + +function saveTCK(filepath, obj) { + const fd = fs.openSync(filepath, 'w'); + const numStreamlines = obj.offsetPt0.length - 1; + const header = buildTckHeader(numStreamlines); + fs.writeSync(fd, header); + + const chunkSize = 16 * 1024 * 1024; // 16MB buffer + const buf = new ArrayBuffer(chunkSize); + const view = new DataView(buf); + const u8View = new Uint8Array(buf); + + let bufOffset = 0; + const pts = obj.pts; + const offsets = obj.offsetPt0; + + for (let i = 0; i < numStreamlines; i++) { + const start = offsets[i]; + const end = offsets[i+1]; + + for (let j = start; j < end; j++) { + if (bufOffset + 12 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setFloat32(bufOffset, pts[j*3], true); + view.setFloat32(bufOffset + 4, pts[j*3 + 1], true); + view.setFloat32(bufOffset + 8, pts[j*3 + 2], true); + bufOffset += 12; + } + if (bufOffset + 12 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setFloat32(bufOffset, NaN, true); + view.setFloat32(bufOffset + 4, NaN, true); + view.setFloat32(bufOffset + 8, NaN, true); + bufOffset += 12; + } + + if (bufOffset + 12 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setFloat32(bufOffset, Infinity, true); + view.setFloat32(bufOffset + 4, Infinity, true); + view.setFloat32(bufOffset + 8, Infinity, true); + bufOffset += 12; + + if (bufOffset > 0) { + fs.writeSync(fd, u8View, 0, bufOffset); + } + fs.closeSync(fd); +} + +function saveTRK(filepath, obj, originalFilename) { + const fd = fs.openSync(filepath, 'w'); + const headerBytes = new Uint8Array(1000); + const view = new DataView(headerBytes.buffer); + + headerBytes.set([84, 82, 65, 67, 75], 0); + + let dim = [256, 256, 256]; + let voxelSize = [1, 1, 1]; + let voxToRas = [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1] + ]; + + if (obj.header) { + if (obj.header.DIMENSIONS) dim = obj.header.DIMENSIONS; + if (obj.header.VOXEL_TO_RASMM) { + voxToRas = obj.header.VOXEL_TO_RASMM; + voxelSize = [ + Math.sqrt(voxToRas[0][0]**2 + voxToRas[1][0]**2 + voxToRas[2][0]**2), + Math.sqrt(voxToRas[0][1]**2 + voxToRas[1][1]**2 + voxToRas[2][1]**2), + Math.sqrt(voxToRas[0][2]**2 + voxToRas[1][2]**2 + voxToRas[2][2]**2) + ]; + } + } + + view.setInt16(6, dim[0], true); + view.setInt16(8, dim[1], true); + view.setInt16(10, dim[2], true); + + view.setFloat32(12, voxelSize[0], true); + view.setFloat32(16, voxelSize[1], true); + view.setFloat32(20, voxelSize[2], true); + + let off = 440; + for (let r = 0; r < 4; r++) { + for (let c = 0; c < 4; c++) { + view.setFloat32(off, voxToRas[r][c], true); + off += 4; + } + } + + headerBytes.set([82, 65, 83, 0], 948); + + let invMat = null; + if (obj.header && obj.header.VOXEL_TO_RASMM) { + let mat = mat4.fromValues( + voxToRas[0][0], voxToRas[1][0], voxToRas[2][0], voxToRas[3][0], + voxToRas[0][1], voxToRas[1][1], voxToRas[2][1], voxToRas[3][1], + voxToRas[0][2], voxToRas[1][2], voxToRas[2][2], voxToRas[3][2], + voxToRas[0][3], voxToRas[1][3], voxToRas[2][3], voxToRas[3][3] + ); + invMat = mat4.create(); + mat4.invert(invMat, mat); + } + + const numStreamlines = obj.offsetPt0.length - 1; + view.setInt32(988, numStreamlines, true); + view.setInt32(992, 2, true); + view.setInt32(996, 1000, true); + + fs.writeSync(fd, headerBytes); + + const chunkSize = 16 * 1024 * 1024; + const buf = new ArrayBuffer(chunkSize); + const payloadView = new DataView(buf); + const u8View = new Uint8Array(buf); + + let bufOffset = 0; + const pts = obj.pts; + const offsets = obj.offsetPt0; + + for (let i = 0; i < numStreamlines; i++) { + const start = offsets[i]; + const end = offsets[i+1]; + const n_pts = end - start; + + if (bufOffset + 4 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + payloadView.setInt32(bufOffset, n_pts, true); + bufOffset += 4; + + for (let j = start; j < end; j++) { + if (bufOffset + 12 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + + let x = pts[j*3]; + let y = pts[j*3 + 1]; + let z = pts[j*3 + 2]; + let vx = x, vy = y, vz = z; + + if (invMat) { + let cx = x * invMat[0] + y * invMat[4] + z * invMat[8] + invMat[12]; + let cy = x * invMat[1] + y * invMat[5] + z * invMat[9] + invMat[13]; + let cz = x * invMat[2] + y * invMat[6] + z * invMat[10] + invMat[14]; + + vx = (cx + 0.5) * voxelSize[0]; + vy = (cy + 0.5) * voxelSize[1]; + vz = (cz + 0.5) * voxelSize[2]; + } else { + vx = x + 0.5; + vy = y + 0.5; + vz = z + 0.5; + } + + payloadView.setFloat32(bufOffset, vx, true); + payloadView.setFloat32(bufOffset + 4, vy, true); + payloadView.setFloat32(bufOffset + 8, vz, true); + bufOffset += 12; + } + } + + if (bufOffset > 0) { + fs.writeSync(fd, u8View, 0, bufOffset); + } + fs.closeSync(fd); +} + +function saveVTK(filepath, obj) { + const fd = fs.openSync(filepath, 'w'); + const numStreamlines = obj.offsetPt0.length - 1; + const numPoints = obj.pts.length / 3; + + const header = `# vtk DataFile Version 3.0\nvtk output\nBINARY\nDATASET POLYDATA\nPOINTS ${numPoints} float\n`; + fs.writeSync(fd, header); + + const pts = obj.pts; + const offsets = obj.offsetPt0; + + const chunkSize = 16 * 1024 * 1024; + const buf = new ArrayBuffer(chunkSize); + const view = new DataView(buf); + const u8View = new Uint8Array(buf); + + let bufOffset = 0; + for (let i = 0; i < numPoints; i++) { + if (bufOffset + 12 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setFloat32(bufOffset, pts[i*3], false); // Big endian + view.setFloat32(bufOffset + 4, pts[i*3 + 1], false); + view.setFloat32(bufOffset + 8, pts[i*3 + 2], false); + bufOffset += 12; + } + if (bufOffset > 0) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + + const cellArraySize = numStreamlines + numPoints; + const linesHeader = `LINES ${numStreamlines} ${cellArraySize}\n`; + fs.writeSync(fd, linesHeader); + + for (let i = 0; i < numStreamlines; i++) { + const start = offsets[i]; + const end = offsets[i+1]; + const n_pts = end - start; + + if (bufOffset + 4 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setInt32(bufOffset, n_pts, false); // Big endian + bufOffset += 4; + + for (let j = start; j < end; j++) { + if (bufOffset + 4 > chunkSize) { + fs.writeSync(fd, u8View, 0, bufOffset); + bufOffset = 0; + } + view.setInt32(bufOffset, j, false); + bufOffset += 4; + } + } + if (bufOffset > 0) { + fs.writeSync(fd, u8View, 0, bufOffset); + } + fs.closeSync(fd); +} + +function saveTRX(filepath, obj, originalFilename) { + let dtype = obj.positions_dtype || "float32"; + let ptsData = obj.pts; + if (ptsData instanceof Float64Array) { + dtype = "float64"; + } + + if (originalFilename.includes("f16") || dtype === "float16") { + dtype = "float16"; + ptsData = float32ToFloat16(obj.pts); + } else if (originalFilename.includes("f64")) { + dtype = "float64"; + ptsData = new Float64Array(obj.pts); + } + + const numStreamlines = obj.offsetPt0.length - 1; + const numPoints = obj.pts.length / 3; + + let header = { + "VOXEL_TO_RASMM": [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1] + ], + "DIMENSIONS": [256, 256, 256], + "NB_STREAMLINES": numStreamlines, + "NB_VERTICES": numPoints + }; + + if (obj.header) { + header.VOXEL_TO_RASMM = obj.header.VOXEL_TO_RASMM || header.VOXEL_TO_RASMM; + header.DIMENSIONS = obj.header.DIMENSIONS || header.DIMENSIONS; + } + + const zipObj = {}; + zipObj["header.json"] = fflate.strToU8(JSON.stringify(header, null, 4)); + zipObj[`positions.3.${dtype}`] = new Uint8Array(ptsData.buffer, ptsData.byteOffset, ptsData.byteLength); + + let offsetDtype = "uint32"; + let offsetData = obj.offsetPt0.subarray(0, numStreamlines + 1); + if (originalFilename.includes("ui64")) { + offsetDtype = "uint64"; + const u64Bytes = new Uint8Array((numStreamlines + 1) * 8); + const view = new DataView(u64Bytes.buffer); + for (let i = 0; i <= numStreamlines; i++) { + view.setUint32(i * 8, offsetData[i], true); + view.setUint32(i * 8 + 4, 0, true); + } + zipObj[`offsets.${offsetDtype}`] = u64Bytes; + } else { + zipObj[`offsets.${offsetDtype}`] = new Uint8Array(offsetData.buffer, offsetData.byteOffset, offsetData.byteLength); + } + + function getDtypeExt(vals) { + if (vals instanceof Float64Array) return 'float64'; + if (vals instanceof Float32Array) return 'float32'; + if (vals instanceof Uint32Array) return 'uint32'; + if (vals instanceof Int32Array) return 'int32'; + if (vals instanceof Uint16Array) return 'uint16'; + if (vals instanceof Int16Array) return 'int16'; + if (vals instanceof Uint8Array) return 'uint8'; + if (vals instanceof Int8Array) return 'int8'; + return 'float32'; + } + + function getCorrectFname(prop) { + let name = prop.fname || prop.id; + let parts = name.split('.'); + let ext = getDtypeExt(prop.vals); + if (parts.length >= 2) { + parts[parts.length - 1] = ext; + return parts.join('.'); + } + return name + '.' + ext; + } + + if (obj.dpv) { + for (let prop of obj.dpv) { + zipObj[`dpv/${getCorrectFname(prop)}`] = new Uint8Array(prop.vals.buffer, prop.vals.byteOffset, prop.vals.byteLength); + } + } + if (obj.dps) { + for (let prop of obj.dps) { + zipObj[`dps/${getCorrectFname(prop)}`] = new Uint8Array(prop.vals.buffer, prop.vals.byteOffset, prop.vals.byteLength); + } + } + if (obj.dpg) { + for (let prop of obj.dpg) { + zipObj[`dpg/${getCorrectFname(prop)}`] = new Uint8Array(prop.vals.buffer, prop.vals.byteOffset, prop.vals.byteLength); + } + } + if (obj.groups) { + for (let prop of obj.groups) { + zipObj[`groups/${getCorrectFname(prop)}`] = new Uint8Array(prop.vals.buffer, prop.vals.byteOffset, prop.vals.byteLength); + } + } + + const zipped = fflate.zipSync(zipObj); + fs.writeFileSync(filepath, zipped); +} From 012d04a2d0c89d28e4552d976c01ff2b237ae4ab Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 22 Jun 2026 15:28:43 -0400 Subject: [PATCH 4/4] ->Implemented robust NIfTI reference header parsing and strict buffer bounds checking --- streamlineIO.mjs | 204 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 156 insertions(+), 48 deletions(-) diff --git a/streamlineIO.mjs b/streamlineIO.mjs index 71c5525..d9ca13c 100644 --- a/streamlineIO.mjs +++ b/streamlineIO.mjs @@ -1,7 +1,7 @@ //Install dependencies // npm install gl-matrix fflate fzstd -export { readTRK, readTCK, readVTK, readTRX, readTT, saveTCK, saveTRK, saveVTK, saveTRX }; +export { readTRK, readTCK, readVTK, readTRX, readTT, saveTCK, saveTRK, saveVTK, saveTRX, readNiftiHeader }; import { mat3, mat4, vec3, vec4 } from "gl-matrix"; //for trk import * as fs from "fs"; import * as fflate from "fflate"; @@ -1069,12 +1069,13 @@ function saveTCK(filepath, obj) { fs.closeSync(fd); } -function saveTRK(filepath, obj, originalFilename) { +function saveTRK(filepath, obj, originalFilename, refHeader = null) { const fd = fs.openSync(filepath, 'w'); const headerBytes = new Uint8Array(1000); const view = new DataView(headerBytes.buffer); - headerBytes.set([84, 82, 65, 67, 75], 0); + headerBytes.set([84, 82, 65, 67, 75], 0); // 'TRACK' + headerBytes.set([82, 65, 83, 0], 948); // 'RAS\0' voxel_order let dim = [256, 256, 256]; let voxelSize = [1, 1, 1]; @@ -1085,17 +1086,16 @@ function saveTRK(filepath, obj, originalFilename) { [0, 0, 0, 1] ]; - if (obj.header) { - if (obj.header.DIMENSIONS) dim = obj.header.DIMENSIONS; - if (obj.header.VOXEL_TO_RASMM) { - voxToRas = obj.header.VOXEL_TO_RASMM; - voxelSize = [ - Math.sqrt(voxToRas[0][0]**2 + voxToRas[1][0]**2 + voxToRas[2][0]**2), - Math.sqrt(voxToRas[0][1]**2 + voxToRas[1][1]**2 + voxToRas[2][1]**2), - Math.sqrt(voxToRas[0][2]**2 + voxToRas[1][2]**2 + voxToRas[2][2]**2) - ]; - } - } + const resolvedHeader = (obj.header && obj.header.VOXEL_TO_RASMM) ? obj.header : refHeader; + if (!resolvedHeader) throw new Error("TCK/VTK → TRK requires a reference NIfTI header (pass refHeader)"); + + dim = resolvedHeader.DIMENSIONS; + voxToRas = resolvedHeader.VOXEL_TO_RASMM; + voxelSize = [ + Math.sqrt(voxToRas[0][0]**2 + voxToRas[1][0]**2 + voxToRas[2][0]**2), + Math.sqrt(voxToRas[0][1]**2 + voxToRas[1][1]**2 + voxToRas[2][1]**2), + Math.sqrt(voxToRas[0][2]**2 + voxToRas[1][2]**2 + voxToRas[2][2]**2) + ]; view.setInt16(6, dim[0], true); view.setInt16(8, dim[1], true); @@ -1112,20 +1112,35 @@ function saveTRK(filepath, obj, originalFilename) { off += 4; } } + // Removed hardcoded RAS voxel order to prevent coordinate flipping by Nibabel + // Reconstruct the exact readTRK scaling matrix + let zoomMat = mat4.fromValues( + 1 / voxelSize[0], 0, 0, -0.5, + 0, 1 / voxelSize[1], 0, -0.5, + 0, 0, 1 / voxelSize[2], -0.5, + 0, 0, 0, 1 + ); - headerBytes.set([82, 65, 83, 0], 948); - - let invMat = null; - if (obj.header && obj.header.VOXEL_TO_RASMM) { - let mat = mat4.fromValues( - voxToRas[0][0], voxToRas[1][0], voxToRas[2][0], voxToRas[3][0], - voxToRas[0][1], voxToRas[1][1], voxToRas[2][1], voxToRas[3][1], - voxToRas[0][2], voxToRas[1][2], voxToRas[2][2], voxToRas[3][2], - voxToRas[0][3], voxToRas[1][3], voxToRas[2][3], voxToRas[3][3] - ); - invMat = mat4.create(); - mat4.invert(invMat, mat); - } + // Map TRK header to a row-major array matching readTRK behavior + let mat = mat4.create(); + mat[0] = voxToRas[0][0]; mat[1] = voxToRas[0][1]; mat[2] = voxToRas[0][2]; mat[3] = voxToRas[0][3]; + mat[4] = voxToRas[1][0]; mat[5] = voxToRas[1][1]; mat[6] = voxToRas[1][2]; mat[7] = voxToRas[1][3]; + mat[8] = voxToRas[2][0]; mat[9] = voxToRas[2][1]; mat[10] = voxToRas[2][2]; mat[11] = voxToRas[2][3]; + mat[12] = voxToRas[3][0]; mat[13] = voxToRas[3][1]; mat[14] = voxToRas[3][2]; mat[15] = voxToRas[3][3]; + + let vox2mmMat = mat4.create(); + mat4.mul(vox2mmMat, zoomMat, mat); + + // Transpose to column-major for gl-matrix inversion + let colMajorVox2mm = mat4.create(); + mat4.transpose(colMajorVox2mm, vox2mmMat); + + let colMajorMm2trk = mat4.create(); + mat4.invert(colMajorMm2trk, colMajorVox2mm); + + // Transpose back to row-major for straightforward vector application + let mm2trkMat = mat4.create(); + mat4.transpose(mm2trkMat, colMajorMm2trk); const numStreamlines = obj.offsetPt0.length - 1; view.setInt32(988, numStreamlines, true); @@ -1134,7 +1149,7 @@ function saveTRK(filepath, obj, originalFilename) { fs.writeSync(fd, headerBytes); - const chunkSize = 16 * 1024 * 1024; + const chunkSize = 16 * 1024; const buf = new ArrayBuffer(chunkSize); const payloadView = new DataView(buf); const u8View = new Uint8Array(buf); @@ -1164,21 +1179,10 @@ function saveTRK(filepath, obj, originalFilename) { let x = pts[j*3]; let y = pts[j*3 + 1]; let z = pts[j*3 + 2]; - let vx = x, vy = y, vz = z; - - if (invMat) { - let cx = x * invMat[0] + y * invMat[4] + z * invMat[8] + invMat[12]; - let cy = x * invMat[1] + y * invMat[5] + z * invMat[9] + invMat[13]; - let cz = x * invMat[2] + y * invMat[6] + z * invMat[10] + invMat[14]; - - vx = (cx + 0.5) * voxelSize[0]; - vy = (cy + 0.5) * voxelSize[1]; - vz = (cz + 0.5) * voxelSize[2]; - } else { - vx = x + 0.5; - vy = y + 0.5; - vz = z + 0.5; - } + + let vx = x * mm2trkMat[0] + y * mm2trkMat[1] + z * mm2trkMat[2] + mm2trkMat[3]; + let vy = x * mm2trkMat[4] + y * mm2trkMat[5] + z * mm2trkMat[6] + mm2trkMat[7]; + let vz = x * mm2trkMat[8] + y * mm2trkMat[9] + z * mm2trkMat[10] + mm2trkMat[11]; payloadView.setFloat32(bufOffset, vx, true); payloadView.setFloat32(bufOffset + 4, vy, true); @@ -1256,7 +1260,7 @@ function saveVTK(filepath, obj) { fs.closeSync(fd); } -function saveTRX(filepath, obj, originalFilename) { +function saveTRX(filepath, obj, originalFilename, refHeader = null) { let dtype = obj.positions_dtype || "float32"; let ptsData = obj.pts; if (ptsData instanceof Float64Array) { @@ -1286,10 +1290,10 @@ function saveTRX(filepath, obj, originalFilename) { "NB_VERTICES": numPoints }; - if (obj.header) { - header.VOXEL_TO_RASMM = obj.header.VOXEL_TO_RASMM || header.VOXEL_TO_RASMM; - header.DIMENSIONS = obj.header.DIMENSIONS || header.DIMENSIONS; - } + const resolvedHeader = (obj.header && obj.header.VOXEL_TO_RASMM) ? obj.header : refHeader; + if (!resolvedHeader) throw new Error("TCK/VTK → TRX requires a reference NIfTI header (pass refHeader)"); + header.VOXEL_TO_RASMM = resolvedHeader.VOXEL_TO_RASMM; + header.DIMENSIONS = resolvedHeader.DIMENSIONS; const zipObj = {}; zipObj["header.json"] = fflate.strToU8(JSON.stringify(header, null, 4)); @@ -1353,7 +1357,111 @@ function saveTRX(filepath, obj, originalFilename) { zipObj[`groups/${getCorrectFname(prop)}`] = new Uint8Array(prop.vals.buffer, prop.vals.byteOffset, prop.vals.byteLength); } } + // Optionally save extra header files + if (obj.header && obj.header.extraFiles) { + for (const [fname, content] of Object.entries(obj.header.extraFiles)) { + zipObj[fname] = content; + } + } const zipped = fflate.zipSync(zipObj); fs.writeFileSync(filepath, zipped); } + +function readNiftiHeader(niftiPath) { + if (!fs.existsSync(niftiPath)) { + throw new Error("NIfTI file not found: " + niftiPath); + } + const buffer = fs.readFileSync(niftiPath); + if (buffer.byteLength < 4) { + throw new Error("Invalid NIfTI file: too small to contain header size"); + } + const view = new DataView(buffer.buffer, buffer.byteOffset, buffer.byteLength); + const sizeof_hdr = view.getInt32(0, true); + + let isNifti1 = false; + let isNifti2 = false; + if (sizeof_hdr === 348) { + if (buffer.byteLength < 348) throw new Error("Invalid NIfTI file: file smaller than expected NIfTI-1 header"); + isNifti1 = true; + } else if (sizeof_hdr === 540) { + if (buffer.byteLength < 540) throw new Error("Invalid NIfTI file: file smaller than expected NIfTI-2 header"); + isNifti2 = true; + } else { + throw new Error("Invalid NIfTI file: unknown sizeof_hdr " + sizeof_hdr); + } + + let dim = []; + let pixdim = []; + let qform_code, sform_code; + let quatern_b, quatern_c, quatern_d; + let qoffset_x, qoffset_y, qoffset_z; + let srow_x = [], srow_y = [], srow_z = []; + + if (isNifti1) { + for(let i=0; i<8; i++) dim.push(view.getInt16(40 + i*2, true)); + for(let i=0; i<8; i++) pixdim.push(view.getFloat32(76 + i*4, true)); + qform_code = view.getInt16(252, true); + sform_code = view.getInt16(254, true); + quatern_b = view.getFloat32(256, true); + quatern_c = view.getFloat32(260, true); + quatern_d = view.getFloat32(264, true); + qoffset_x = view.getFloat32(268, true); + qoffset_y = view.getFloat32(272, true); + qoffset_z = view.getFloat32(276, true); + for(let i=0; i<4; i++) srow_x.push(view.getFloat32(280 + i*4, true)); + for(let i=0; i<4; i++) srow_y.push(view.getFloat32(296 + i*4, true)); + for(let i=0; i<4; i++) srow_z.push(view.getFloat32(312 + i*4, true)); + } else { + for(let i=0; i<8; i++) dim.push(Number(view.getBigInt64(16 + i*8, true))); + for(let i=0; i<8; i++) pixdim.push(view.getFloat64(80 + i*8, true)); + qform_code = view.getInt32(344, true); + sform_code = view.getInt32(348, true); + quatern_b = view.getFloat64(352, true); + quatern_c = view.getFloat64(360, true); + quatern_d = view.getFloat64(368, true); + qoffset_x = view.getFloat64(376, true); + qoffset_y = view.getFloat64(384, true); + qoffset_z = view.getFloat64(392, true); + for(let i=0; i<4; i++) srow_x.push(view.getFloat64(400 + i*8, true)); + for(let i=0; i<4; i++) srow_y.push(view.getFloat64(432 + i*8, true)); + for(let i=0; i<4; i++) srow_z.push(view.getFloat64(464 + i*8, true)); + } + + let DIMENSIONS = [dim[1], dim[2], dim[3]]; + let VOXEL_TO_RASMM = []; + + if (sform_code > 0) { + VOXEL_TO_RASMM = [srow_x, srow_y, srow_z, [0, 0, 0, 1]]; + } else if (qform_code > 0) { + let b = quatern_b; + let c = quatern_c; + let d = quatern_d; + let a = Math.sqrt(Math.max(0, 1.0 - (b*b + c*c + d*d))); + let qfac = pixdim[0] === 0 ? 1 : pixdim[0]; + let dx = pixdim[1]; + let dy = pixdim[2]; + let dz = pixdim[3]; + + let R00 = a*a + b*b - c*c - d*d; + let R01 = 2*(b*c - a*d); + let R02 = 2*(b*d + a*c); + let R10 = 2*(b*c + a*d); + let R11 = a*a + c*c - b*b - d*d; + let R12 = 2*(c*d - a*b); + let R20 = 2*(b*d - a*c); + let R21 = 2*(c*d + a*b); + let R22 = a*a + d*d - c*c - b*b; + + VOXEL_TO_RASMM = [ + [R00 * dx, R01 * dy, R02 * qfac * dz, qoffset_x], + [R10 * dx, R11 * dy, R12 * qfac * dz, qoffset_y], + [R20 * dx, R21 * dy, R22 * qfac * dz, qoffset_z], + [0, 0, 0, 1] + ]; + } else { + throw new Error("NIfTI file has no valid spatial transform"); + } + + return { DIMENSIONS, VOXEL_TO_RASMM }; +}