Process Large GeoTIFF File as Blocked Image

This example shows how to efficiently apply and visualize image processing operations on a large GeoTIFF file.

The goal of this example is to detect boats on an image of a river. Applying the boat detection algorithm directly to a large file requires a lot of memory and potentially wastes computational resources. The example mitigates these challenges by working with the GeoTIFF file as a blocked image. The blocked image approach can scale to very large files because it loads and processes one block of data at a time. The blocked image approach is also computationally efficient because it can omit blocks of nonmeaningful data, such as regions of land, for the boat detection algorithm.

The steps to apply the boat detection algorithm to a blocked image include:

• Read the GeoTIFF file into a blockedImage object.

• Create a low-resolution version of the image.

• Create a mask of the water regions of interest (ROIs) at the low resolution level.

• Apply the boat detection algorithm, at the high resolution level, to blocks that overlap with the mask.

Create Multiresolution Blocked Image

Create a blockedImage object that reads the GeoTIFF file. This image is of Boston and the surrounding region, including the Charles River and Boston Harbor. The image has one resolution level.

bim = blockedImage("boston.tif");

Create an adapter that reads TIFF files and applies JPEG compression.

adapter = images.blocked.TIFF;
adapter.Compression = Tiff.Compression.JPEG;

Create a two-level blockedImage object by using the makeMultiLevel2D function. Specify the block size as 1024-by-1024 pixels, which is generally appropriate for various image sizes and machine configurations. Specify the scale factors as the two-element vector [1 1/8]. The image at the first resolution level consists of the original large image in bim. The image at the second resolution level is 1/8 the size of the original image in each dimension.

timestamp = datetime("now",Format="yyyy-MM-dd-HH-mm-ss");
outputFile = "boston_multiLevel_" + string(timestamp) + ".tif";
bim2 = makeMultiLevel2D(bim,BlockSize=[1024 1024],Scales=[1 1/8], ...
    OutputLocation=outputFile,Adapter=adapter);

Display the multiresolution blocked image using the bigimageshow function. As you zoom in and out of the displayed image, the bigimageshow function switches between the resolution levels.

bigimageshow(bim2)

Create Mask of Water ROI

You can process large images more quickly by applying the operations to only regions of interest (ROIs). For this example, because the ROIs consist of bodies of water, you want the boat detection algorithm to run on only regions of water and not on regions of land.

For this image, regions of water are darker than regions of land. Therefore, you can identify the water ROIs by applying a threshold to the pixel values. Alternatively, if you have a shapefile containing geological or political boundaries, you can create a mask from the shapefile. For an example, see Crop and Mask Large GeoTIFF File Using Shapefile.

Specify a pixel value threshold of 50. Then, apply the threshold to the image at the coarse resolution level and clean up the resulting mask by using the detectWater helper function, which is defined at the end of the example. Call the detectWater helper function on all blocks in the blocked image by using the apply function. Ensure that the first block is full by specifying a block size of 256-by-256 pixels.

maxWaterPixelValue = 50;
lvl = bim2.NumLevels; 
waterMask = apply(bim2,@(bs)detectWater(bs,maxWaterPixelValue), ...
    BlockSize=[256 256],Level=lvl);

Display the mask in blue over a grayscale version of the original image.

hb = bigimageshow(bim2);
showlabels(hb,waterMask,Colormap=[1 1 1; 0 0 1])

Detect Boats Within ROI

Identify the blocks at the high resolution level that are within the ROI by using the selectBlockLocations function and specifying the Masks name-value argument. To include all blocks with at least one pixel within the ROI, also specify the InclusionThreshold name-value argument as 0.

blockSet = selectBlockLocations(bim,BlockSize=[1024 1024], ...
    Masks=waterMask,InclusionThreshold=0);

Specify the minimum pixel values and areas for the boat detection algorithm. For the pixel values, use a minimum threshold of 150, which corresponds to bright pixels. For the areas, specify a minimum area of 5 pixels.

minBoatPixelValue = 150;
minBoatArea = 5;

Apply the boat detection algorithm by using the detectBoats helper function, which is defined at the end of the example. Call the detectBoats helper function on all blocks of the high resolution image by using the apply function. To limit processing to the identified blocks within the ROI, specify the BlockLocationSet name-value argument. To improve the accuracy of detecting boats that span two blocks, specify the BorderSize name-value argument.

detections = apply(bim, ...
    @(block,waterMask)detectBoats(block,waterMask,minBoatPixelValue,minBoatArea), ...
    BlockLocationSet=blockSet,BorderSize=[10 10],ExtraImages=waterMask);

Load all detections from the resulting blockedImage output into an array.

detections = gather(detections);
centroidsXY = vertcat(detections.centroidsXY);

