Aplicar campo de deformación sin afectar las intensidades.

Estoy trabajando con Volúmenes (matrices 3D), y tengo campos de transformación geométrica y desplazamiento de la función malformado . Mi imagen en movimiento y la imagen fija son imágenes binarias (0s y 1s). Ahora quiero aplicar la misma transformación geométrica (o el campo de desplazamiento) a otra versión (esta versión no es una imagen binaria) de la imagen en movimiento. No quiero usar enano A medida que cambia los valores en la imagen. Por lo tanto, me gustaría aplicar la transformación sin afectar los valores en la imagen.

He escrito la siguiente función de MATLAB para hacer lo mismo usando el campo Desplazamiento:

    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

También he intentado usar el forma . A continuación está el código para eso.

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 funciones devuelven la misma salida que se espera. Pero no replican lo enano hace. Sería de gran ayuda si alguien puede señalar mi error o guiarme de otra manera.

Gracias.