Spatial Referencing in bigimage

bigimage works with multiresolution images, in which image data of a scene is stored as a set of images at different resolution levels. bigimage assumes 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 Big 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(tempdir,'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 bigimage with the default spatial referencing information. By default, bigimage 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.

bim = bigimage(fullfile(imageDir,'tumor_091.tif'));

Display the spatial referencing information at the finest level. The image size (specified by the ImageSize 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 = bim.FinestResolutionLevel;
finestLevelInfo = bim.SpatialReferencing(finestLevel)

finestLevelInfo = 

  imref2d with properties:

           XWorldLimits: [0.5000 6.1441e+04]
           YWorldLimits: [0.5000 5.3761e+04]
              ImageSize: [53760 61440]
    PixelExtentInWorldX: 1
    PixelExtentInWorldY: 1
    ImageExtentInWorldX: 61440
    ImageExtentInWorldY: 53760
       XIntrinsicLimits: [0.5000 6.1441e+04]
       YIntrinsicLimits: [0.5000 5.3761e+04]

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.CoarsestResolutionLevel;
disp(bim.SpatialReferencing(coarsestLevel))

  imref2d with properties:

           XWorldLimits: [0.5000 6.1441e+04]
           YWorldLimits: [0.5000 5.3761e+04]
              ImageSize: [512 512]
    PixelExtentInWorldX: 120
    PixelExtentInWorldY: 105
    ImageExtentInWorldX: 61440
    ImageExtentInWorldY: 53760
       XIntrinsicLimits: [0.5000 512.5000]
       YIntrinsicLimits: [0.5000 512.5000]

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.LevelSizes(:,1),bim.LevelSizes(:,2), ...
    bim.LevelSizes(:,1)./bim.LevelSizes(:,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 bigimage by using the bigimageshow function. Display the coarsest resolution level on the left side of a figure window.

figure
subplot(1,2,1);
hl = bigimageshow(bim,'ResolutionLevel',coarsestLevel);
title('Coarsest Resolution Level (8)')

Display image data at the default resolution level on the right side of a 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 Big 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 bigimage.

binfo = imfinfo(bim.DataSource);
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.SpatialReferencing(8).ImageSize.*pixelDims;

Update the spatial referencing of level 8.

bim.SpatialReferencing(8).XWorldLimits = [0.5 worldExtents(2)];
bim.SpatialReferencing(8).YWorldLimits = [0.5 worldExtents(1)];

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;
hr.CData = bim;
hr.ResolutionLevel = 1;