matlab-process-large-images

name: matlab-process-large-images description: "Patterns for using blockedImage to process large images, harness parallel compute for image processing, and write custom adapters. Use when writing code that creates, processes, or visualizes blockedImage objects, when implementing images.blocked.Adapter subclasses, or when a user needs help with large image data. Always use this skill when working with TIFF,GeoTIFF, .svs, .ndpi, .czi or other WSI, satellite imagery or microscopy volume image formats." license: MathWorks BSD-3-Clause metadata: author: MathWorks version: "1.0"

When to Use

Writing code that creates, processes, or visualizes blockedImage objects
Implementing images.blocked.Adapter subclasses for custom file formats
User needs help with large image data that doesn't fit in memory
User wants to leverage parallel computing for image processing tasks
Working with whole-slide images, geospatial rasters, 3D microscopy volumes, or tiled mosaics
Any time a user is working with a TIFF file — load this skill and check the file first with imfinfo. Use blockedImage if the file is large (dimensions > 10,000 pixels in any axis), has multiple channels beyond RGB, or contains multiple IFDs (could be pyramid levels, time series, or Z-stack slices)

When NOT to Use

Image fits comfortably in memory and no parallel processing is needed — use standard imread/imwrite workflows
General image display without blockedImage — use the matlab-display-image skill instead
Deep learning model architecture or training loop design — this skill covers data preparation with blockedImageDatastore, not network design

Large Image Processing with blockedImage

Reference files (read on demand when the topic comes up):

references/visualization.md — imageshow, overlays, synced views, overview+detail, annotations
references/training-and-labeling.md — Image Labeler, labeled blockedImages, training datastores, DL inference
references/3d-volumes.md — 3D blocked images, 3D inference, 3D training data
references/geospatial.md — GeoTIFF, shapefiles, coordinate systems, spatial referencing
references/data-discovery.md — 7-step data discovery workflow, custom adapters, format conversion adapter selection, virtual composite adapters, adapter performance tips
references/specialized-workflows.md — writing block-by-block, coordinate systems, mask-based block selection, ExtraImages for pixel-level masking, in-place resume
references/whole-slide-images.md — importing WSI files via OpenSlide/Bio-Formats, choosing between libraries, display, downstream processing (segmentation, classification)

When to use blockedImage

blockedImage is not only for data that exceeds RAM. Even when the data fits in memory, wrapping it as a blockedImage lets you harness a parallel cluster: each worker can independently read and write different blocks, enabling parallel partial reads and parallel writes into separate file blocks.

Always ask the user what their intent with the data is. Ask about visualization, processing, DL training, archival, and whether they want to leverage parallel compute — even if the data fits in RAM.

When the user points you at a directory or multi-file dataset

ALWAYS read references/data-discovery.md and follow its 7-step workflow before writing any code. The workflow covers: inspecting file layout, parsing naming conventions, understanding user intent, choosing between single blockedImage vs. collection, assessing random access capabilities, checking for gaps/overlaps, and recommending an approach (built-in adapter, format conversion, or custom adapter).

Construction

% Read-only (auto-selects adapter; reads only metadata, fast)
bim = blockedImage("large_image.tif");

% Multiple images at once
[bim1, bim2] = blockedImage(["img1.tif", "img2.tif"]);

% Writable
bim = blockedImage(dest, imageSize, blockSize, initVal, Mode="w");

% Override adapter
bim = blockedImage(source, Adapter=myAdapter);

BlockSize

