Set Up Spatial Referencing for Blocked Images

This example shows how to set up and verify the spatial referencing information of a blockedImage object.

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 Camelyon16 image.

imageDir = fullfile('I:\','Camelyon16');
if ~exist(imageDir,'dir')
    mkdir(imageDir);
end

To download the image, go to the Camelyon17 website and click the first "CAMELYON16 data set" link. Open the "training" then "tumor" directory. Download the "tumor_091.tif" file and move the file to the directory specified by the imageDir variable.

Explore Default Spatial Referencing

Create a 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.

linkaxes([hl.Parent,hr.Parent]);

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.

Get the ImageDescription metadata field at the finest resolution level of the blockedImage object.

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" &quot

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