Aplique o campo de deformação sem afetar as intensidades

Question

Estou trabalhando com Volumes (matrizes 3D) e tenho campos de transformação e deslocamento geométricos da função imregtform . Minha imagem em movimento e imagem fixa são imagens binárias (0s e 1s). Agora eu quero aplicar a mesma transformação geométrica (ou o campo de deslocamento) para outra versão (esta versão não é imagem binária) da imagem em movimento. Eu não quero usar imwarp , pois altera os valores na imagem. Então, eu gostaria de aplicar a transformação sem afetar os valores na imagem.

Eu escrevi a seguinte função do MATLAB para fazer o mesmo usando o campo Displacement:

    function [ Volumeafterd ] = ApplyDeformation2Volume( Volume, Deformation , OPFilename)
%   Detailed explanation goes here

%% Check the Size
SizeD = size(Deformation);
SizeV = size(Volume);
Volumeafterd =  zeros(SizeV(1),SizeV(2),SizeV(3));
if SizeV(1) == SizeD(1) & SizeV(2) == SizeD(2) & SizeV(3) == SizeD(3)
    %% Apply the Deformation
    fnd = find(Volume~=0);
    counterd1=0;
    for it=1:length(fnd)
        [x,y,z] = ind2sub(SizeV,fnd(it));
        d1 = squeeze(Deformation(x,y,z,:));
        x1=round(x+d1(1));y1=round(y+d1(2));z1=round(z+d1(1));
        if x1<1 | y1<1|z1<1
            fprintf('Non positive index (D1) for x %d, y %d, z %d adding %d,%d,%d to get x %d,y %d,z %d\n'...
                ,x,y,z,d1(1),d1(2),d1(3),x1,y1,z1);
        end
        if (x1>=1 & y1>=1 & z1>=1)
            if (Volumeafterd(x1,y1,z1)~=0)
                %fprintf('Over writing %d %d %d\n',x1,y1,z1);
                counterd1 = counterd1+1;
            end
            Volumeafterd(x1,y1,z1) = Volume(x,y,z);
        end
    end
    fprintf("Overwritten counter : %d\n",counterd1);
        else
    fprintf("The size of Deformation field and the volume does not match\nThe output will be empty Volume \n");
end
end

Eu também tentei usar o tform . A seguir está o código para isso.

function [ Volumeafterd ] = Applytform2Volume( Volume, TForm, OPFilename)

szV = size(Volume);
Volumeafterd = zeros(szV(1),szV(2),szV(3));

counter=0;

%find the non zeros labelled voxels
fndnz = find(Volume~=0);
for i = 1:length(fndnz)
    %for  i1 = 1:szV(1)
    %    for  i2= 1:szV(2)
    %       for  i3 = 1:szV(3)
    [i1,i2,i3] = ind2sub(size(Volume),fndnz(i));
    [nx,ny,nz ] =   transformPointsForward(TForm,i1,i2,i3);
    nr = round([nx,ny,nz ]);
    nx=nr(1);  ny=nr(2); nz=nr(3);
    %if nx<1 | ny<1|nz<1
    %fprintf('Non positive index (D1) for x %d, y %d, z %d to get x %d,y %d,z %d\n'...
    %,i1,i2,i3,nx,ny,nz);
    %else
    try
        Volumeafterd(nx,ny,nz) = Volume(i1,i2,i3);
    catch
        % fprintf('Non positive index (D1) for x %d, y %d, z %d to get x %d,y %d,z %d\n'...
        %,i1,i2,i3,nx,ny,nz);
        counter = counter+1;
    end
end
%end
%       end
%  end
%end

fprintf('Counter %d for %s \n',counter,OPFilename);
end

Ambas as funções retornam a mesma saída esperada. Mas eles não replicam o que o imwarp faz. Seria de grande ajuda se alguém pudesse apontar meu erro ou me orientar de outra maneira.

Obrigado.