Warp Big Image

This example shows how to apply a geometric transformation (warping) to a big image.

Applying a geometric transformation to an image is a key step in many image processing applications like image registration. You can use imwarp to warp coarse images that fit in-memory. For high-resolution big images that do not fit in memory, warp blocks of the image. Set the spatial referencing of the warped image to preserve characteristics of the image such as pixel extents.

Create a bigimage. This examples uses a modified version of image "tumor_091.tif" from the CAMELYON16 data set. The original image is a training image of a lymph node containing tumor tissue. The original image has eight resolution levels, and the finest level has resolution 53760-by-61440. The modified image has only three coarse resolution levels. The spatial referencing of the modified image has been adjusted to enforce a consistent aspect ratio and to register features at each level.

bim = bigimage('tumor_091R.tif');

Apply Geometric Transformation to Coarse Image

Create an affine2d object that stores information about an affine geometric transformation. This transformation applies translation and shear.

tform = affine2d([
    0.99 0.01 0
    0.17 0.98 0
    120 -30 1]);

Get the image at the coarsest resolution level.

imCoarse = getFullLevel(bim);

Warp the coarse image by using imwarp. Display the original image and the warped image in a montage.

imCoarseWarped = imwarp(imCoarse,tform);
montage({imCoarse,imCoarseWarped})

Create Spatial Referencing for Warped Fine Image

Before applying the geometric transformation to the image at a fine resolution level, calculate the spatial referencing of the big image after the warping. Use this spatial referencing when transforming blocks

Get the pixel extent of the original image from its spatial referencing information.

srcLevel = bim.FinestResolutionLevel;
inRef = bim.SpatialReferencing(srcLevel);
inPixelExtent = [inRef.PixelExtentInWorldX,inRef.PixelExtentInWorldY];

Calculate the output horizontal and vertical spatial limits when the transformation is applied.

[xout,yout] = outputLimits(tform,inRef.XWorldLimits,inRef.YWorldLimits);

Calculate the size of the output image that preserves the pixel extent. Specify the image size in the format [numrows, numcols].

outImgSize = [ceil(diff(yout)/inPixelExtent(2)),ceil(diff(xout)/inPixelExtent(1))];

Create an imref2d object that stores the spatial referencing information of the warped image. Set the world limits and image size of the warped image.

outRef = imref2d(outImgSize,xout,yout)
outRef = 
  imref2d with properties:

           XWorldLimits: [120.5800 6275]
           YWorldLimits: [-29.5050 4.9241e+03]
              ImageSize: [4954 6155]
    PixelExtentInWorldX: 0.9999
    PixelExtentInWorldY: 0.9999
    ImageExtentInWorldX: 6.1544e+03
    ImageExtentInWorldY: 4.9536e+03
       XIntrinsicLimits: [0.5000 6.1555e+03]
       YIntrinsicLimits: [0.5000 4.9545e+03]

Specify the origin of the output as the top left pixel.

outOrigin = [xout(1),yout(1)];

Calculate the corresponding output pixel dimensions.

outPixelExtent = [outRef.PixelExtentInWorldX,outRef.PixelExtentInWorldY];
halfPixWidth = outPixelExtent/2;

Apply Block-Wise Warping to Fine Image

Create a writable bigimage by specifying the output spatial referencing information. Specify a block size that is large enough to use memory efficiently.

outBlockSize = [1024 1024];
bwarped = bigimage(outRef,bim.Channels,bim.ClassUnderlying, ...
    'BlockSize',outBlockSize);

Loop through the output image, one block at a time. For each output block:

  1. Find the coordinates of the four corners of the output block.

  2. Inverse map these coordinates back to the input to get the input (source) region.

  3. Read the contents of the input region.

  4. Create spatial referencing describing the input region.

  5. Calculate the output block content by using imwarp.

  6. Write the output block to the output image by using the setBlock function.

If you have Parallel Computing Toolbox™, then you can replace the outer for statement with a parfor statement to run the loop in parallel.

colStarts = 1:outBlockSize(2):outRef.ImageSize(2);
for cInd = 1:numel(colStarts)
    colStart = colStarts(cInd);
    for rStart = 1:outBlockSize(1):outRef.ImageSize(1)
        
        % Center of top left pixel of this block in world units
        xyStart = [colStart,rStart].* outPixelExtent;
        xyStart = xyStart + outOrigin;
        
        % Center of bottom right pixel of this block in world units
        bsize = fliplr(bwarped.BlockSize); % (r,c) -> (x,y)
        xyEnd = ([colStart,rStart] + (bsize-1)).*outPixelExtent;
        xyEnd = xyEnd + outOrigin;
        
        % Output bounding box
        outbbox = [xyStart
            xyStart(1) xyEnd(2)
            xyEnd(1) xyStart(2)
            xyEnd];
        
        % Corresponding spatial reference which describes this rectangle
        outRegionRef = imref2d(bsize);
        outRegionRef.XWorldLimits = [xyStart(1)-halfPixWidth(1),xyEnd(1)+halfPixWidth(1)];
        outRegionRef.YWorldLimits = [xyStart(2)-halfPixWidth(2),xyEnd(2)+halfPixWidth(2)];
        
        % Get corresponding input region. The region is not rectangular if the transformation
        % includes shear
        inRegion = transformPointsInverse(tform,outbbox);   
        
        % Clamp to image extents (working with pixel centers)
        xcenterLims = inRef.XWorldLimits + [halfPixWidth(1)-halfPixWidth(1)];
        ycenterLims = inRef.YWorldLimits + [halfPixWidth(2)-halfPixWidth(2)];
        inRegion(:,1) = max(inRegion(:,1),xcenterLims(1));
        inRegion(:,1) = min(inRegion(:,1),xcenterLims(2));
        inRegion(:,2) = max(inRegion(:,2),ycenterLims(1));
        inRegion(:,2) = min(inRegion(:,2),ycenterLims(2));
        
        % Find the input bounding box (still using pixel centers)
        inbboxStart = [min(inRegion(:,1)) min(inRegion(:,2))];
        inbboxEnd = [max(inRegion(:,1)) max(inRegion(:,2))];
        
        % Read corresponding input region 
        inputRegion = getRegion(bim,srcLevel,inbboxStart,inbboxEnd);
        
        % Input region's spatial referencing
        inRegionRef = imref2d(size(inputRegion));
        % Convert to pixel edges from pixel centers
        inRegionRef.XWorldLimits = [inbboxStart(1)-halfPixWidth(1),inbboxEnd(1)+halfPixWidth(2)];
        inRegionRef.YWorldLimits = [inbboxStart(2)-halfPixWidth(1),inbboxEnd(2)+halfPixWidth(2)];
        
        warpedBlock = imwarp(inputRegion,inRegionRef,tform,'OutputView',outRegionRef);
        
        % Set the block data in the output bigimage
        setBlock(bwarped,1,xyStart,warpedBlock);
    end
end

Display the warped image.

bigimageshow(bwarped)

See Also

| | | | | |