Skip to content

Fix black previews of off-center masks and >4 GiB .nii.gz expansion #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marcoduering merged 1 commit into from
Sep 27, 2026
Merged

Fix black previews of off-center masks and >4 GiB .nii.gz expansion #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 28 additions & 26 deletions MIQ.Core/Parsing/MiqParser.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,9 @@ public static MiqImage Parse(string filePath, CancellationToken ct = default)
/// <c>.nii.gz</c> (decompresses only volume 0), and to uncompressed <c>.nii</c>
/// above <see cref="PartialLoadThreshold"/> (or <see cref="MaxArrayBytes"/>, in
/// which case it's permanent — <see cref="MiqImage.ExpansionBlocked"/>).
/// Falls back to a full <see cref="Parse"/> for 3-D / small files, or when
/// ISIZE is unreliable; those return a complete image (IsPartial = false).
/// Falls back to a full <see cref="Parse"/> for 3-D / small files, or when the
/// data is shorter than the header implies; those return a complete image
/// (IsPartial = false).
/// </summary>
public static MiqImage ParsePartial(string filePath)
{
Expand All @@ -91,10 +92,6 @@ public static MiqImage ParsePartial(string filePath)
if (kind != MiqFileKind.NiiGz)
return Parse(filePath);

var isize = ReadGzipIsize(filePath);
if (isize == 0)
return Parse(filePath); // unknown / >4 GB — let full path handle it

// Probe: decompress just enough to read the NIfTI header (max 1024 B;
// NIfTI-2 header is 540 B so 1024 covers both variants with margin).
const int probeSize = 1024;
Expand All @@ -121,12 +118,17 @@ public static MiqImage ParsePartial(string filePath)
+ (long)header.Width * header.Height * header.Depth
* header.Datatype.BytesPerVoxel();

if (budget >= (long)isize || budget > int.MaxValue)
return Parse(filePath); // volume 0 is the whole file
if (budget > int.MaxValue)
return Parse(filePath); // volume 0 alone can't be held — full path reports it

// Full size from the header, NOT the gzip ISIZE footer: ISIZE is the size
// mod 2^32, so a >4 GiB series can report a small ISIZE and be offered an
// expansion that can never fit a byte[].
var blocked = NiftiFileExtent(header) > MaxArrayBytes;

// Too few volumes to pay off (see PartialGzipMinVolumes) — unless the full
// data can't fit a byte[], which must stay on the blocked vol-0 view.
if (header.Volumes < PartialGzipMinVolumes && isize <= (ulong)MaxArrayBytes)
if (header.Volumes < PartialGzipMinVolumes && !blocked)
return Parse(filePath);

byte[] partial;
Expand All @@ -137,8 +139,10 @@ public static MiqImage ParsePartial(string filePath)
}
catch { return Parse(filePath); }

if (partial.Length < header.VoxOffset)
return Parse(filePath); // decompression shorter than expected
// Stream ended inside volume 0: the file is shorter than its header says,
// so it is small enough for the full path, which reports the truncation.
if (partial.Length < budget)
return Parse(filePath);

return new MiqImage
{
Expand All @@ -147,10 +151,22 @@ public static MiqImage ParsePartial(string filePath)
PayloadOffset = header.VoxOffset,
IsPartial = true,
// Full data won't fit a byte[] → never expand; show volume 0 + notice.
ExpansionBlocked = isize > (ulong)MaxArrayBytes,
ExpansionBlocked = blocked,
};
}

/// Uncompressed NIfTI file size implied by the header: vox_offset plus every
/// dim[1..7] (dim[5+] also hold data) × bytes-per-voxel. ParseHeader's
/// ValidateDimensionExtent keeps the product within a long; the sum saturates.
internal static long NiftiFileExtent(MiqHeader header)
{
long payload = header.Datatype.BytesPerVoxel();
foreach (var d in header.Dimensions) payload *= Math.Max(1, d);
return payload > long.MaxValue - header.VoxOffset
? long.MaxValue
: payload + header.VoxOffset;
}

/// Loads only volume 0 of an uncompressed multi-volume NIfTI. Reads the header,
/// then exactly header + volume-0 bytes from disk. Falls back to a full
/// <see cref="Parse"/> for a single-volume file (nothing to defer) or when the
Expand Down Expand Up @@ -214,20 +230,6 @@ private static byte[] Trim(byte[] buf, int length)
return t;
}

