From 2c8fcafcb1f5a9c24392c5c7de1be558c08cfc8d Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 27 Sep 2026 12:49:12 +0000 Subject: [PATCH 1/3] Support 64-bit integer datatypes (int64/uint64) NIfTI datatypes 1024/1280, MIF Int64/UInt64 and NRRD 64-bit integer types were rejected, which broke numpy-generated label maps (int64 is NumPy's default integer dtype). Add them to MiqDatatype, decode them through double -> float like MIQCore, make them label-detection candidates, and give ScanVolume0's binary-mask confirm a 64-bit arm. Port of MIQ 2cf1764. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LYbaAXasvPKzEgLcoESA9s --- MIQ.Core/Parsing/MifParser.cs | 2 + MIQ.Core/Parsing/MiqDatatype.cs | 6 +- MIQ.Core/Parsing/NrrdParser.cs | 7 +- MIQ.Core/Rendering/MiqVolume.cs | 11 ++- MIQ.Tests/Int64DatatypeTests.cs | 170 ++++++++++++++++++++++++++++++++ 5 files changed, 190 insertions(+), 6 deletions(-) create mode 100644 MIQ.Tests/Int64DatatypeTests.cs diff --git a/MIQ.Core/Parsing/MifParser.cs b/MIQ.Core/Parsing/MifParser.cs index 836ca04..1884fde 100644 --- a/MIQ.Core/Parsing/MifParser.cs +++ b/MIQ.Core/Parsing/MifParser.cs @@ -408,6 +408,8 @@ private static (MiqDatatype datatype, bool littleEndian, bool isBit) ParseDataty "int32" or "int32_t" => MiqDatatype.Int32, "float32" => MiqDatatype.Float32, "float64" => MiqDatatype.Float64, + "uint64" or "uint64_t" => MiqDatatype.Uint64, + "int64" or "int64_t" => MiqDatatype.Int64, _ => throw new MiqException($"MIF: unsupported datatype '{s}'."), }; diff --git a/MIQ.Core/Parsing/MiqDatatype.cs b/MIQ.Core/Parsing/MiqDatatype.cs index 88eb8f6..22695ee 100644 --- a/MIQ.Core/Parsing/MiqDatatype.cs +++ b/MIQ.Core/Parsing/MiqDatatype.cs @@ -14,6 +14,8 @@ public enum MiqDatatype : short Uint32 = 768, Float32 = 16, Float64 = 64, + Int64 = 1024, + Uint64 = 1280, Rgb24 = 128, Rgba32 = 2304, } @@ -27,7 +29,7 @@ public static class MiqDatatypeExtensions MiqDatatype.Int8 or MiqDatatype.Uint8 => 1, MiqDatatype.Int16 or MiqDatatype.Uint16 => 2, MiqDatatype.Int32 or MiqDatatype.Uint32 or MiqDatatype.Float32 => 4, - MiqDatatype.Float64 => 8, + MiqDatatype.Float64 or MiqDatatype.Int64 or MiqDatatype.Uint64 => 8, MiqDatatype.Rgb24 => 3, MiqDatatype.Rgba32 => 4, _ => throw new ArgumentOutOfRangeException(nameof(dt)), @@ -43,6 +45,8 @@ public static class MiqDatatypeExtensions MiqDatatype.Uint32 => "uint32", MiqDatatype.Float32 => "float32", MiqDatatype.Float64 => "float64", + MiqDatatype.Int64 => "int64", + MiqDatatype.Uint64 => "uint64", MiqDatatype.Rgb24 => "rgb24", MiqDatatype.Rgba32 => "rgba32", _ => throw new ArgumentOutOfRangeException(nameof(dt)), diff --git a/MIQ.Core/Parsing/NrrdParser.cs b/MIQ.Core/Parsing/NrrdParser.cs index f840d4c..26d49c5 100644 --- a/MIQ.Core/Parsing/NrrdParser.cs +++ b/MIQ.Core/Parsing/NrrdParser.cs @@ -388,10 +388,9 @@ private static MiqDatatype ParseDatatype(string value) => "float" => MiqDatatype.Float32, "double" => MiqDatatype.Float64, "int64" or "int64_t" or "longlong" or "long long" - or "signed long long" or "signed long long int" - or "uint64" or "uint64_t" or "ulonglong" - or "unsigned long long" or "unsigned long long int" => throw new MiqException( - "NRRD 64-bit integer types are not supported; convert to float or int32 before previewing."), + or "signed long long" or "signed long long int" => MiqDatatype.Int64, + "uint64" or "uint64_t" or "ulonglong" + or "unsigned long long" or "unsigned long long int" => MiqDatatype.Uint64, _ => throw new MiqException($"Unrecognised NRRD type '{value}'."), }; diff --git a/MIQ.Core/Rendering/MiqVolume.cs b/MIQ.Core/Rendering/MiqVolume.cs index d3e7763..52e0d7f 100644 --- a/MIQ.Core/Rendering/MiqVolume.cs +++ b/MIQ.Core/Rendering/MiqVolume.cs @@ -627,6 +627,11 @@ static long Rd64(byte[] a, long o, bool le) for (var o = p; o < end; o += 4) { int v = Rd32(s, o, le); if (v != 0 && v != label) return Vol0LabelShape.MultiLabel; } break; + case MiqDatatype.Int64: + case MiqDatatype.Uint64: // a uint64 above long.MaxValue reads negative: still MultiLabel + for (var o = p; o < end; o += 8) + { long v = Rd64(s, o, le); if (v != 0 && v != label) return Vol0LabelShape.MultiLabel; } + break; case MiqDatatype.Float32: for (var o = p; o < end; o += 4) { @@ -674,7 +679,7 @@ private Vol0LabelShape ScanVolume0PerVoxel(int label) // Floats included: label maps are often re-saved as float by downstream // tools. Continuous float intensity fails the integrality check. MiqDatatype.Int8 or MiqDatatype.Uint8 or MiqDatatype.Int16 or MiqDatatype.Uint16 - or MiqDatatype.Int32 or MiqDatatype.Uint32 + or MiqDatatype.Int32 or MiqDatatype.Uint32 or MiqDatatype.Int64 or MiqDatatype.Uint64 or MiqDatatype.Float32 or MiqDatatype.Float64 => true, _ => false, }; @@ -929,6 +934,10 @@ private float RawVoxelValue(int byteOffset) MiqDatatype.Uint32 => MiqBinaryReader.Uint32(s, abs, le), MiqDatatype.Float32 => MiqCompat.Int32BitsToSingle((int)MiqBinaryReader.Uint32(s, abs, le)), MiqDatatype.Float64 => (float)MiqCompat.Int64BitsToDouble((long)MiqBinaryReader.Uint64(s, abs, le)), + // Via double to match MIQCore's Float(Double(x)); a direct cast rounds + // once and can differ from it in the last bit. + MiqDatatype.Int64 => (float)(double)MiqBinaryReader.Int64(s, abs, le), + MiqDatatype.Uint64 => (float)(double)MiqBinaryReader.Uint64(s, abs, le), // RGB normally takes ReadRgb; luminance fallback just in case. MiqDatatype.Rgb24 or MiqDatatype.Rgba32 => 0.299f * _image.Byte(byteOffset) + 0.587f * _image.Byte(byteOffset + 1) diff --git a/MIQ.Tests/Int64DatatypeTests.cs b/MIQ.Tests/Int64DatatypeTests.cs new file mode 100644 index 0000000..867ca43 --- /dev/null +++ b/MIQ.Tests/Int64DatatypeTests.cs @@ -0,0 +1,170 @@ +using System.Buffers.Binary; +using System.Text; +using MIQ.Parsing; +using MIQ.Rendering; +using Xunit; + +namespace MIQ.Tests; + +// 64-bit integer voxels (NIfTI 1024/1280, MIF Int64/UInt64, NRRD int64/uint64) +// used to be rejected, which broke numpy-generated label maps (int64 is NumPy's +// default integer dtype). Port of MIQ@2cf1764. +public class Int64DatatypeTests +{ + const int N1MinOffset = 352; + + /// NIfTI-1, little-endian, 8 bytes per voxel; voxel i holds valueAt(i). + static byte[] Nifti1(MiqDatatype datatype, int w, int h, int d, Func valueAt) + { + var n = w * h * d; + var file = new byte[N1MinOffset + n * 8]; + var hdr = file.AsSpan(); + BinaryPrimitives.WriteInt32LittleEndian(hdr.Slice(0), 348); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(40), 3); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(42), (short)w); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(44), (short)h); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(46), (short)d); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(70), (short)datatype); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(72), 64); + for (var i = 0; i < 4; i++) + BinaryPrimitives.WriteSingleLittleEndian(hdr.Slice(76 + i * 4), 1f); + BinaryPrimitives.WriteSingleLittleEndian(hdr.Slice(108), N1MinOffset); + for (var i = 0; i < n; i++) + BinaryPrimitives.WriteUInt64LittleEndian(hdr.Slice(N1MinOffset + i * 8), valueAt(i)); + return file; + } + + /// Same volume as int16, the reference the 64-bit label maps must match. + static byte[] Nifti1Int16(int w, int h, int d, Func valueAt) + { + var n = w * h * d; + var file = new byte[N1MinOffset + n * 2]; + var hdr = file.AsSpan(); + BinaryPrimitives.WriteInt32LittleEndian(hdr.Slice(0), 348); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(40), 3); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(42), (short)w); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(44), (short)h); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(46), (short)d); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(70), (short)MiqDatatype.Int16); + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(72), 16); + for (var i = 0; i < 4; i++) + BinaryPrimitives.WriteSingleLittleEndian(hdr.Slice(76 + i * 4), 1f); + BinaryPrimitives.WriteSingleLittleEndian(hdr.Slice(108), N1MinOffset); + for (var i = 0; i < n; i++) + BinaryPrimitives.WriteInt16LittleEndian(hdr.Slice(N1MinOffset + i * 2), valueAt(i)); + return file; + } + + // ── Decode ────────────────────────────────────────────────────────────── + + [Fact] + public void Nifti_Int64_DecodesBeyondInt32Range() + { + // -5e9 is below int32.MinValue and exact in float, so a 32-bit read can't pass. + const long v = -5_000_000_000L; + var image = NiftiParser.Parse(Nifti1(MiqDatatype.Int64, 2, 2, 2, _ => unchecked((ulong)v))); + Assert.Equal(MiqDatatype.Int64, image.Header.Datatype); + Assert.Equal("int64", image.Header.Datatype.Label()); + Assert.Equal(-5e9f, new MiqVolume(image).VoxelValue(1, 1, 1, 0)); + } + + [Fact] + public void Nifti_Uint64_DecodesAboveInt64Max() + { + // 2^63 reads as long.MinValue if the unsigned read is missed. + const ulong v = 1UL << 63; + var image = NiftiParser.Parse(Nifti1(MiqDatatype.Uint64, 2, 2, 2, _ => v)); + Assert.Equal(MiqDatatype.Uint64, image.Header.Datatype); + Assert.Equal(9.223372E18f, new MiqVolume(image).VoxelValue(1, 1, 1, 0)); + } + + [Theory] + [InlineData("Int64LE", MiqDatatype.Int64)] + [InlineData("UInt64BE", MiqDatatype.Uint64)] + public void Mif_64BitTypes_Parse(string mifType, MiqDatatype expected) + { + var le = mifType.EndsWith("LE"); + var header = Encoding.ASCII.GetBytes( + $"mrtrix image\ndim: 2,2,2\nvox: 1,1,1\nlayout: +0,+1,+2\ndatatype: {mifType}\nfile: .\nEND\n"); + var payload = new byte[8 * 8]; + for (var i = 0; i < 8; i++) + { + var span = payload.AsSpan(i * 8); + if (le) BinaryPrimitives.WriteInt64LittleEndian(span, 6_000_000_000L + i); + else BinaryPrimitives.WriteInt64BigEndian(span, 6_000_000_000L + i); + } + var image = MifParser.Parse([.. header, .. payload]); + Assert.Equal(expected, image.Header.Datatype); + Assert.Equal(6e9f, new MiqVolume(image).VoxelValue(0, 0, 0, 0)); + } + + [Theory] + [InlineData("int64", MiqDatatype.Int64)] + [InlineData("signed long long int", MiqDatatype.Int64)] + [InlineData("uint64_t", MiqDatatype.Uint64)] + [InlineData("unsigned long long", MiqDatatype.Uint64)] + public void Nrrd_64BitTypes_Parse(string nrrdType, MiqDatatype expected) + { + var header = Encoding.ASCII.GetBytes( + $"NRRD0004\ntype: {nrrdType}\ndimension: 3\nsizes: 2 2 2\nencoding: raw\nendian: little\n\n"); + var payload = new byte[8 * 8]; + for (var i = 0; i < 8; i++) + BinaryPrimitives.WriteInt64LittleEndian(payload.AsSpan(i * 8), 42 + i); + var image = NrrdParser.Parse([.. header, .. payload]); + Assert.Equal(expected, image.Header.Datatype); + Assert.Equal(49f, new MiqVolume(image).VoxelValue(1, 1, 1, 0)); + } + + // ── Label detection ───────────────────────────────────────────────────── + + static readonly MiqRenderingOptions Auto = new(Segmentation: MiqSegmentationColoring.Auto); + + // FreeSurfer-style labels in contiguous 4³ blocks (piecewise constant). + static readonly int[] LabelsFS = [0, 2, 3, 41, 42, 10, 11, 17, 251]; + const int N = 16; + + static int BlockLabel(int i) + { + int x = i % N, y = i / N % N, z = i / (N * N); + return LabelsFS[(x / 4 + y / 4 * 4 + z / 4 * 16) % LabelsFS.Length]; + } + + [Theory] + [InlineData(MiqDatatype.Int64)] + [InlineData(MiqDatatype.Uint64)] + public void LabelMap_DetectedSameAsInt16Equivalent(MiqDatatype datatype) + { + var reference = new MiqVolume(NiftiParser.Parse(Nifti1Int16(N, N, N, i => (short)BlockLabel(i)))) + .BuildSegmentationLut(Auto); + Assert.NotNull(reference); + + var lut = new MiqVolume(NiftiParser.Parse(Nifti1(datatype, N, N, N, i => (ulong)BlockLabel(i)))) + .BuildSegmentationLut(Auto); + Assert.NotNull(lut); + Assert.Equal(reference!.IsFreeSurfer, lut!.IsFreeSurfer); + Assert.Equal(reference.IsMonochromeWhite, lut.IsMonochromeWhite); + Assert.True(lut.IsFreeSurfer); + } + + static bool InBox(int i) + { + int x = i % N, y = i / N % N, z = i / (N * N); + return x is >= 4 and < 12 && y is >= 4 and < 12 && z is >= 4 and < 12; + } + + [Fact] + public void Int64BinaryMask_IsConfirmedByFullScan() + { + // The center slices see one label, so ScanVolume0's 64-bit arm decides. + var mask = new MiqVolume(NiftiParser.Parse( + Nifti1(MiqDatatype.Int64, N, N, N, i => InBox(i) ? 5UL : 0UL))); + Assert.True(mask.BuildSegmentationLut(Auto)!.IsMonochromeWhite); + + // A second label at (0,0,0), off every center plane, must be found by it. + var twoLabels = new MiqVolume(NiftiParser.Parse( + Nifti1(MiqDatatype.Int64, N, N, N, i => i == 0 ? 7UL : InBox(i) ? 5UL : 0UL))); + var lut = twoLabels.BuildSegmentationLut(Auto); + Assert.NotNull(lut); + Assert.False(lut!.IsMonochromeWhite); + } +} From 93c4fc351f68097b1e12c76516e0a82186b88cb7 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 27 Sep 2026 12:49:12 +0000 Subject: [PATCH 2/3] Window tiny-unit float data over its non-zero voxels IntensityWindow.GetBounds preferred the |v| > 1e-6 subset, so data in tiny units (an SI-unit ADC map, ~1e-9) had no foreground subset and was windowed over tissue plus background zeros. When the floored subset is too small, exactly-non-zero values now stand in. Ordinary data is unchanged (pinned bit-identical against the previous implementation). Port of MIQ bcf3e18. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LYbaAXasvPKzEgLcoESA9s --- MIQ.Core/Rendering/IntensityWindow.cs | 41 +++++++-- MIQ.Tests/TinyUnitWindowTests.cs | 115 ++++++++++++++++++++++++++ 2 files changed, 151 insertions(+), 5 deletions(-) create mode 100644 MIQ.Tests/TinyUnitWindowTests.cs diff --git a/MIQ.Core/Rendering/IntensityWindow.cs b/MIQ.Core/Rendering/IntensityWindow.cs index 97a49e3..0a7b918 100644 --- a/MIQ.Core/Rendering/IntensityWindow.cs +++ b/MIQ.Core/Rendering/IntensityWindow.cs @@ -10,6 +10,10 @@ public readonly struct Bounds(float low, float high) public float High { get; } = high; } + /// Magnitude at or below which a value counts as background for the preferred + /// foreground subset (see ). + internal const float NonZeroFloor = 1e-6f; + /// Derives window bounds from a pooled value set. Returns null if no finite /// values are present. public static Bounds? GetBounds(IReadOnlyList values, double lowerPercentile, double upperPercentile) @@ -19,12 +23,39 @@ public readonly struct Bounds(float low, float high) if (MiqCompat.IsFinite(v)) finite.Add(v); if (finite.Count == 0) return null; - // Prefer a non-zero subset if substantial; the /20 ratio guards against - // rejecting legitimate dim regions when most voxels are background. - var nonZero = new List(finite.Count); + // Window over the foreground when it's substantial; the /20 ratio keeps a + // dim region from being rejected when most voxels are background. + // + // Two tiers, tried in order. The NonZeroFloor also drops near-zero + // interpolation residue, and wins whenever it leaves enough voxels (every + // ordinary image, unchanged from before). Only when it doesn't (data in + // tiny units: an SI-unit ADC map sits around 1e-9, entirely below the + // floor) do exactly-non-zero values stand in, so such a map windows over + // its tissue instead of over tissue plus background zeros. Deliberately + // not a floor scaled by the data's maximum: one outlier voxel (1e30 from a + // failed fit) would lift that above real tissue. Port of MIQCore's bounds. + var minimumSubset = Math.Max(64, finite.Count / 20); + var aboveFloor = new List(finite.Count); + var nonZeroCount = 0; foreach (var v in finite) - if (Math.Abs(v) > 1e-6f) nonZero.Add(v); - var source = nonZero.Count >= Math.Max(64, finite.Count / 20) ? nonZero : finite; + { + if (v == 0f) continue; + nonZeroCount++; + if (Math.Abs(v) > NonZeroFloor) aboveFloor.Add(v); + } + + List source; + if (aboveFloor.Count >= minimumSubset) + source = aboveFloor; + else if (nonZeroCount >= minimumSubset) + { + // Tiny-unit fallback: the only case that pays a second pass. + source = new List(nonZeroCount); + foreach (var v in finite) + if (v != 0f) source.Add(v); + } + else + source = finite; source.Sort(); var lower = Percentile(source, (float)lowerPercentile / 100f); diff --git a/MIQ.Tests/TinyUnitWindowTests.cs b/MIQ.Tests/TinyUnitWindowTests.cs new file mode 100644 index 0000000..a82962c --- /dev/null +++ b/MIQ.Tests/TinyUnitWindowTests.cs @@ -0,0 +1,115 @@ +using MIQ.Rendering; +using Xunit; + +namespace MIQ.Tests; + +// IntensityWindow.GetBounds preferred a foreground subset of |v| > 1e-6. Data in +// tiny units (an SI-unit ADC map, ~1e-9) sits entirely below that floor, so its +// window was computed over tissue plus the background zeros. Exactly-non-zero +// values now stand in when the floored subset is too small. Port of MIQ@bcf3e18. +public class TinyUnitWindowTests +{ + /// The previous GetBounds, verbatim: the reference for data the change must + /// leave alone. + static IntensityWindow.Bounds? LegacyBounds(IReadOnlyList values, double lo, double hi) + { + var finite = values.Where(float.IsFinite).ToList(); + if (finite.Count == 0) return null; + var nonZero = finite.Where(v => Math.Abs(v) > 1e-6f).ToList(); + var source = nonZero.Count >= Math.Max(64, finite.Count / 20) ? nonZero : finite; + source.Sort(); + float Pct(float p) + { + if (source.Count == 1) return source[0]; + var pos = Math.Max(0f, Math.Min(1f, p)) * (source.Count - 1); + int l = (int)Math.Floor((double)pos), u = (int)Math.Ceiling((double)pos); + if (l == u) return source[l]; + var f = pos - l; + return source[l] * (1 - f) + source[u] * f; + } + var lower = Pct((float)lo / 100f); + var upper = Pct((float)hi / 100f); + return lower < upper + ? new IntensityWindow.Bounds(lower, upper) + : new IntensityWindow.Bounds(finite.Min(), finite.Max()); + } + + /// A center slice in tiny units: 60% exact-zero background, tissue spread over + /// 0.5e-9 … 3e-9. Every value sits below the 1e-6 floor. + static float[] TinyUnitSlice(float? outlier = null) + { + var values = new float[10_000]; + for (var i = 0; i < values.Length; i++) + values[i] = i % 5 < 3 ? 0f : 0.5e-9f + (i % 997) / 996f * 2.5e-9f; + if (outlier is { } o) values[1] = o; + return values; + } + + [Fact] + public void TinyUnitData_WindowsOverItsTissue() + { + var b = IntensityWindow.GetBounds(TinyUnitSlice(), 2, 98); + Assert.NotNull(b); + // Lower bound comes from the tissue, not from the background zeros. + Assert.True(b!.Value.Low > 0.5e-9f, $"low {b.Value.Low}"); + Assert.True(b.Value.High < 3e-9f, $"high {b.Value.High}"); + + // Before, the zeros pulled the lower bound down to 0. + Assert.Equal(0f, LegacyBounds(TinyUnitSlice(), 2, 98)!.Value.Low); + } + + [Fact] + public void TinyUnitWindow_IgnoresAnOutlier() + { + // A failed-fit voxel must not drag the window: the fallback is a plain + // non-zero test, not a floor scaled by the data's maximum. + var b = IntensityWindow.GetBounds(TinyUnitSlice(outlier: 1e30f), 2, 98); + Assert.NotNull(b); + Assert.True(b!.Value.Low > 0.5e-9f, $"low {b.Value.Low}"); + Assert.True(b.Value.High < 3e-9f, $"high {b.Value.High}"); + } + + [Fact] + public void Floor_StillDropsNearZeroResidue_OnOrdinaryData() + { + // Ordinary-unit signal plus interpolation residue in (0, 1e-6]: the floor + // tier qualifies, so the residue stays out of the window, as before. + var values = new float[10_000]; + for (var i = 0; i < values.Length; i++) + values[i] = i % 2 == 0 ? (i % 7 + 1) * 1e-8f : 100 + i % 500; + var b = IntensityWindow.GetBounds(values, 2, 98); + Assert.NotNull(b); + Assert.True(b!.Value.Low >= 100, $"low {b.Value.Low}"); + AssertSame(LegacyBounds(values, 2, 98), b); + } + + public static TheoryData OrdinaryShapes => new() { "anatomical", "sparse", "all-zero", "few-voxels" }; + + [Theory] + [MemberData(nameof(OrdinaryShapes))] + public void OrdinaryData_IsBitIdenticalToBefore(string shape) + { + var rng = new Random(7); + float[] values = shape switch + { + // Mostly foreground in ordinary units. + "anatomical" => Enumerable.Range(0, 20_000) + .Select(i => i % 4 == 0 ? 0f : (float)(rng.NextDouble() * 1500)).ToArray(), + // Foreground below the /20 ratio: both before and after, all finite values. + "sparse" => Enumerable.Range(0, 20_000) + .Select(i => i % 50 == 0 ? (float)(rng.NextDouble() * 10) : 0f).ToArray(), + "all-zero" => new float[4096], + // Fewer than 64 non-zero voxels. + _ => Enumerable.Range(0, 500).Select(i => i < 30 ? 3e-9f * (i + 1) : 0f).ToArray(), + }; + AssertSame(LegacyBounds(values, 2, 98), IntensityWindow.GetBounds(values, 2, 98)); + } + + static void AssertSame(IntensityWindow.Bounds? expected, IntensityWindow.Bounds? actual) + { + Assert.Equal(expected.HasValue, actual.HasValue); + if (expected is not { } e) return; + Assert.Equal(e.Low, actual!.Value.Low); + Assert.Equal(e.High, actual.Value.High); + } +} From a78ad5254c22633204023497e9daf415d00331d2 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 27 Sep 2026 12:49:12 +0000 Subject: [PATCH 3/3] Scan MIF volume 0 contiguously when its layout is dense ScanVolume0 and CollectVolume0Labels took the per-voxel walk whenever ElementStrides was set, and MifParser always sets it, so every MIF mask paid for it. Key off Volume0IsContiguous instead: volume 0 is the first W*H*D elements in some order, which is all these order-independent scans need. A 256^3 MIF mask's BuildSegmentationLut drops from ~750 ms to ~30 ms; layouts with the volume axis faster than a spatial one still walk. Port of the mask-scan part of MIQ 8323423. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01LYbaAXasvPKzEgLcoESA9s --- MIQ.Core/Rendering/MiqVolume.cs | 42 +++++++++-- MIQ.Tests/BinaryMaskLayoutTests.cs | 108 +++++++++++++++++++++++++++++ 2 files changed, 143 insertions(+), 7 deletions(-) create mode 100644 MIQ.Tests/BinaryMaskLayoutTests.cs diff --git a/MIQ.Core/Rendering/MiqVolume.cs b/MIQ.Core/Rendering/MiqVolume.cs index 52e0d7f..887efa6 100644 --- a/MIQ.Core/Rendering/MiqVolume.cs +++ b/MIQ.Core/Rendering/MiqVolume.cs @@ -423,9 +423,9 @@ bool Add(float v) return !(labels.Add(label) && labels.Count > SegmentationLut.MaxLabels); } - if (_image.ElementStrides is not null) + if (!Volume0IsContiguous) { - // MIF custom strides may interleave volumes: per-voxel walk. + // Volumes interleaved with space: per-voxel walk. for (var z = 0; z < Depth; z++) for (var y = 0; y < Height; y++) for (var x = 0; x < Width; x++) @@ -433,7 +433,7 @@ bool Add(float v) } else { - // Row-major: volume 0 is the first W·H·D payload elements. + // Volume 0 is the first W·H·D payload elements (order is irrelevant here). var bpv = H.Datatype.BytesPerVoxel(); var elems = (int)Math.Min((long)Width * Height * Depth, (long)_image.PayloadCount / bpv); for (var i = 0; i < elems; i++) @@ -579,11 +579,11 @@ private static int GetSingle(HashSet set) // binary mask scans to completion. private Vol0LabelShape ScanVolume0(int label) { - // MIF custom strides may interleave volumes: use the per-voxel walk. - if (_image.ElementStrides is not null) + // A layout that interleaves volumes with space: use the per-voxel walk. + if (!Volume0IsContiguous) return ScanVolume0PerVoxel(label); - // Row-major formats: volume 0 is the first N payload elements. Only which + // Volume 0 is the first N payload elements, in some order. Only which // values occur matters, not where, so scan the raw buffer sequentially with // the datatype switch hoisted out of the loop. var s = _image.Storage; @@ -658,7 +658,35 @@ static long Rd64(byte[] a, long o, bool le) return Vol0LabelShape.Binary; } - // Layout-agnostic fallback for MIF custom strides. + /// Whether volume 0 occupies exactly the first W·H·D payload elements, in any + /// order. The volume-0 label scans don't depend on voxel order, so this (not + /// the absence of custom strides) is what their contiguous path needs. + /// + /// True for row-major layouts and for every dense permuted one: the spatial + /// strides sorted ascending are 1, d₀, d₀·d₁. MifParser sets ElementStrides + /// even for an identity layout, so keying off null sent every MIF mask down + /// the per-voxel walk. The time stride never matters: volume 0 is t = 0. A + /// layout with the volume axis faster than a spatial one fails the check, as + /// it must, since its first N elements mix timepoints. Internal so tests can + /// pin which path a layout takes. Port of MIQCore's volumeZeroIsContiguous. + internal bool Volume0IsContiguous + { + get + { + if (_image.ElementStrides is not { } s) return true; + if (s.Length < 3 || _image.BaseElementIndex != 0) return false; + // Ties (a size-1 axis shares its stride with the next) sort by size so + // the size-1 axis comes first and the check still holds. A mis-ordered + // tie would only cost the fast path, never correctness. + var axes = new[] { (stride: s[0], size: Width), (stride: s[1], size: Height), (stride: s[2], size: Depth) }; + Array.Sort(axes, (a, b) => a.stride != b.stride ? a.stride.CompareTo(b.stride) : a.size.CompareTo(b.size)); + return axes[0].stride == 1 + && axes[1].stride == axes[0].size + && axes[2].stride == (long)axes[0].size * axes[1].size; + } + } + + // Layout-agnostic fallback for layouts that interleave volumes with space. private Vol0LabelShape ScanVolume0PerVoxel(int label) { for (var z = 0; z < Depth; z++) diff --git a/MIQ.Tests/BinaryMaskLayoutTests.cs b/MIQ.Tests/BinaryMaskLayoutTests.cs new file mode 100644 index 0000000..c701bc2 --- /dev/null +++ b/MIQ.Tests/BinaryMaskLayoutTests.cs @@ -0,0 +1,108 @@ +using System.Text; +using MIQ.Parsing; +using MIQ.Rendering; +using Xunit; + +namespace MIQ.Tests; + +// The volume-0 label scans take the contiguous path whenever volume 0 is the first +// W·H·D payload elements in *any* order, not only when ElementStrides is null. +// MifParser always sets strides, so keying off null sent every MIF mask down the +// per-voxel walk (~8× slower on a 256³ mask in MIQCore). These pin that dense MIF +// layouts take the fast path with the same verdict, and that a volume-interleaved +// layout still does not. Port of MIQ@8323423's BinaryMaskLayoutTests. +public class BinaryMaskLayoutTests +{ + static readonly MiqRenderingOptions Auto = new(Segmentation: MiqSegmentationColoring.Auto); + + // Distinct sizes so a stride mix-up can't cancel out. The box covers all three + // center planes (x = 6, y = 5, z = 4); voxel (0,0,0) lies on none of them. + static readonly int[] Dims = [12, 10, 8]; + + static bool InBox(int x, int y, int z) => x is >= 3 and < 9 && y is >= 3 and < 7 && z is >= 2 and < 6; + + /// A uint8 MIF whose voxel (x, y, z, t), in `dim:` order, holds value(x, y, z, t), + /// stored in the given layout (MRtrix storage ranks, all ascending). + static MiqVolume MakeMif(int[] dims, int[] layout, Func value) + { + var n = dims.Length; + var byRank = Enumerable.Range(0, n).OrderBy(a => layout[a]).ToArray(); + var strides = new int[n]; + var stride = 1; + foreach (var axis in byRank) { strides[axis] = stride; stride *= dims[axis]; } + + var volumes = n > 3 ? dims[3] : 1; + var payload = new byte[stride]; + for (var t = 0; t < volumes; t++) + for (var z = 0; z < dims[2]; z++) + for (var y = 0; y < dims[1]; y++) + for (var x = 0; x < dims[0]; x++) + { + var index = x * strides[0] + y * strides[1] + z * strides[2] + (n > 3 ? t * strides[3] : 0); + payload[index] = value(x, y, z, t); + } + + var header = Encoding.ASCII.GetBytes( + "mrtrix image\n" + + $"dim: {string.Join(",", dims)}\n" + + $"vox: {string.Join(",", dims.Select(_ => "1"))}\n" + + $"layout: {string.Join(",", layout.Select(r => "+" + r))}\n" + + "datatype: UInt8\nfile: .\nEND\n"); + return new MiqVolume(MifParser.Parse([.. header, .. payload])); + } + + [Theory] + [InlineData(new[] { 0, 1, 2 })] + [InlineData(new[] { 2, 0, 1 })] + [InlineData(new[] { 1, 2, 0 })] + [InlineData(new[] { 1, 0, 2 })] + public void DenseMifMask_TakesContiguousScan(int[] layout) + { + var mask = MakeMif(Dims, layout, (x, y, z, _) => InBox(x, y, z) ? (byte)5 : (byte)0); + Assert.True(mask.Volume0IsContiguous); + Assert.True(mask.BuildSegmentationLut(Auto)!.IsMonochromeWhite); + + // A second label off the center planes must still be found by the full scan. + var twoLabels = MakeMif(Dims, layout, (x, y, z, _) => + x == 0 && y == 0 && z == 0 ? (byte)7 : InBox(x, y, z) ? (byte)5 : (byte)0); + var lut = twoLabels.BuildSegmentationLut(Auto); + Assert.NotNull(lut); + Assert.False(lut!.IsMonochromeWhite); + } + + [Fact] + public void VolumeSlowest4DMif_ScansOnlyVolumeZero() + { + // Volume 0 is a clean mask; volume 1 carries another label. The contiguous + // scan stops at W·H·D elements, so volume 1 must not leak into the verdict. + var volume = MakeMif([.. Dims, 2], [0, 1, 2, 3], (x, y, z, t) => + !InBox(x, y, z) ? (byte)0 : t == 0 ? (byte)5 : (byte)9); + Assert.True(volume.Volume0IsContiguous); + Assert.True(volume.BuildSegmentationLut(Auto)!.IsMonochromeWhite); + } + + [Fact] + public void VolumeFastest4DMif_FallsBackToPerVoxelWalk() + { + // DWI-style layout: the volume axis varies fastest, so the first W·H·D + // elements interleave both timepoints. A contiguous scan would see volume + // 1's label and misreport a multi-label volume. + var volume = MakeMif([.. Dims, 2], [1, 2, 3, 0], (x, y, z, t) => + !InBox(x, y, z) ? (byte)0 : t == 0 ? (byte)5 : (byte)9); + Assert.False(volume.Volume0IsContiguous); + Assert.True(volume.BuildSegmentationLut(Auto)!.IsMonochromeWhite); + } + + [Fact] + public void RowMajorImage_IsContiguous() + { + var header = new MiqHeader + { + LittleEndian = true, Dimensions = [2, 2, 2], Pixdim = [1f, 1f, 1f, 1f], + Datatype = MiqDatatype.Uint8, VoxOffset = 0, SclSlope = 0f, SclInter = 0f, + QformCode = 0, SformCode = 0, SrowX = [1, 0, 0], SrowY = [0, 1, 0], SrowZ = [0, 0, 1], + }; + var vol = new MiqVolume(new MiqImage { Header = header, Storage = new byte[8], PayloadOffset = 0 }); + Assert.True(vol.Volume0IsContiguous); + } +}