Process Large GeoTIFF File as Blocked Image

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:

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:

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