private static uint ReadGzipIsize(string filePath)
{
try
{
using var fs = new FileStream(filePath, FileMode.Open, FileAccess.Read, FileShare.Read);
if (fs.Length < 18) return 0;
var footer = new byte[4];
fs.Seek(-4, SeekOrigin.End);
if (fs.Read(footer, 0, 4) != 4) return 0;
return (uint)(footer[0] | (footer[1] << 8) | (footer[2] << 16) | (footer[3] << 24));
}
catch { return 0; }
}

/// Rejects a header whose declared voxel extent (product of every axis ×
/// bytes-per-voxel) can't be represented in a long. C# multiplies unchecked, so
/// a wrapped product could slip past a `data.Length &lt; required` check.
Expand Down
140 changes: 126 additions & 14 deletions MIQ.Core/Rendering/MiqVolume.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,20 +176,12 @@ public IReadOnlyDictionary<SlicePlane, CenterSlice> CenterSlices(
() => prepared[2] = PrepareSlice(planes[2]));

// Label detection reuses the prepared gray arrays — no extra decode.
SegmentationLut? lut = null;
if (SegmentationEligible(options) && CollectLabels(prepared) is { } labels
&& (IsPiecewiseConstant(prepared) ?? true))
lut = FinishSegmentationLut(labels, options);
var lut = DetectSegmentationLut(prepared, options);

// Label volumes map through the LUT, not the intensity window.
IntensityWindow.Bounds? window = null;
if (lut is null)
{
var pooled = new List<float>();
foreach (var p in prepared)
if (p.Gray is { } g) pooled.AddRange(g); // RGB slices bypass pooling
window = IntensityWindow.GetBounds(pooled, options.LowerPercentile, options.UpperPercentile);
}
window = WindowFromSlices(prepared, options, timepoint: 0);

var result = new Dictionary<SlicePlane, CenterSlice>();
for (var i = 0; i < planes.Length; i++)
Expand Down Expand Up @@ -283,11 +275,71 @@ internal readonly struct PreparedSlice(float[]? gray, byte[]? rgb, SliceConfig c
/// does not flicker the brightness. Pass a <paramref name="timepoint"/>
/// other than 0 to compute the window for a specific 4-D volume.
public IntensityWindow.Bounds? SharedWindow(MiqRenderingOptions options, int timepoint = 0)
{
var prepared = new[]
{
PrepareSlice(SlicePlane.Coronal, timepoint: timepoint),
PrepareSlice(SlicePlane.Sagittal, timepoint: timepoint),
PrepareSlice(SlicePlane.Axial, timepoint: timepoint),
};
return WindowFromSlices(prepared, options, timepoint);
}

// Pooled center-slice window. When the center slices hold no finite value, or
// only one value (e.g. a mask whose foreground lies elsewhere), that window is
// null or degenerate and would render every slice black; fall back to a
// window over the whole volume at this timepoint.
private IntensityWindow.Bounds? WindowFromSlices(
PreparedSlice[] prepared, MiqRenderingOptions options, int timepoint)
{
var pooled = new List<float>();
foreach (var plane in new[] { SlicePlane.Coronal, SlicePlane.Sagittal, SlicePlane.Axial })
if (PrepareSlice(plane, timepoint: timepoint).Gray is { } g) pooled.AddRange(g);
return IntensityWindow.GetBounds(pooled, options.LowerPercentile, options.UpperPercentile);
foreach (var p in prepared)
if (p.Gray is { } g) pooled.AddRange(g); // RGB slices bypass pooling
var window = IntensityWindow.GetBounds(pooled, options.LowerPercentile, options.UpperPercentile);
if (window is { } w && w.High > w.Low) return window;
if (pooled.Count == 0) return window; // RGB: no intensity window

float? constant = window is { } d ? d.Low : null;
return VolumeWindow(options, timepoint, constant) ?? window;
}

// Upper bound on voxels VolumeWindow sorts; larger sets are strided evenly.
private const long MaxVolumeWindowSamples = 1_000_000;

// Window over every voxel of one timepoint. With a <paramref name="constant"/>
// (the value the center slices were filled with) only the other values are
// sampled, then the window is widened to include the constant, so a sparse
// mask on a zero background maps to [0, label]. Returns null when nothing
// (else) finite is present. Two passes: count, then an evenly strided sample.
private IntensityWindow.Bounds? VolumeWindow(
MiqRenderingOptions options, int timepoint, float? constant)
{
bool Keep(float v) => MiqCompat.IsFinite(v) && (constant is not { } c || v != c);

long count = 0;
for (var z = 0; z < Depth; z++)
for (var y = 0; y < Height; y++)
for (var x = 0; x < Width; x++)
if (Keep(Voxel(x, y, z, timepoint))) count++;
if (count == 0) return null;

var stride = (count + MaxVolumeWindowSamples - 1) / MaxVolumeWindowSamples;
var samples = new List<float>((int)Math.Min(count, MaxVolumeWindowSamples));
long seen = 0;
for (var z = 0; z < Depth; z++)
for (var y = 0; y < Height; y++)
for (var x = 0; x < Width; x++)
{
var v = Voxel(x, y, z, timepoint);
if (Keep(v) && seen++ % stride == 0) samples.Add(v);
}

if (IntensityWindow.GetBounds(samples, options.LowerPercentile, options.UpperPercentile)
is not { } b)
return null;
return constant is { } k
? new IntensityWindow.Bounds(Math.Min(b.Low, k), Math.Max(b.High, k))
: b;
}

/// Extract a single slice at an arbitrary index using a precomputed window
Expand Down Expand Up @@ -333,11 +385,65 @@ internal CenterSlice FinalizeInteractive(
PrepareSlice(SlicePlane.Sagittal),
PrepareSlice(SlicePlane.Axial),
};
return DetectSegmentationLut(prepared, options);
}

