粗い解像度レベルおよび細かい解像度レベルでブロック化されたイメージのワープ

この例では、ブロック化されたイメージに幾何学的変換 (ワーピング) を適用する方法を説明します。

イメージへの幾何学的変換の適用は、イメージレジストレーションなどのような多数のイメージ処理アプリケーションの重要な手順です。メモリに収まる粗いイメージをワープさせるには、imwarpを使用します。大きいイメージ (メモリに収まらない高解像度イメージ) の場合は、ブロック化されたイメージを使用します。ワープされたイメージの空間参照を、ピクセル範囲などのイメージの特性が維持されるように設定します。

CAMELYON16 データセットのイメージ "tumor_091.tif" の変更したバージョンを使用してブロック化されたイメージを作成します。元のイメージは、腫瘍組織が含まれるリンパ節の学習イメージです。元のイメージには 8 つの解像度レベルがあり、最も細かいレベルの解像度は 53760 x 61440 です。変更したイメージには、3 つの粗い解像度レベルのみが含まれています。変更したイメージの空間参照は、縦横比が一定に維持され、各レベルで特徴がレジストレーションされるように調整されています。

bim = blockedImage("tumor_091R.tif");

粗いイメージへの幾何学的変換の適用

アフィン幾何学的変換に関する情報を格納する affinetform2d オブジェクトを作成します。この変換は平行移動とせん断を適用します。

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

最も粗い解像度レベルのイメージを取得します。

imCoarse = gather(bim);

imwarp を使用して粗いイメージをワープさせます。元のイメージとワープされたイメージをモンタージュに表示します。

imCoarseWarped = imwarp(imCoarse,tform);
figure
imshow(imCoarseWarped)

ワープされた細かいイメージの空間参照の作成

細かい解像度レベルのイメージに幾何学的変換を適用する前に、ワーピング後のブロック化されたイメージの空間参照を計算します。ブロックを変換するときに、この空間参照を使用します。

その空間参照情報から元のイメージのピクセル範囲を取得します。

inPixelExtent = (bim.WorldEnd(1,:)-bim.WorldStart(1,:))./bim.Size(1,:);

変換を適用するときの出力の水平方向と垂直方向の空間範囲を計算します。

yWorldLimits = [bim.WorldStart(1,1),bim.WorldEnd(1,1)];
xWorldLimits = [bim.WorldStart(1,2),bim.WorldEnd(1,2)];
[xout,yout] = outputLimits(tform,xWorldLimits,yWorldLimits);

ピクセル範囲を維持する出力イメージのサイズを計算します。[numrows, numcols] の形式でイメージサイズを指定します。

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

ワープされたイメージの空間参照情報を格納します。ワープされたイメージのワールド座標範囲とイメージサイズを設定します。

outWorldStart = [yout(1) xout(1)];
outWorldEnd = [yout(2) xout(2)];

対応する出力ピクセルの次元を計算します。

outPixelExtent = (outWorldEnd-outWorldStart)./outImgSize;
halfPixelWidth = outPixelExtent/2;

細かいイメージへのブロック単位のワーピングの適用

出力の空間参照情報を指定して、書き込み可能なブロック化されたイメージを作成します。メモリを効率的に使用できる十分な大きさのブロックサイズを指定します。

outBlockSize = [1024 1024 3];
bwarped = blockedImage([],[outImgSize 3],outBlockSize,uint8(0), ...
    WorldStart=outWorldStart,WorldEnd=outWorldEnd,Mode="w");

出力イメージを一度に 1 ブロックずつループ処理します。各出力ブロックに対して以下を行います。

出力ブロックの 4 つの隅の座標を求めます。
これらの座標を入力に逆マッピングして、入力 (ソース) 領域を取得します。
入力領域のコンテンツを読み取ります。
入力領域を記述する空間参照を作成します。
imwarp を使用して出力ブロックのコンテンツを計算します。
関数 setBlock を使用して出力ブロックを出力イメージに書き込みます。

Parallel Computing Toolbox™ がある場合、外側の for ステートメントを parfor ステートメントに置き換えて、このループを並列実行することができます。

inYWorldLimits = [bim.WorldStart(1,1), bim.WorldEnd(1,1)];
inXWorldLimits = [bim.WorldStart(1,2), bim.WorldEnd(1,2)];

for rBlockInd = 1:bwarped.SizeInBlocks(1)
    for cBlockInd = 1:bwarped.SizeInBlocks(2)
        
        % Get the indices of the block
        blockSub = [rBlockInd cBlockInd 1];
        
        % Convert the block indices to pixel subscripts to get the 
        % subscript of the top left pixel. Based on the block size, 
        % get the subscripts of the bottom right pixel
        blockStartSub = blocksub2sub(bwarped,blockSub);
        blockEndSub = blockStartSub + outBlockSize - 1;
        
        % Convert the pixel indices to world coordinates. The world 
        % coordinates indicate the center of the top left and bottom 
        % right pixels of the block in world units
        blockStartWorld = sub2world(bwarped,blockStartSub(1:2));
        blockEndWorld = sub2world(bwarped,blockEndSub(1:2));
        
        % Spatial referencing information for this block (Note: spatial
        % referencing is in x-y order, while blockStart etc are in y-x
        % order).
        outRegionRef = imref2d(fliplr(outBlockSize(1:2)));
        % Expand the region outwards by half a pixel to align with the
        % outer edge of the block
        outRegionRef.YWorldLimits = [blockStartWorld(1)-halfPixelWidth(1),...
                                     blockEndWorld(1)+halfPixelWidth(1)];        
        outRegionRef.XWorldLimits = [blockStartWorld(2)-halfPixelWidth(2),...
                                     blockEndWorld(2)+halfPixelWidth(2)];
        
        % Output bounding box in world coordinates in x-y order
        outbbox = [
            fliplr(blockStartWorld) % top left
            blockStartWorld(2) blockEndWorld(1) % bottom left
            blockEndWorld(2) blockStartWorld(1) % top right
            fliplr(blockEndWorld)   % bottom right
            ];
        
        % Get corresponding input region. Note: This region need NOT be
        % rectangular if the transformation includes shear
        inRegion = transformPointsInverse(tform,outbbox);   
        
        % Clamp region to image extents
        inRegion(:,2) = max(inRegion(:,2),inYWorldLimits(1));
        inRegion(:,2) = min(inRegion(:,2),inYWorldLimits(2));       
        inRegion(:,1) = max(inRegion(:,1),inXWorldLimits(1));
        inRegion(:,1) = min(inRegion(:,1),inXWorldLimits(2));
        
        % Find the corresponding input bounding box
        inbboxStart = [min(inRegion(:,1)) min(inRegion(:,2))];
        inbboxEnd   = [max(inRegion(:,1)) max(inRegion(:,2))];
        
        % Move to y-x (row-col) order
        inbboxStart = fliplr(inbboxStart);
        inbboxEnd = fliplr(inbboxEnd);
        
        % Convert to pixel subscripts
        inbboxStartSub = world2sub(bim,inbboxStart);
        inbboxEndSub = world2sub(bim,inbboxEnd);
        
        % Read corresponding input region 
        inputRegion = getRegion(bim,inbboxStartSub,inbboxEndSub);
               
        % Get the input region's spatial referencing
        inRegionRef = imref2d(size(inputRegion));
        
        % Convert the actual region pixel's centers back to world
        % coordinates
        inbboxStart = sub2world(bim,inbboxStartSub);
        inbboxEnd = sub2world(bim,inbboxEndSub);
        
        % Convert to pixel edges from pixel centers
        inRegionRef.YWorldLimits = [inbboxStart(1)-halfPixelWidth(1),...
                                    inbboxEnd(1)+halfPixelWidth(2)];
        inRegionRef.XWorldLimits = [inbboxStart(2)-halfPixelWidth(1),...
                                    inbboxEnd(2)+halfPixelWidth(2)];
        
        % Warp this block
        warpedBlock = imwarp(inputRegion,inRegionRef,tform,OutputView=outRegionRef);
        
        % Set the block data in the output blocked image
        setBlock(bwarped,blockSub,warpedBlock);
        
    end
end

ワープされたイメージを表示します。

bwarped.Mode = "r";
figure
bigimageshow(bwarped)

参考