GeoTIFF files have raster reference information embedded in the metadata. This information is required to convert the coordinates of a detected object from raster coordinates to map coordinates. Get the raster reference object for the image by using the georasterinfo function, then convert the raster coordinates to map coordinates.

ginfo = georasterinfo(bim.Source);
R = ginfo.RasterReference;
[centroidsWorldX,centroidsWorldY] = intrinsicToWorld(R,centroidsXY(:,1),centroidsXY(:,2));
centroidsWorld = [centroidsWorldX centroidsWorldY];

Display Results

Raster Coordinates

Display the results in raster coordinates by using the bigimageshow function. Plot the centroids of the detected boats over the image.

figure
bigimageshow(bim2);
hold on
plot(centroidsXY(:,1),centroidsXY(:,2),"ro",MarkerSize=10,LineWidth=2)

When you zoom in, bigimageshow switches to higher available resolutions.

xlim([1380 1900])
ylim([1150 1480])

Map Coordinates

To display the results in map coordinates, use the low-resolution image from the blocked image. Update the corresponding raster reference object to match the size of this level.

imLowRes = gather(bim2);
RLowRes = R;
RLowRes.RasterSize = size(imLowRes);

Display the image and detections in map coordinates by using the mapshow function.

figure
mapshow(imLowRes,RLowRes)
mapshow(centroidsWorld(:,1),centroidsWorld(:,2),DisplayType="point")

Helper Functions

The detectWater helper function detects large bodies of water in a blocked image. Within the function:

• Identify dark pixels, which likely correspond to water, by using intensity thresholding. ROIs consist of pixels with a value less than the maxWaterPixelValue argument.

• Clean up the ROIs by filling holes using the imfill function and by shrinking the ROI using the imerode function.

• Retain only large ROIs that correspond to bodies of water by using the bwareafilt function. This example defines a large ROI as one containing at least 500 pixels.

function bw = detectWater(bs,maxWaterPixelValue)
    bw = bs.Data(:,:,1) < maxWaterPixelValue;
    bw = imfill(bw,"holes");
    bw = imerode(bw,ones(3));
    bw = bwareafilt(bw,[500 Inf]);
end

The detectBoats helper function detects boats on bodies of water using intensity thresholding. Within the function:

• Set all block pixels outside the mask of water ROIs to 0.

• Identify bright pixels, which likely correspond to boats, by using intensity thresholding. Detections consist of pixels with a value greater than the minBoatPixelValue argument.

• Measure the area and centroid of all detections using the regionprops function.

• Remove detections with an area smaller than the minBoatArea argument.

• Remove detections on the border by using the clearBorderData helper function, which is defined at the end of this example.

• Return the centroids of the detections.

function detections = detectBoats(block,waterMaskBlock, ...
    minBoatPixelValue,minBoatArea)
    
    % Resize the mask block to match the size of the image block
    waterMaskBlock = imresize(waterMaskBlock,size(block.Data,[1 2]),"nearest");
    
    % Set pixels outside the ROI to 0 and detect objects using intensity
    % thresholding
    imMasked = im2gray(block.Data);
    imMasked(~waterMaskBlock) = 0;
    bw = imMasked > minBoatPixelValue;

    % Measure the area and centroid of detected objects
    stats = regionprops(bw,["Area" "Centroid"]);

    % Remove detections with small area
    stats([stats.Area] < minBoatArea) = [];
    
    % Remove detections on the border
    stats = clearBorderData(stats,block);
    
    if isempty(stats)
        detections.centroidsXY = [];
    else
        blockCentroidsXY = vertcat(stats.Centroid);
        % Convert coordinates from row-column to X-Y order
        blockOffsetXY = fliplr(block.Start(1:2));
        detections.centroidsXY = blockCentroidsXY + blockOffsetXY;
    end
end

The clearBorderData helper function removes detections whose centroids are in the border area of a block.

function stats = clearBorderData(stats,block)

    % Convert coordinates from row-column to X-Y order
    lowEdgeXY = fliplr(block.BorderSize(1:2) + 1);
    highEdgeXY = fliplr(block.BlockSize(1:2) - block.BorderSize(1:2));
    
    % Find indices of centroids that are inside the current block, excluding 
    % the border region
    centroidsXY = reshape([stats.Centroid],2,[]);
    centroidsXY = round(centroidsXY);
    inBlock = centroidsXY(1,:)>=lowEdgeXY(1) & centroidsXY(1,:)<=highEdgeXY(1) ...
        & centroidsXY(2,:)>=lowEdgeXY(2) & centroidsXY(2,:)<=highEdgeXY(2);
    
    % Filter stats to retain only inside block detections
    stats = stats(inBlock);
end