Problem with rectifying perspective distortion

Good evening,

I have the image below and its vanishing points. Now the white paper on the image has to be transformed as if one was looking on it from above, in other words, the white paper has to be rectified, resulting in horizontal lines of the white paper.

I used the command ‘projective’, but the result is not that what it has to be, it still appears afterwards as a rectangle instead of a rectified image.

The theory says I have to fill in a 9x9 matrix using the homogenious vanishing point coordinates (points.png), invert that matrix and use the resulting H homography matrix to rectify the image. What comes out with the code below: (result.png)

clear variables; close all;

im_rgb = im2double(imread('DSC_1699.jpg'));
figure;
imshow(im_rgb,[],'InitialMagnification','fit');

%[xx,yy] = getpts
xx = [595.761146496815;2323.14968152866;4045.44267515924;5390.66560509554];
yy = [1755.91401273885;634.894904458599;3483.30254777070;1725.34076433121];

%   2
%       4
%1
%    3

p1 = [yy(1),xx(1)];
p2 = [yy(2),xx(2)];
p3 = [yy(3),xx(3)];
p4 = [yy(4),xx(4)];

hold on
plot(p1(2),p1(1),'*r');
plot(p2(2),p2(1),'*r');
plot(p3(2),p3(1),'*y');
plot(p4(2),p4(1),'*y');

line1 = cross([p1,1],[p2,1]);
line2 = cross([p3,1],[p4,1]);

line3 = cross([p4,1],[p2,1]);
line4 = cross([p3,1],[p1,1]);

vp1 = cross(line1,line2);
vp2 = cross(line3,line4);

plot(vp1(2)/vp1(3),vp1(1)/vp1(3),'*k');
plot(vp2(2)/vp2(3),vp2(1)/vp2(3),'*k');

Matrix = zeros(9,9);
Matrix(1,1) = vp1(1)/vp1(3);
Matrix(1,2) = vp1(2)/vp1(3);
Matrix(1,3) = vp1(3)/vp1(3);

Matrix(2,4) = vp1(1)/vp1(3);
Matrix(2,5) = vp1(2)/vp1(3);
Matrix(2,6) = vp1(3)/vp1(3);

Matrix(3,7) = vp1(1)/vp1(3);
Matrix(3,8) = vp1(2)/vp1(3);
Matrix(3,9) = vp1(3)/vp1(3);

Matrix(4,1) = vp2(1)/vp2(3);
Matrix(4,2) = vp2(2)/vp2(3);
Matrix(4,3) = vp2(3)/vp2(3);

Matrix(5,4) = vp2(1)/vp2(3);
Matrix(5,5) = vp2(2)/vp2(3);
Matrix(5,6) = vp2(3)/vp2(3);

Matrix(6,7) = vp2(1)/vp2(3);
Matrix(6,8) = vp2(2)/vp2(3);
Matrix(6,9) = vp2(3)/vp2(3);

Matrix(7,3) = 1;
Matrix(8,6) = 1;
Matrix(9,9) = 1;

u = transpose([1,0,0,0,1,0,0,0,0.0001]);

H = Matrix\u;

%Construct H matrix

Hmatrix = [H(1),H(2),H(3)
           H(4),H(5),H(6)
           H(7),H(8),H(9)];

tform = projective2d(Hmatrix);

imout = imwarp(im_rgb,tform);

figure;
imshow(imout,[],'InitialMagnification','fit');