Defaults to a factor of IOBlockSize (adapter's native read unit), targeting ~1024 or the minimum dimension size.
The object's BlockSize is the default for all downstream processing and visualization.
Override per-operation via the BlockSize name-value argument (e.g., when a DL network requires a fixed input size).
For best performance, keep BlockSize as a multiple of IOBlockSize and as large as practical.

Reading Data

getBlock(bim, blocksub) — read one block by block subscript. Returns [data, blockinfo].
getRegion(bim, pixelStart, pixelEnd) — read an arbitrary pixel region.
gather(bim) — load entire image into memory (only if it fits in RAM). Use Level to select resolution.

Processing with `apply()`

apply(bim, @fcn) calls fcn on each block one at a time. Scales to arbitrarily large data. Supports UseParallel=true for cluster execution.

Key name-value arguments: BlockSize, BorderSize, Level, ExtraImages, PadPartialBlocks, PadMethod, OutputLocation, BatchSize, Resume, BlockLocationSet, DisplayWaitbar.

Block struct (`bs`) fields

The user function receives a block struct bs with:

bs.Data — pixel data for this block (includes border if BorderSize > 0)
bs.Start — 1-by-N pixel subscript of the block's top-left corner (excluding border)
bs.End — 1-by-N pixel subscript of the block's bottom-right corner (excluding border)
bs.BorderSize — 1-by-N border size used
bs.Level — resolution level
bs.BlockSub — block subscript (index into SizeInBlocks)
bs.ImageNumber — index of the source image (when processing an array)

Struct / non-image output

When the output per block is not image data (e.g., detection results, counts), use images.blocked.MATBlocks adapter to store structs:

bResults = apply(bim, @countNuclei, ...
    BlockLocationSet=bls, BorderSize=[r r], ...
    Adapter=images.blocked.MATBlocks, OutputLocation="results_dir");
allResults = gather(bResults);  % load all structs into memory

Limit work to real data

Crop first. crop(bim, cstart, cend) is virtual — no data copy, executes instantly. Use to restrict to a spatial ROI or single channel before processing.
Use masks with selectBlockLocations. Build a mask from a coarse level, then process only blocks overlapping the ROI. See references/specialized-workflows.md for the full pattern including InclusionThreshold and block-size tuning.

Key patterns

Single pass. Perform all pre/post-processing inside fcn to go through data only once.
World coordinates. fcn receives block metadata. Use sub2world() to convert block-local positions to world coordinates before returning detection results.
BorderSize for overlap. Set BorderSize when the processing function needs context from neighboring blocks (e.g., DL network with wider input than output field of view).
Unique detections via centroids. Set BorderSize >= half the max object diameter. In the user function, discard detections whose centroid falls in the border region, then convert remaining centroids to global coordinates using bs.Start:
```
onBorder = centroids(:,1) < bs.BorderSize(1) ...
         | centroids(:,1) > size(bs.Data,1) - 2*bs.BorderSize(1) ...
         | centroids(:,2) < bs.BorderSize(2) ...
         | centroids(:,2) > size(bs.Data,2) - 2*bs.BorderSize(2);
centroids(onBorder,:) = [];
centroids = centroids + bs.Start(1:2);  % local -> global
```
For inconsistent detectors, add a second apply() pass with BorderSize=BlockSize to deduplicate by bounding-box overlap in global coordinates.
Multi-resolution output. Use makeMultiLevel2D(bResult, Scales=[1 0.5 0.3 0.1]) to generate a pyramid from a single-level result. Use "nearest" interpolation for label/mask data. Or return multiple levels from fcn and combine with concatenateLevels.
Small block output. When per-block output is tiny (e.g., a count), use gather() to load all results into memory, or write() to rewrite with a larger block size.
Resume. Set Resume=true to continue a previously interrupted apply() from where it left off. For in-place resume with overwrite, see references/specialized-workflows.md.
Multi-output apply. When the processing function returns N outputs, apply returns N separate blockedImages. Use this to split channels in a single pass:
```
[o1, o2, o3, o4] = apply(bim, @(bs) splitChannels(bs), ...
    BlockSize=[10 2048 2048 4], UseParallel=true, OutputLocation=outDir);
```
ExtraImages for pixel-level masking. Pass a second blockedImage via ExtraImages to access it in the block function alongside the primary image. See references/specialized-workflows.md for the full pattern including resolution mismatch handling.

Pyramidal TIFF in one pass. Generate multi-resolution levels and write as a single pyramidal TIFF using LevelImages:

[b1,b2,b3,b4] = apply(bim, @resizeBlocks);
write(bim, "pyramid.tif", LevelImages=[b1 b2 b3 b4], BlockSize=[2048 2048]);
function [b1,b2,b3,b4] = resizeBlocks(bs)
    b1 = imresize(bs.Data, 0.5, 'nearest');
    b2 = imresize(b1, 0.5, 'nearest');
    b3 = imresize(b2, 0.5, 'nearest');
    b4 = imresize(b3, 0.5, 'nearest');
end

BatchSize. Process multiple blocks at once by setting BatchSize=N. The block function receives bs.Data with an extra trailing dimension of size N and bs.BatchSize indicating the count. Useful when per-block overhead is high.

Iterative development workflow

Processing full-resolution blocked images is slow. Shorten the feedback cycle by prototyping on smaller data first:

Prototype on a coarse level. Use gather(bim, Level=bim.NumLevels) to load the coarsest level into memory. Develop and tune parameters on this small image, then apply to the full image:
```
imCoarse = gather(bim, Level=bim.NumLevels);
thresh = graythresh(im2gray(imCoarse));
bResult = apply(bim, @(bs) ~imbinarize(im2gray(bs.Data), thresh));
```

Prototype on an ROI. Extract a small region with representative features using getRegion, tune parameters, then apply:

roi = getRegion(bim, [900 2400], [1700 3300], Level=1);
thresh = graythresh(im2gray(roi));
bResult = apply(bim, @(bs) ~imbinarize(im2gray(bs.Data), thresh));

Parallel processing

When writing results to disk, check for Parallel Computing Toolbox and use it if available:

outDir = tempname;
usepar = ~isempty(ver("parallel"));
bResult = apply(bim, @fcn, UseParallel=usepar, OutputLocation=outDir);

Always include UseParallel when specifying OutputLocation — set it based on toolbox availability.

Visualization

Choose 2D vs 3D based on the data, not convenience:

3D volumes (Z-stacks, CT, MRI, microscopy stacks): Default to volshow + viewer3d. Use OverlayData for multi-channel overlay, SlicePlanes for interactive browsing, and linkviewers for overview+detail. Read references/3d-volumes.md for full patterns. Do NOT flatten 3D data into 2D slice-by-slice viewers — that discards spatial context and ignores the 3D tooling.
2D images (whole-slide, satellite, mosaics): Use imageshow + viewer2d. Read references/visualization.md for overlays, synced views, overview+detail, and annotations.

If the data has a Z/depth/slice dimension, it is 3D — use the 3D path even if each slice is a separate file.

Format Conversion

Use write(bim, destination) to convert between formats or re-chunk data:

write(bim, "output.tif", BlockSize=[512 512 3]);

For adapter selection guidance when converting (JPEG vs PNG vs H5Blocks), see "Adapter Preference for Format Conversion" in references/data-discovery.md.