Set Up Spatial Referencing for Blocked Images
This example shows how to set up and verify the spatial referencing information of a
Spatial Referencing in Blocked Images
Blocked images work with multiresolution images, in which image data of a scene is stored as a set of images at different resolution levels. Blocked images assume that the spatial extents of each level are the same, in other words, that all levels cover the same physical area in the real world. The first step in working with a large multiresolution image is to validate this assumption.
Download Blocked Image Data
This example uses one image from the Camelyon16 data set. This data set contains 400 whole-slide images (WSIs) of lymph nodes, stored as multiresolution TIF files that are too large to be loaded into memory.
Create a directory to store the image.
imageDir = fullfile(tempdir,'Camelyon17');
To download the image
tumor_091.tif for use in this example, go to the GigaDB page for the data set. On the GigaDB page, scroll down to the table and open the Files tab. Navigate the table to find the entry for
tumor_091.tif, and click on the file name to download the image.
After you download the image, move the file to the directory specified by the
Explore Default Spatial Referencing
blockedImage object with the default spatial referencing information. By default, a blocked image sets the spatial referencing of each level to have the same world extents as the finest layer. The finest layer is the layer that has the highest resolution and the most pixels.
fileName = fullfile(imageDir,'tumor_091.tif'); bim = blockedImage(fileName);
Display the spatial referencing information at the finest level. The image size (specified by the
Size property) matches the extents in world coordinates. Notice that the default image coordinate system puts the center of the first pixel at (1,1). Because the pixel extents are 1 unit wide in each dimension, the left edge of the first pixel starts at (0.5,0.5).
finestLevel = 1; finestStart = bim.WorldStart(finestLevel,:) finestEnd = bim.WorldEnd(finestLevel,:)
finestStart = 0.5000 0.5000 0.5000 finestEnd = 1.0e+04 * 5.3761 6.1441 0.0003
Display the spatial referencing information at the coarsest level. The world extents are the same as the finest level, but the coarse image size is only 512-by-512 pixels. Effectively, each pixel in this coarse level corresponds to a 105-by-120 block of pixels in the finest resolution.
coarsestLevel = bim.NumLevels; coarsestStart = bim.WorldStart(coarsestLevel,:) coarsestEnd = bim.WorldEnd(coarsestLevel,:)
coarsestStart = 0.5000 0.5000 0.5000 coarsestEnd = 1.0e+04 * 5.3761 6.1441 0.0003
Verify Aspect Ratio
Display the image size and aspect ratio at each level. The aspect ratio is not consistent, which indicates that levels do not all span the same world area. Therefore, the default assumption is incorrect for this image.
t = table((1:8)',bim.Size(:,1),bim.Size(:,2), ... bim.Size(:,1)./bim.Size(:,2), ... 'VariableNames',["Level" "Height" "Width" "Aspect Ratio"]); disp(t)
Level Height Width Aspect Ratio _____ ______ _____ ____________ 1 53760 61440 0.875 2 27136 30720 0.88333 3 13824 15360 0.9 4 7168 7680 0.93333 5 3584 4096 0.875 6 2048 2048 1 7 1024 1024 1 8 512 512 1
Display Layers to Compare Spatial Extents
Display the blocked image by using the
bigimageshow function. Display the coarsest resolution level.
figure subplot(1,2,1); hl = bigimageshow(bim,'ResolutionLevel',coarsestLevel); title('Coarsest Resolution Level (8)') %
Display image data at the default resolution level in the same figure window. By default,
bigimageshow selects the level to display based on the screen resolution and the size of the displayed region.
subplot(1,2,2); hr = bigimageshow(bim); title('Default Resolution Level') %
Ensure that both displays show the same extents.
Check Default Spatial Referencing
Zoom in on a feature.
xlim([45000 50000]); ylim([12000 17000]);
Change the resolution level of the image on the right side of the figure window. At level 6, the features look aligned with the coarsest level.
hr.ResolutionLevel = 6; title('Level 6'); snapnow %
At level 1, the features are not aligned. Therefore, level 1 and level 8 do not span the same world extents.
hr.ResolutionLevel = 1; title('Level 1'); snapnow %
Get Spatial Extents from Blocked Image Metadata
Usually the original source of the data has spatial referencing information encoded in its metadata. For the Camelyon16 data set, spatial referencing information is stored as XML content in the
ImageDescription metadata field at the finest resolution level. The XML content has a DICOM_PIXEL_SPACING attribute for each resolution level that specifies the pixel extents.
ImageDescription metadata field at the finest resolution level of the
binfo = imfinfo(bim.Source); binfo = binfo(1).ImageDescription;
Search the content for the string "DICOM_PIXEL_SPACING". There are nine matches found. The second instance of the attribute corresponds to the pixel spacing of the finest level. The last instance of the attribute corresponds to the pixel spacing of the coarsest level.
indx = strfind(binfo,"DICOM_PIXEL_SPACING");
Store the pixel spacing at the finest level. To extract the value of the pixel spacing from the XML text, visually inspect the text following the second instance of the "DICOM_PIXEL_SPACING" attribute.
disp(binfo(indx(2):indx(2)+100)) pixelSpacing_L1 = 0.000227273;
DICOM_PIXEL_SPACING" Group="0x0028" Element="0x0030" PMSVR="IDoubleArray">"0.000227273" &qu
Similarly, store the pixel spacing at the coarsest level. To extract the value of the pixel spacing from the XML text, visually inspect the text following the last instance of the "DICOM_PIXEL_SPACING" attribute.
disp(binfo(indx(end):indx(end)+100)) pixelSpacing_L8 = 0.0290909;
DICOM_PIXEL_SPACING" Group="0x0028" Element="0x0030" PMSVR="IDoubleArray">"0.0290909" "
Set Spatial Extents
Calculate the relative pixel width between level 8 and level 1.
pixelDims = pixelSpacing_L8/pixelSpacing_L1;
The finest level has the reference spatial extent. Calculate the corresponding extents for level 8 with respect to the extents of level 1.
worldExtents = bim.Size(8,1:2).*pixelDims;
Update the spatial referencing of level 8.
bim.WorldEnd(8,1:2) = worldExtents(2);
Verify Alignment with Custom Spatial Referencing
Redisplay the data to confirm alignment of the key feature. Show level 8 on the left side and level 1 on the right side.
hl.CData = bim; hl.ResolutionLevel = 8; snapnow hr.CData = bim; hr.ResolutionLevel = 1; snapnow