// Shared by BuildSegmentationLut and CenterInteractiveState so the initial
// and post-expansion LUTs are decided identically.
private SegmentationLut? DetectSegmentationLut(PreparedSlice[] prepared, MiqRenderingOptions options)
{
if (!SegmentationEligible(options)) return null;
if (CollectLabels(prepared) is not { } labels) return null;

// Center slices are pure background (a lesion or single-ROI mask often
// misses all three): no evidence either way, so collect from all of
// volume 0 instead of falling back to a window of zeros.
if (labels.Count == 0)
return CollectVolume0Labels() is { Count: > 0 } all
? FinishSegmentationLut(all, options, volume0Scanned: true)
: null;

if (!(IsPiecewiseConstant(prepared) ?? true)) return null;
return FinishSegmentationLut(labels, options);
}

// Distinct foreground labels over ALL of volume 0 (present even on partial
// loads). Null on the first fractional value or once more than MaxLabels
// distinct values are seen, mirroring CollectLabels. Only reached when the
// center slices are all background, so a real label volume is sparse here.
// Eligibility guarantees identity scaling, so raw values are the labels.
private HashSet<int>? CollectVolume0Labels()
{
var labels = new HashSet<int>();
bool Add(float v)
{
if (!MiqCompat.IsFinite(v)) return true;
var label = MiqCompat.RoundToInt(v);
if (Math.Abs(v - label) > 1e-3f) return false; // fractional → intensity
return !(labels.Add(label) && labels.Count > SegmentationLut.MaxLabels);
}

if (_image.ElementStrides is not null)
{
// MIF custom strides may interleave volumes: per-voxel walk.
for (var z = 0; z < Depth; z++)
for (var y = 0; y < Height; y++)
for (var x = 0; x < Width; x++)
if (!Add(Voxel(x, y, z, 0))) return null;
}
else
{
// Row-major: volume 0 is the first W·H·D payload elements.
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++)
if (!Add(RawVoxelValue(i * bpv))) return null;
}

labels.Remove(0);
return labels;
}

