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
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
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.
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
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
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 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])
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
figure mapshow(imLowRes,RLowRes) mapshow(centroidsWorld(:,1),centroidsWorld(:,2),DisplayType="point")
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
Clean up the ROIs by filling holes using the
imfillfunction and by shrinking the ROI using the
Retain only large ROIs that correspond to bodies of water by using the
bwareafiltfunction. 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
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
Measure the area and centroid of all detections using the
Remove detections with an area smaller than the
Remove detections on the border by using the
clearBorderDatahelper 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
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
makeMultiLevel2D(Image Processing Toolbox) |
selectBlockLocations(Image Processing Toolbox) |
regionprops(Image Processing Toolbox) |
bigimageshow(Image Processing Toolbox) |
apply(Image Processing Toolbox)
blockedImage(Image Processing Toolbox)