Самостоятельная функция преобразования подобия изображения

Question

Это домашняя задача, в которой мы пишем самоопределяемые функции, используя лежащую в основе математику для получения прогнозируемого результата. Возьмите набор соответствующих пиксельных координат из двух изображений, где одно изображение («движущееся») является повернутой, масштабированной и переведенной версией первого. Оцените параметры преобразования и примените преобразование ко всем пикселям первого изображения. Это должно создать фиксированное изображение, из которого были взяты соответствующие пиксельные координаты.

Я предоставляю наш код, а также движущиеся, фиксированные и выходные изображения. Встроенная функция и вывод нашего кода должны быть одинаковыми. Визуально кажется, что масштабирование и перемещение отключены, а также, возможно, центр вращения. Мой партнер и я озадачены и ищем подсказки относительно того, как исправить код, чтобы получить правильный результат. Любые рекомендации или предложения приветствуются.

      moving_pts=[104,85;135,84;107,145;154,136;110,22;183,148;62,124;99,186];
fixed_pts=[95,97;127,102;88,158;138,158;113,37;161,174;49,128;72,198];
moving_image = imread('BrainProtonDensitySliceR10X13Y17.png');
output=xformed_image(fixed_pts, moving_pts, moving_image)


function output=xformed_image(fixed_pts, moving_pts, moving_image)
   %get the size for moving image
   [x,y,z]=size(moving_image);
  
   %pre-allocate the image for output
   output=zeros(x,y,z);
   %get the translation parameters
   theta=getTheta(moving_pts, fixed_pts);
  
   %rotate the image and get the intensity value
   output=getRotate(moving_image, theta);
end

function theta=getTheta(g, f)
   % g is the points from moving image
   % f is the correspondences points from fixed image
   %get the size
   [row,col]=size(g);
   % pre-allocate the size for theta
   theta=zeros(4,1);
   g=g+1;
   f=f+1;
  
   %pre-allocate the size for X
   X=zeros(4,4);
   for i=1:row
       Xi=[1,0,g(i,1),g(i,2);
           0,1,g(i,2),-g(i,1)]; 
       X=X+Xi'*Xi;
   end
   %pre-allocate the size for Y
   Y=zeros(4,1);
   for i=1:row
       Xi=[1,0,g(i,1),g(i,2);
           0,1,g(i,2),-g(i,1)]; 
       Y=Y+Xi'*f(i,:)';
   end
   theta=inv(X)*Y;
end

function R=getRotate(m, theta)
   %m is the moving image
   %theta is the translation parameters
  
   %get the size of moving image
   [x,y,z]=size(m);
    %pre-allocate the size for rotation ans shifting
   R=zeros(x, y, z);
  
   R=uint8(R);
  
   %shifting matrix
   shiftMatrix=[x/2;y/2];
  
   %scaling matrix for Similarity
   scaleMatrix=[theta(1); theta(2)];
  
   %rotate matrix
   rotateMatrix=[theta(3),theta(4); -theta(4), theta(3)];
  
   %iterate each channel
   for d=1:z
       for i=1:x
           for j=1:y
               %new image point set
               newImagePoint=[i;j];
              
               %Step1: shift to the origin
               shiftPoint=newImagePoint-shiftMatrix;
              
               %Step2: rotate based on the translation parameters
               rotatePoint=scaleMatrix+rotateMatrix'*shiftPoint;
              
               %Step3: shifted the image back to the original center
               originalImagePoint=rotatePoint+shiftMatrix; %final (x,y)
               originalImageX=uint8(originalImagePoint(1));
               originalImageY=uint8(originalImagePoint(2));
 
               %Step4: interpolate to get the intensity value
               if(originalImageX>=1 && originalImageX<x && originalImageY>=1 && originalImageY<y)
                   R(i,j,d)=getInterpolate(m, d, originalImagePoint);
                   %R(i, j, d)=m(originalImageX,originalImageY,d);
               end
           end
       end
   end
end

function intensity=getInterpolate(m,channel,originalImagePoint)
   %m is the moving image
   %rotatedBackPoint is the rotated back x,y value
  
   x=uint8(originalImagePoint(1));
   y=uint8(originalImagePoint(2));
   lambda=originalImagePoint(1)-double(x);
   mu=originalImagePoint(2)-double(y);
  
   f1=mu*double(m(x,y+1,channel))+(1-mu)*double(m(x,y,channel));
   f2=mu*double(m(x+1,y+1,channel))+(1-mu)*double(m(x+1,y,channel));
  
   intensity=lambda*f2+(1-lambda)*f1;
  
   intensity=uint8(intensity);
end

Входное изображение

Фиксированное изображение, которое должно выглядеть на выходе, выглядит следующим образом:

Наш вывод такой: