Image Preprocessing And Tiling

Image Preprocessing and Tiling

Failure modes here are silent, not exceptions: a -9999 nodata sentinel shifts a channel mean by orders of magnitude, a mask gets bilinear-interpolated into non-integer class ids, augmentation randomness is duplicated across eight workers, inference normalization differs from training by a factor of 255. The model still trains. It just does not work.

1. Tiling: choosing tile size and stride

Why overlap at inference: a pixel 3 px from a tile border saw almost no context on one side, so its prediction is systematically worse. With stride == tile you keep exactly those bad predictions along every seam - that is the visible grid pattern. Overlap lets you replace edge predictions with center predictions from the neighbour. Naive overwrite is not enough: writing one tile over another just moves the discontinuity to the write boundary. You need a smooth weight.

2. Reassembly with feathered (Hann) blending

import numpy as np

def hann2d(th, tw, eps=1e-6):
    # np.hanning(N) is EXACTLY 0 at both endpoints. If a border pixel is covered
    # by only one tile, its accumulated weight would be 0 -> 0/0 -> NaN.
    # Sampling the interior of a length-(N+2) window avoids the zeros entirely.
    wy = np.hanning(th + 2)[1:-1]
    wx = np.hanning(tw + 2)[1:-1]
    return np.maximum(np.outer(wy, wx), eps).astype(np.float32)

def tile_starts(size, tile, stride):
    """Start indices covering [0, size). Final tile is SHIFTED to end exactly at
    `size` instead of padding - so no tile ever contains synthetic pixels."""
    if size <= tile:
        return [0]
    starts = list(range(0, size - tile + 1, stride))
    if starts[-1] != size - tile:
        starts.append(size - tile)
    return starts

def predict_scene(image, model_fn, tile=512, stride=256, n_classes=1):
    """image: (C, H, W) float32, already normalized. model_fn: (C,t,t)->(K,t,t)."""
    C, H, W = image.shape
    if H < tile or W < tile:                         # scene smaller than one tile:
        pad = ((0, 0), (0, max(0, tile - H)), (0, max(0, tile - W)))
        big = np.pad(image, pad, mode="symmetric")   # pad up, predict, crop back
        return predict_scene(big, model_fn, tile, stride, n_classes)[:, :H, :W]
    acc  = np.zeros((n_classes, H, W), np.float32)   # weighted sum of logits/probs
    wacc = np.zeros((1, H, W), np.float32)           # accumulated weight
    win  = hann2d(tile, tile)[None]                  # (1,t,t) broadcasts over K

    for y in tile_starts(H, tile, stride):
        for x in tile_starts(W, tile, stride):
            patch = image[:, y:y+tile, x:x+tile]
            pred  = model_fn(patch)                  # (K, t, t)
            acc[:,  y:y+tile, x:x+tile] += pred * win
            wacc[:, y:y+tile, x:x+tile] += win
    return acc / wacc

Non-obvious points in that code:

• Divide by accumulated weight, not by a count. The shifted final tile overlaps its neighbour by an irregular amount; some pixels get 2 tiles, some 4, some 3. Only acc / wacc is correct for all of them, and it makes the irregular last-row/last-column stride a non-issue.
• Blend probabilities or logits, not argmax labels. Averaging class indices is meaningless. Prefer averaging softmax probabilities (bounded, well-behaved); averaging logits is also fine and slightly sharper, but do not mix the two across tiles.
• acc must be float32 minimum. float16 accumulation over 4 overlapping tiles of a 20k x 20k scene loses precision visibly in low-probability regions.
• Memory. acc for a 20000x20000 scene with 10 classes in float32 is 16 GB. Use np.lib.format.open_memmap, or process the scene in row-blocks.
• The torch equivalent is F.unfold/F.fold: F.fold sums overlaps, so fold pred * win and separately fold win broadcast to the same shape, then divide. But unfold materializes every tile at once and will OOM on a large scene.
• 3D volumes: do not hand-roll this. monai.inferers.sliding_window_inference(inputs, roi_size, sw_batch_size, predictor, overlap=0.5, mode="gaussian") is the same algorithm with a Gaussian instead of a Hann window, batched patches, and device=/sw_device= so the accumulator can live on CPU while the model runs on GPU. mode="constant" is the seam-producing default trap.
• Round-trip assertion (section 10): with model_fn = lambda p: p, the reassembled output must equal the input to ~1e-5. If it does not, your seams are a tiling bug, not a model bug.

3. Edges and remainders

Two valid strategies for H % tile != 0:

If you pad, pad by reflection, not zeros:

pad_h, pad_w = (-H) % tile, (-W) % tile
padded = np.pad(image, ((0,0), (0,pad_h), (0,pad_w)), mode="reflect")
# then crop back: out = out[:, :H, :W]

• Zero padding creates a hard step from real radiance to exactly 0. Convolutions near the border see an edge that exists in every padded sample at the same place, so the network learns a border-specific response - a bright or dark 5-20 px rim on predictions and worse metrics near scene edges. 0 is also a legitimate reflectance value and a common nodata sentinel.
• np.pad(mode="reflect") tolerates pad_width >= dim (it reflects repeatedly), but the same request raises in torch.nn.functional.pad(mode="reflect") and in cv2.copyMakeBorder, both of which require pad < dim (pad < dim - 1 for BORDER_REFLECT_101). Padding a 3-px strip by 100 therefore works offline in numpy and blows up inside the model or an albumentations transform. In numpy, mode="symmetric" repeats the edge sample under no length constraint - the robust default. Either way, padding a strip by more than its own width fabricates most of the tile; drop such slivers instead of padding them.
• In OpenCV/albumentations the equivalent is cv2.BORDER_REFLECT_101 (= numpy reflect) vs cv2.BORDER_REFLECT (= numpy symmetric). BORDER_CONSTANT is the zero-padding trap.
• Pad the mask with an ignore label, not with class 0. Reflect-padding the image but zero-padding the mask teaches the model that reflected content is background. Use fill_mask=255 (or your ignore_index) and set ignore_index=255 in the loss.

4. Normalization

The four rules

1. Statistics come from the training split only. Train+val+test is leakage - small, but exactly what reviewers ask about, and not small at all with a distribution-shifted test set.
2. Per channel. A 12-band Sentinel-2 stack has band means spanning an order of magnitude; one scalar mean/std flattens the informative bands into noise.
3. ImageNet statistics are only valid for 8-bit RGB natural photographs. For multispectral, SAR (dB, can be negative), thermal, DEM (metres), or 16-bit medical data they are meaningless. Even for 3-band RGB satellite imagery they are wrong - the radiometry is different.
4. Save the statistics next to the checkpoint and reload them at inference. Recomputing them from the inference scene is the most common source of train/inference skew, and it presents exactly as a model that "doesn't generalize".

Computing stats without poisoning them

import json, numpy as np, rasterio

C = 12
n = np.zeros(C, np.float64); s = np.zeros(C, np.float64); ss = np.zeros(C, np.float64)
SHIFT = None  # provisional mean, prevents catastrophic cancellation for large values

for path in TRAIN_PATHS:                        # TRAIN SPLIT ONLY
    with rasterio.open(path) as src:
        a = src.read(masked=True)               # masked=True honours the nodata tag
    a = a.astype(np.float64)
    valid = ~np.ma.getmaskarray(a)
    # Also drop known sentinels the file failed to declare:
    valid &= np.isfinite(a) & (a != -9999) & (a != 0)   # tune per dataset!
    if SHIFT is None:
        SHIFT = np.array([a[c][valid[c]].mean() if valid[c].any() else 0.0
                          for c in range(C)])
    d = np.where(valid, a - SHIFT[:, None, None], 0.0)
    n  += valid.sum(axis=(1, 2), dtype=np.float64)
    s  += d.sum(axis=(1, 2), dtype=np.float64)
    ss += (d * d).sum(axis=(1, 2), dtype=np.float64)

assert (n > 0).all(), n          # a channel with 0 valid pixels -> NaN stats, silently
mean = SHIFT + s / n
std  = np.sqrt(np.maximum(ss / n - (s / n) ** 2, 0.0))
json.dump({"mean": mean.tolist(), "std": std.tolist(), "n": n.tolist()},
          open("norm_stats.json", "w"))

Gotchas encoded above:

• src.read() returns nodata pixels as their raw sentinel value. masked=True gives a MaskedArray using the file's declared nodata - but many GeoTIFFs declare nothing, or declare 0 while also using -9999 in some tiles. Print src.nodata, then histogram one scene and look for a spike at 0 / -9999 / 65535 before trusting it. A -9999 sentinel in 1% of pixels drags a band with mean 1500 / std 300 to mean 1385 and std 1182 - the mean moves 8%, the std nearly 4x. Every downstream z-score is then compressed toward zero.
• Always pass dtype=np.float64 to .sum(). Integer reductions accumulate in the platform's default integer, which is 32-bit on Windows under numpy < 2.0 - uint16 sums overflow silently past 4.3e9 (a single 20k x 20k band of mid-range reflectance). NumPy >= 2.0 made the Windows default 64-bit, so this bug appears and disappears with the numpy version; pin the accumulator.
• ss/n - mean^2 is numerically fragile for uint16 reflectance (values ~1e4, N ~1e10). The provisional-mean shift above costs nothing and removes the problem; np.maximum(..., 0) guards the residual negative variance.
• Do not compute stats over the tiles your foreground-biased sampler emits - that biases the statistics toward the rare class. Compute over the full training scenes.

Percentile clipping

For imagery with heavy tails (specular water glint, clouds, sensor spikes, saturated pathology whites), z-scoring alone leaves 99.9% of pixels squeezed into a narrow band.

# Sample ~1e7 valid pixels per channel from the TRAIN split, then:
p1, p99 = np.percentile(sampled, [1, 99], axis=0)      # shape (C,)
x = np.clip(x, p1[:, None, None], p99[:, None, None])
x = (x - p1[:, None, None]) / (p99 - p1)[:, None, None]   # -> [0, 1]

Order matters: clip first, then scale. Compute the percentiles globally over the training set, never per-scene at inference - per-scene percentiles make the same physical radiance map to different network inputs depending on what else is in the scene, destroying cross-scene comparability. (Per-scene is acceptable only for visualization, or when sensor gain genuinely varies per acquisition and you have decided to normalize it away.) The right percentile is data dependent: 1/99 is fine on cloud-free agricultural scenes; 2/98 is often needed with water glint or cloud.

Domain notes

• CT is calibrated in Hounsfield units. Window first (soft tissue ~ center 40 / width 400; lung ~ -600 / 1500; bone ~ 400 / 1800), then scale to [0,1]. Global stats are meaningful.
• MRI intensities are in arbitrary units that vary by scanner, coil and sequence. Global mean/std across a dataset is wrong. Normalize per volume (z-score over the brain mask, or scale by the 99th percentile inside the mask), and do it after skull-stripping / bias-field correction if you are doing those.
• SAR in linear power is log-normal - convert to dB (10*log10) before z-scoring, and handle zeros before the log.
• Voxel spacing before patching (CT/MRI). A 96x96x96 patch is a different physical volume in every scan until you resample to a fixed spacing (e.g. 1x1x1 mm, or the dataset median as nnU-Net does). CT series are routinely 0.7x0.7x5 mm - patching that anisotropic grid directly means the z-axis receptive field covers 10x the anatomy the in-plane one does. Resample the image with linear/spline interpolation and the label map with nearest neighbour, and re-check the label set afterwards. Also apply the DICOM RescaleSlope/RescaleIntercept before treating values as HU; pydicom does not apply them for you (apply_modality_lut does).
• Adapting a pretrained RGB encoder to N bands. Do not feed 12 bands through an ImageNet-stat normalization. Compute your own N-channel stats, then inflate the stem weights: w_new = w.mean(1, keepdim=True).repeat(1, N, 1, 1) * (3 / N) - the 3/N keeps the stem's output magnitude where the pretrained downstream layers expect it. Without the rescale, a 12-band stem emits ~4x the activation scale and the first epochs are spent undoing it.

5. Augmentation with albumentations

import cv2, numpy as np, albumentations as A
from albumentations.pytorch import ToTensorV2

MEAN, STD = stats["mean"], stats["std"]   # loaded from norm_stats.json

train_tf = A.Compose([
    # --- geometric: applied identically to image AND mask ---
    A.PadIfNeeded(min_height=512, min_width=512,
                  border_mode=cv2.BORDER_REFLECT_101),
    A.RandomCrop(height=512, width=512),
    A.HorizontalFlip(p=0.5),
    A.VerticalFlip(p=0.5),        # overhead imagery ONLY - see section 6
    A.RandomRotate90(p=0.5),      # overhead imagery ONLY
    A.Affine(scale=(0.9, 1.1), translate_percent=(-0.05, 0.05),
             rotate=(-15, 15), p=0.5),
    # --- photometric: image only, albumentations handles this automatically ---
    A.RandomBrightnessContrast(brightness_limit=0.2, contrast_limit=0.2, p=0.3),
    # --- normalization + tensor ---
    A.Normalize(mean=MEAN, std=STD, max_pixel_value=1.0),
    ToTensorV2(),
])

out = train_tf(image=img_hwc, mask=mask_hw)   # img (H,W,C), mask (H,W) int
x, y = out["image"], out["mask"]               # (C,H,W) float32, (H,W) int64-able

The traps, in order of how often they bite:

• A.Normalize(max_pixel_value=255.0) is the default. If your array is already float in [0,1], or is 16-bit reflectance you scaled yourself, the default divides by 255 again and your inputs land near 0.004. The model trains, slowly, to a mediocre score. Pass max_pixel_value=1.0 whenever the input is not raw uint8.

• ToTensorV2 does NOT divide by 255 (unlike torchvision.transforms.ToTensor); it only transposes HWC->CHW. Porting a torchvision pipeline without adding an explicit scale is a silent 255x error.

• ToTensorV2 does not transpose the mask unless you pass transpose_mask=True. A (H,W,K) one-hot mask stays (H,W,K) and your loss silently broadcasts wrong.

• Masks are interpolated with nearest neighbour by albumentations' geometric transforms, so class ids stay integral - but only if you pass the array as mask=. If you sneak the mask through as a 4th image channel to "save a call", it gets bilinear-interpolated and you get class id 3.7. Assert set(np.unique(y)) <= set(range(K)) | {ignore_index} after augmentation.

• additional_targets is how you get a second mask, a second image, or a change-detection pair transformed with the same random parameters:

  tf = A.Compose([...], additional_targets={"image_t2": "image", "mask_prev": "mask"})
  out = tf(image=a, image_t2=b, mask=m, mask_prev=mp)
  

  Without it, calling the transform twice draws two independent random states and image and mask are no longer aligned - which presents as a model that cannot learn, not as a bug.

• Photometric augs must not touch the mask. Albumentations enforces this for its built-ins (apply_to_mask is identity for brightness/contrast/blur/noise). For a custom transform, subclass A.ImageOnlyTransform, not A.DualTransform.

• Version rename (albumentations 2.0): the border/fill kwargs were inconsistent before 2.0 and were unified in it. Pre-2.0: most geometric transforms (PadIfNeeded, Rotate, ShiftScaleRotate) took border_mode= / value= / mask_value=, while A.Affine took mode= / cval= / cval_mask=. From 2.0 everything is border_mode= / fill= / fill_mask=. Print A.__version__ and inspect.signature(A.Affine.__init__) rather than guessing; an unknown kwarg is silently ignored in some versions and raises in others, so a wrong name can mean you were zero-padding all along.

• >3 channels: geometric transforms and Normalize handle arbitrary channel counts. Colour-space transforms (HueSaturationValue, RGBShift, CLAHE, ToGray, most RandomFog-style ones) assume 3 or 1 channels and will raise or silently mangle a 12-band stack. Verify on one sample. Leave is_check_shapes=True (default) on - it catches image/mask shape mismatch.

6. Augmentations that are wrong for your domain

7. Class imbalance in dense prediction

Uniformly sampling tile origins from a scene where the target covers 0.5% of pixels gives you tiles that are ~85-95% completely empty. The gradient is dominated by background, the loss drops fast, and the model converges to predicting all-background - which scores 99.5% pixel accuracy and 0 IoU.

Precompute a tile index with foreground fraction, then sample it.

# Offline, once - tile_starts() is from section 2:
index = []
for scene_id, mask in enumerate(train_masks):            # each (H, W) uint8
    H, W = mask.shape
    for y in tile_starts(H, TILE, TILE):
        for x in tile_starts(W, TILE, TILE):
            fg = float((mask[y:y+TILE, x:x+TILE] > 0).mean())
            index.append({"scene": scene_id, "y": y, "x": x, "fg": fg})

# In the Dataset: 50/50 mix
fg_pool = [i for i, r in enumerate(index) if r["fg"] > 0.01]
bg_pool = [i for i, r in enumerate(index) if r["fg"] <= 0.01]

def __getitem__(self, i):
    pool = fg_pool if (i % 2 == 0) else bg_pool
    rec = index[random.choice(pool)]   # global RNG: PyTorch reseeds it per worker.
    ...                                # An RNG built in __init__ would not be - section 9.

A.CropNonEmptyMaskIfExists(height, width, p=1.0) crops around a non-zero mask region when one exists and falls back to a random crop otherwise - a fast baseline, but it centres crops on foreground, biasing object position toward the tile centre. Mix it at p=0.5 with plain RandomCrop to avoid that.

• Oversampling foreground changes the class prior, so predicted probabilities are miscalibrated relative to the real scene. Tune any decision threshold on a validation set sampled the way deployment samples (full scenes), not on the oversampled tiles.
• Do not stack mechanisms: aggressive oversampling plus class-weighted loss plus focal loss overshoots into over-prediction. Pick one primary (usually sampling) and one mild secondary (Dice).
• Keep genuinely-empty tiles in the mix. A model that never sees pure background produces false positives everywhere on the real scene.

8. DataLoader performance

Decoding is almost always the bottleneck - not the GPU, not the disk. A 512x512 JPEG decode is ~2-5 ms; a compressed-tile GeoTIFF window read 10-50 ms; an OpenSlide read_region at level 0 50-200 ms. 8 workers x 3 ms is ~2600 img/s (more than most models consume); 8 workers x 50 ms is 160 img/s and your GPU sits idle.

loader = DataLoader(
    ds, batch_size=16, shuffle=True,
    num_workers=8,              # start at physical cores; measure, do not guess
    pin_memory=True,            # only helps if you then use .to(dev, non_blocking=True)
    persistent_workers=True,    # essential on Windows: avoids respawning every epoch
    prefetch_factor=4,          # batches queued per worker; raise if timings are spiky
    drop_last=True,
)

Storage format decision table

Rule of thumb: convert to memmap or WebDataset when a profiling run shows GPU utilization below ~70% and worker count is already at core count.

Whole-slide images (OpenSlide)

• Windows install. Python 3.8+ ignores PATH when resolving extension DLLs, so import openslide fails with a DLL load error until you point at the unpacked binaries first:

  import os
  with os.add_dll_directory(r"C:\openslide\bin"):   # before the import, every process
      import openslide
  

  Because workers spawn, this must run at module import, not once in main().

• read_region((x, y), level, (w, h)) takes (x, y) in the level-0 frame no matter which level you pass, while (w, h) is in that level's pixels. Scaling the origin by slide.level_downsamples[level] "to be consistent" is the classic bug: you read a region downsample^2 away from where you meant, and every tile is misaligned with its annotation.

• It returns RGBA. .convert("RGB") explicitly - the alpha channel is 0 outside the scanned area and naive np.array(region)[..., :3] leaves those regions black, which your tissue filter then happily accepts as "dark = tissue".

• Filter glass before tiling. 70-90% of a slide is background. Otsu-threshold the saturation channel of a low-resolution level (slide.get_thumbnail(...) or level_count - 1), then keep only tiles whose tissue fraction exceeds ~0.1. Doing this offline into a tile index (section 7) is the difference between a 3-hour and a 30-hour epoch.

• Normalize magnification, not pixels. slide.properties[openslide.PROPERTY_NAME_MPP_X] is 0.25 um/px at 40x and 0.5 at 20x, and it is missing on some NDPI/scanner exports. A fixed 512-px tile therefore covers 2x different tissue across a multi-site cohort - resample to a target MPP and fail loudly when the property is absent rather than assuming 40x.

Windows specifics

• Windows uses spawn, not fork. num_workers > 0 requires your training entry point to be guarded, or each worker re-executes the script and you get an infinite spawn storm (usually presenting as a RuntimeError about the current process finishing bootstrapping, or as the machine simply freezing):

  if __name__ == "__main__":
      main()
  

• Everything crossing the process boundary must be picklable: no lambdas in worker_init_fn, no local closures in collate_fn, no partial over a nested function. Module-level functions only.

• Spawn re-imports your module in every worker, so module-level heavy work (loading a big index, importing torch, opening a DB) is paid num_workers times per spawn. Worker startup is ~1-3 s each on Windows vs ~50 ms on Linux, so num_workers=16 can cost 30 s per epoch in pure startup; 4-8 is the sweet spot even on a 16-core machine, and persistent_workers=True (pay once per run, not once per epoch) is often the single biggest wall-clock win.

• Open file handles lazily, inside the worker. h5py.File, rasterio.open, lmdb.open and OpenSlide handles created in Dataset.__init__ are either unpicklable (spawn crashes with a confusing pickling error) or shared unsafely. Standard pattern:

  def __getitem__(self, i):
      if self._h5 is None:              # set to None in __init__
          self._h5 = h5py.File(self.path, "r")
      ...
  

• Add cv2.setNumThreads(0) at module level (and OMP_NUM_THREADS=1). OpenCV's thread pool times 8 worker processes oversubscribes the CPU and can make the loader slower than num_workers=0.

9. Determinism in the data path

import random, numpy as np, torch

def seed_worker(worker_id):                  # must be module-level (Windows pickling)
    s = torch.initial_seed() % 2**32         # per-worker, derived from base_seed
    np.random.seed(s)
    random.seed(s)
    tf = torch.utils.data.get_worker_info().dataset.transform
    if hasattr(tf, "set_random_seed"):       # albumentations >= 2.0, see below
        tf.set_random_seed(s)

g = torch.Generator()
g.manual_seed(1337)                          # controls the shuffle order
loader = DataLoader(ds, ..., worker_init_fn=seed_worker, generator=g)

