% исходные данные
clear all;
M=400; N=500; H=120; L=100; Alfa=0;
Tetaz=-30; Tetay=-30; Tetax=0;
Tz=10; Ty=30; Tx=35;
Sz=0.8; Sy=1.2; Sx=1.5;

% формированаие исходного изображения
% Im2D=im2double(imread('c:\АБВ.bmp'));
% [M, N]=size(Im2D);
Im2D=zeros(M,N);
Im2D(round(M/2-L/2):round(M/2+L/2),round(N/2-H/2):round(N/2+H/2))=1;
Im2D=imrotate(Im2D,Alfa,'crop'); % в градусах против час стр
figure; imshow(Im2D); title('ISX IZOB');

% матрицы прямого геометрического преобразования
R3DX=[1  0                 0                 0
      0  cos(Tetax*pi/180) sin(Tetax*pi/180) 0
      0 -sin(Tetax*pi/180) cos(Tetax*pi/180) 0
      0  0                 0                 1]; % в рад по час стр
R3DY=[cos(Tetay*pi/180)    0                -sin(Tetay*pi/180) 0
      0                    1                 0                 0
      sin(Tetay*pi/180)    0                 cos(Tetay*pi/180) 0
      0                    0                 0                 1]; % в рад по час стр
R3DZ=[cos(Tetaz*pi/180) sin(Tetaz*pi/180) 0  0
     -sin(Tetaz*pi/180) cos(Tetaz*pi/180) 0  0
      0                 0                 1  0
      0                 0                 0  1]; % в рад по час стр
R3D=R3DZ*R3DY*R3DX;
S3D=[Sx  0   0   0
     0   Sy  0   0
     0   0   Sz  0
     0   0   0   1]; % растяжение раз
T3D=[1   0   0   0
     0   1   0   0
     0   0   1   0
     Tx  Ty  Tz  1]; % перенос вправо и вниз
TFormR3D=maketform('affine',R3D);
TFormS3D=maketform('affine',S3D);
TFormT3D=maketform('affine',T3D);
MP3D=[1   0   0   0
      0   1   0   0
      0   0   0   0
      0   0   0   1]; % ортографическая проекция
TFormMP3D=maketform('projective',MP3D);

% прямое геометрическое преобразование исх изображения
A3D=zeros(M*N,3);
k=0;
for j=1:M
    for i=1:N
        if Im2D(j,i)>=0.5
            k=k+1;
            A3D(k,1)=i-round(N/2)+1; A3D(k,2)=j-round(M/2)+1;
        end;
    end;
end;
B3D=A3D(1:k,:);
C3D1=tformfwd(TFormR3D, B3D);
C3D2=tformfwd(TFormS3D, C3D1);
C3D3=tformfwd(TFormT3D, C3D2);
D3D1=tformfwd(TFormMP3D, C3D1);
D3D2=tformfwd(TFormMP3D, C3D2);
D3D3=tformfwd(TFormMP3D, C3D3);
Im2DTransform1=zeros(M,N);
for i=1:k
    x=round(D3D1(i,1)+round(N/2)-1); y=round(D3D1(i,2)+round(M/2)-1);
    Im2DTransform1(y,x)=1;
end;
Im2DTransform2=zeros(M,N);
for i=1:k
    x=round(D3D2(i,1)+round(N/2)-1); y=round(D3D2(i,2)+round(M/2)-1);
    Im2DTransform2(y,x)=1;
end;
Im2DTransform3=zeros(M,N);
for i=1:k
    x=round(D3D3(i,1)+round(N/2)-1); y=round(D3D3(i,2)+round(M/2)-1);
    Im2DTransform3(y,x)=1;
end;

% вывод преобразованного изображения
figure; subplot(2,2,1); imshow(Im2D); title('Isx Izob');
subplot(2,2,2); imshow(Im2DTransform1); title('Proek Povorot');
subplot(2,2,3); imshow(Im2DTransform2); title('Proek Maschtab');
subplot(2,2,4); imshow(Im2DTransform3); title('Proek Peremesch');

% % матрицы обратного геометрического преобразования
% R3DXInv=[1  0                  0                 0
%          0  cos(Tetax*pi/180) -sin(Tetax*pi/180) 0
%          0  sin(Tetax*pi/180)  cos(Tetax*pi/180) 0
%          0  0                 0                 1]; % в рад по час стр
% R3DYInv=[cos(Tetay*pi/180)    0                 sin(Tetay*pi/180) 0
%          0                    1                 0                 0
%         -sin(Tetay*pi/180)    0                 cos(Tetay*pi/180) 0
%          0                    0                 0                 1]; % в рад по час стр
% R3DZInv=[cos(Tetaz*pi/180) -sin(Tetaz*pi/180) 0  0
%          sin(Tetaz*pi/180)  cos(Tetaz*pi/180) 0  0
%          0                 0                 1  0
%          0                 0                 0  1]; % в рад по час стр
% R3DInv=R3DZInv*R3DYInv*R3DXInv;
% S3DInv=[1/Sx  0     0     0
%         0     1/Sy  0     0
%         0     0     1/Sz  0
%         0     0     0     1]; % растяжение раз
% T3DInv=[1   0   0   0
%         0   1   0   0
%         0   0   1   0
%        -Tx -Ty -Tz  1]; % перенос вправо и вниз
% TFormR3DInv=maketform('affine',R3DInv);
% TFormS3DInv=maketform('affine',S3DInv);
% TFormT3DInv=maketform('affine',T3DInv);
 
% обратное геометрическое преобразование преоб изображения
C3D2Restore=tforminv(TFormT3D, C3D3);
C3D1Restore=tforminv(TFormS3D, C3D2Restore);
B3DRestore=tforminv(TFormR3D, C3D1Restore);
D3D3Restore=tformfwd(TFormMP3D, C3D2Restore);
D3D2Restore=tformfwd(TFormMP3D, C3D1Restore);
D3D1Restore=tformfwd(TFormMP3D, B3DRestore);
Im2DRestore3=zeros(M,N);
for i=1:k
    x=round(D3D3Restore(i,1)+round(N/2)-1); y=round(D3D3Restore(i,2)+round(M/2)-1);
    Im2DRestore3(y,x)=1;
end;
Im2DRestore2=zeros(M,N);
for i=1:k
    x=round(D3D2Restore(i,1)+round(N/2)-1); y=round(D3D2Restore(i,2)+round(M/2)-1);
    Im2DRestore2(y,x)=1;
end;
Im2DRestore1=zeros(M,N);
for i=1:k
    x=round(D3D1Restore(i,1)+round(N/2)-1); y=round(D3D1Restore(i,2)+round(M/2)-1);
    Im2DRestore1(y,x)=1;
end;

% вывод восстановленного изображения
figure; subplot(2,2,1); imshow(Im2DTransform3); title('Preob Izob');
subplot(2,2,2); imshow(Im2DRestore3); title('Obr Peremesch');
subplot(2,2,3); imshow(Im2DRestore2); title('Obr Maschtab');
subplot(2,2,4); imshow(Im2DRestore1); title('Obr Povorot');

% получение координат точек
figure; [XIsx, YIsx, PIsx]=impixel(Im2D) 
[XTransform3, YTransform3, PTransform3]=impixel(Im2DTransform3) 
[XRestore3, YRestore3, PRestore3]=impixel(Im2DRestore1)