[解決済み] Matlabでどのようにimwarp転送ポイント？

質問

Matlabを使って、画像をターゲット画像に変換しているところです。幾何学変換(tform)がある。

例えば、これは私の「tform」です。

    1.0235    0.0022   -0.0607         0
   -0.0276    1.0002    0.0089         0
   -0.0170   -0.0141    1.1685         0
   12.8777    5.0311  -70.0325    1.0000

Matlab2013では、imwarpを使って簡単に行うことができます。

%nii = 3D MR Image
I = nii.img; 
dii=nii.hdr.dime.pixdim(2:4);
Rfixed=imref3d(size(I),dii(2),dii(1),dii(3));    
new_img= imwarp(old_img, Rfixed, tform, 'OutputView', Rfixed);

imwarpで完璧に仕上がりました（画像は赤い肺）。

imwarpがどのように動作しているのか知る必要があるため、自分で関数を書いてみました。

function [new_img] = aff3d(old_img, tform, range_x, range_y, range_z)

   [U, V, W] = ndgrid(range_x, range_y, range_z);
   xyz = [reshape(U,[],1)';reshape(V,[],1)';reshape(W,[],1)'];
   xyz = [xyz; ones(1,size(xyz,2))];


   uvw = tform.T * xyz;
   % avoid homogeneous coordinate  
   uvw = uvw(1:3,:)';


   xi = reshape(uvw(:,1), length(range_x),length(range_y),length(range_z));
   yi = reshape(uvw(:,2), length(range_x),length(range_y),length(range_z));
   zi = reshape(uvw(:,3), length(range_x),length(range_y),length(range_z));

   old_img = single(old_img);
   new_img = interp3(old_img,yi,xi,zi,'linear');

   ii = find(isnan(new_img));
   if(~isempty(ii))
      new_img(ii) = 0;
   end
end

私の関数の結果( 詳細情報 がimwarpの出力と一致しません（赤い肺が正しい位置にない）、どなたか助けてください。

解決方法は？

Ander さんの提案のように、逆変換を掛けてみてください。

Tinv = tform.invert();
TinvMatrix = Tinv.T;

ということで、あなたのコードはこうなります。

function [new_img] = aff3d(old_img, tform, range_x, range_y, range_z)
   [U, V, W] = ndgrid(range_x, range_y, range_z);
   xyz = [reshape(U,[],1)';reshape(V,[],1)';reshape(W,[],1)'];
   xyz = [xyz; ones(1,size(xyz,2))];

   tformInv = invert(tform);
   uvw = tformInv.T * xyz;
   % avoid homogeneous coordinate  
   uvw = uvw(1:3,:)';
   xi = reshape(uvw(:,1), length(range_x),length(range_y),length(range_z));
   yi = reshape(uvw(:,2), length(range_x),length(range_y),length(range_z));
   zi = reshape(uvw(:,3), length(range_x),length(range_y),length(range_z));
   old_img = single(old_img);
   new_img = interp3(old_img,yi,xi,zi,'linear');
   ii = find(isnan(new_img));
   if(~isempty(ii))
      new_img(ii) = 0;
   end
end

あなたのコードでは、old_img の中で補間して、ワープした new_img を見つけようとしています。これは、出力画像空間から入力画像空間へマッピングする逆マッピングを使用したいことを示唆しています。点の順方向マッピングを使って古い画像を補間しているように見えますが、これは正しくありません。

http://blogs.mathworks.com/steve/2006/04/28/spatial-transforms-forward-mapping/ http://blogs.mathworks.com/steve/2006/05/05/spatial-transformations-inverse-mapping/

上記のリンクを使って、フォワードマッピングとインバースマッピングを見直すといいと思います。IMWARPはインバースマッピングを使用しています。

混乱の原因の一つは、人々が幾何学的変換を考えるとき、一般に点が古い_imageから新しい_imageにどのようにマッピングするかという順方向のマッピングで考えるからです。このため、アフィン3d変換のquot;T"特性は、順方向マッピングの観点から表現されています。

幾何学変換をソフトウェアで実装する必要がある場合、逆写像で実装する方がはるかに簡単で良いのです。imwarpがそうであり、imwarpの動作を再現しようとするときに、変換を反転させる必要があるのはそのためです。インバースマッピングのブログ記事を読むと、このアルゴリズムはまさにimwarpがやっていることです。

唯一の問題は、WorldLimitsがピクセルの離散グリッドで均等に分割できない場合に、IMWARPがデフォルトでない座標系で（デフォルトでない空間参照オブジェクトを使用して）何を行うかということです。この動作は任意であり、正しい動作はありません。IMWARPの動作は、要求された解像度（PixelExtentInWorld）を尊重し、この場合、ワールドリミットをわずかに調整することです。