• The generator= argument seeds the sampler (which indices, in what order). worker_init_fn seeds the augmentation RNGs. You need both; they are independent.
• Modern PyTorch does seed random and numpy's global RNG per worker, but a np.random.RandomState or random.Random instance you construct in Dataset.__init__ is copied to every worker identically - so all 8 workers draw the same augmentation sequence. Construct per-worker RNGs lazily inside __getitem__/worker_init_fn, or use the global RNG.
• Albumentations changed RNG model in 2.0 and seed_worker no longer covers it. Pre-2.0 transforms drew from the global random / numpy RNGs, so seeding the globals per worker was enough. From 2.0 (and late 1.4.x) every transform owns a per-instance RNG fixed at construction time, plus A.Compose(..., seed=N). A Compose built in the parent process and shipped to workers therefore carries the same RNG state into all of them - eight workers, one augmentation sequence - and nothing in worker_init_fn that touches globals will change that. Either reseed the pipeline per worker (the set_random_seed lines above) or construct the Compose lazily inside the worker. Symptom if you miss it: the k-th sample produced by every worker gets identical augmentation parameters, so each batch contains num_workers copies of the same flip/rotate/brightness draw applied to different tiles. Verify by transforming one fixed array in each worker and comparing hashes - they must differ.
• To reproduce one exact sample, use A.ReplayCompose and store the returned replay dict.
• Reproducing a run needs the sampler seed, the worker seed scheme, num_workers, the dataset ordering, and library versions - log all of them. Changing num_workers changes which sample gets which seed, so a run is not reproducible across a worker-count change even with identical seeds.
• torch.use_deterministic_algorithms(True) plus CUBLAS_WORKSPACE_CONFIG=:4096:8 handles the model side; it does nothing for the data path.

10. Debugging: how to actually verify the pipeline

Do these before training anything. Each one has caught a real, silent, model-killing bug.

1. Dump a batch and look at it with your eyes.

import torchvision
x, y = next(iter(train_loader))
vis = x[:, :3]                                     # pick RGB bands for multispectral
vis = (vis - vis.amin()) / (vis.amax() - vis.amin() + 1e-8)
torchvision.utils.save_image(torchvision.utils.make_grid(vis, nrow=4), "batch.png")
torchvision.utils.save_image(
    torchvision.utils.make_grid((y.float() / max(1, y.max())).unsqueeze(1), nrow=4),
    "batch_mask.png")

Then overlay (vis[:, 0] = torch.where(y > 0, 1.0, vis[:, 0]), save again). A side-by-side will not reveal a 1-pixel misalignment or a transposed mask; an overlay will.

2. Assert dtype and range after every stage. Put these in the Dataset temporarily:

assert img.dtype == np.float32, img.dtype
assert np.isfinite(img).all()
assert -6 < img.mean() < 6 and 0.2 < img.std() < 5, (img.mean(), img.std())
assert mask.dtype in (np.uint8, np.int64), mask.dtype
assert set(np.unique(mask)).issubset(ALLOWED_IDS), np.unique(mask)

The mask-class assertion catches: bilinear-interpolated masks, an unremapped label file where classes are 0/38/75/113 instead of 0/1/2/3, and a padding fill of 0 colliding with a real class.

3. Round-trip the tiler with an identity model.

img = np.random.rand(3, 1731, 2049).astype(np.float32)   # deliberately not divisible
rec = predict_scene(img, model_fn=lambda p: p, tile=512, stride=256, n_classes=3)
assert np.allclose(rec, img, atol=1e-5), np.abs(rec - img).max()

Use non-square, indivisible dimensions - a 2048x2048 test image passes with almost any broken implementation.

4. Verify inference normalization equals training normalization. Store the stats inside the checkpoint and have the inference script read them from there, so they cannot drift; if they live in a sidecar file, assert equality of the loaded dicts at inference startup.

5. Overfit 4 samples to ~0 loss. With augmentation off. If the model cannot memorize 4 tiles, the bug is in the data path (misaligned mask, wrong loss target dtype, ignore_index eating everything), not in the architecture or LR.

6. Count label pixels across the whole training set, print the per-class fraction. A class at 0.0000 is absent - either genuinely, or because your remap dropped it.

7. Flag constant tiles. img.std() < 1e-6 means an all-nodata tile reached training; those go NaN under per-tile normalization and poison the whole batch.