// Cheap, no-decode eligibility gate: integer/float identity-scaled data in a
// non-Off mode.
private bool SegmentationEligible(MiqRenderingOptions options)
Expand Down Expand Up @@ -431,10 +537,16 @@ private bool SegmentationEligible(MiqRenderingOptions options)
// single-label sample is confirmed against all of volume 0 before committing
// to the (sticky) monochrome LUT. Multi-label volumes get the FreeSurfer
// palette (Auto only) or the random palette.
private SegmentationLut? FinishSegmentationLut(HashSet<int> labels, MiqRenderingOptions options)
// volume0Scanned: the labels already came from all of volume 0, so a single
// label is a confirmed binary mask and needs no second scan.
private SegmentationLut? FinishSegmentationLut(
HashSet<int> labels, MiqRenderingOptions options, bool volume0Scanned = false)
{
if (labels.Count == 0) return null;

if (labels.Count == 1 && volume0Scanned)
return new SegmentationLut(useFreeSurfer: false, monochromeWhite: true);

if (labels.Count == 1)
{
switch (ScanVolume0(GetSingle(labels)))
Expand Down
114 changes: 114 additions & 0 deletions MIQ.Tests/GzipPartialBlockingTests.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
using System.Buffers.Binary;
using System.IO.Compression;
using MIQ.Parsing;
using Xunit;

namespace MIQ.Tests;

// Lowers the process-wide MiqParser.MaxArrayBytes seam, which MifParser and the
// .nii paths also read, so these tests must not overlap any other class.
[CollectionDefinition(nameof(StaticParserSeams), DisableParallelization = true)]
public class StaticParserSeams { }

// ParsePartial's blocked/expandable decision for .nii.gz. The gzip ISIZE footer
// is the uncompressed size mod 2^32, so a >4 GiB series can report a small
// ISIZE; the decision must come from the header. The wrap is simulated by
// patching the footer and lowering MaxArrayBytes (the managed partial gunzip
// stops long before the footer, so the patched value is never validated).
[Collection(nameof(StaticParserSeams))]
public class GzipPartialBlockingTests
{
const int Side = 8;
const int VolumeBytes = Side * Side * Side * 2; // int16
const int Offset = 352;
const int Volumes = 12; // == PartialGzipMinVolumes default
const int FileBytes = Offset + Volumes * VolumeBytes; // 12640

static byte[] Nifti1(int volumes, int payloadBytes)
{
var file = new byte[Offset + payloadBytes];
var h = file.AsSpan();
BinaryPrimitives.WriteInt32LittleEndian(h, 348);
short[] dims = [4, Side, Side, Side, (short)volumes];
for (var i = 0; i < dims.Length; i++)
BinaryPrimitives.WriteInt16LittleEndian(h[(40 + i * 2)..], dims[i]);
BinaryPrimitives.WriteInt16LittleEndian(h[70..], (short)MiqDatatype.Int16);
BinaryPrimitives.WriteInt16LittleEndian(h[72..], 16);
for (var i = 0; i < 8; i++)
BinaryPrimitives.WriteSingleLittleEndian(h[(76 + i * 4)..], 1f);
BinaryPrimitives.WriteSingleLittleEndian(h[108..], Offset);
return file;
}

static byte[] Gzip(byte[] raw, uint? isizeOverride = null)
{
using var ms = new MemoryStream();
using (var gz = new GZipStream(ms, CompressionLevel.Optimal, leaveOpen: true))
gz.Write(raw);
var bytes = ms.ToArray();
if (isizeOverride is { } isize)
BinaryPrimitives.WriteUInt32LittleEndian(bytes.AsSpan(bytes.Length - 4), isize);
return bytes;
}

static MiqImage ParsePartialWith(byte[] gz, long maxArrayBytes)
{
var path = Path.Combine(Path.GetTempPath(), $"miq-block-{Guid.NewGuid():N}.nii.gz");
File.WriteAllBytes(path, gz);
var saved = MiqParser.MaxArrayBytes;
MiqParser.MaxArrayBytes = maxArrayBytes;
try { return MiqParser.ParsePartial(path); }
finally
{
MiqParser.MaxArrayBytes = saved;
File.Delete(path);
}
}

[Fact]
public void WrappedIsize_HeaderTooLarge_IsBlocked()
{
// Header implies 12640 B > limit 10000, but the footer claims 5000 B
// (> volume 0, < limit): the shape of a >4 GiB series whose ISIZE wrapped.
var gz = Gzip(Nifti1(Volumes, Volumes * VolumeBytes), isizeOverride: 5000);
var image = ParsePartialWith(gz, maxArrayBytes: 10_000);
Assert.True(image.IsPartial);
Assert.True(image.ExpansionBlocked);
}

[Fact]
public void HeaderFits_IsExpandable()
{
var gz = Gzip(Nifti1(Volumes, Volumes * VolumeBytes));
var image = ParsePartialWith(gz, maxArrayBytes: FileBytes);
Assert.True(image.IsPartial);
Assert.False(image.ExpansionBlocked);
}

[Fact]
public void FewVolumes_HeaderTooLarge_StaysOnBlockedVolume0View()
{
// Below PartialGzipMinVolumes a fitting file parses fully; one that can't
// fit must still take the blocked partial path, whatever ISIZE says.
const int few = 2;
var gz = Gzip(Nifti1(few, few * VolumeBytes), isizeOverride: 1000);
var image = ParsePartialWith(gz, maxArrayBytes: Offset + VolumeBytes);
Assert.True(image.IsPartial);
Assert.True(image.ExpansionBlocked);
}

[Fact]
public void DataEndsInsideVolume0_FullParseReportsTruncation()
{
// Header claims 12 volumes but only half of volume 0 is present.
var gz = Gzip(Nifti1(Volumes, VolumeBytes / 2));
Assert.Throws<MiqException>(() => ParsePartialWith(gz, maxArrayBytes: 0x7FFFFFC7));
}

[Fact]
public void NiftiFileExtent_CountsEveryDimension()
{
var header = NiftiParser.ParseHeader(Nifti1(Volumes, 0));
Assert.Equal((long)FileBytes, MiqParser.NiftiFileExtent(header));
}
}
Loading
Loading