%————————————————————————–
% anamorphic.m
% Alexander Boxer 2012
%————————————————————————–
%—– Load the image
A = imread(‘~/Desktop/Holbein.jpg’);
%—– Flip the image to an xy coordinate system
A(:,:,1) = flipud(A(:,:,1));
A(:,:,2) = flipud(A(:,:,2));
A(:,:,3) = flipud(A(:,:,3));
%—– Set the painting size-parameters
Lx = 2095; % the length of the painting in mm
Ly = 2070; % the height of the painting in mm
x0 = 1083.2; % the location of the plumb-line
yEye = Ly/2; % vertical height of viewer’s eye
yFocus = yEye; % vertical height of eye’s focus on the painting
y0 = yEye; % central y-coordinate of the painting
R = 0.8449*Lx; % radial distance from (x0,y0) to the eye
alpha = 81.7108; % angle between the eye and the z-axis
beta = atand((yEye-yFocus)/R); % angle between eye-level and focus
resolution = 1; % 0.2 is highest resolution (and slowest)
%————————————————————————–
%—– Map the image onto an xy mesh-grid and mark key points
%————————————————————————–
[X1,Y1] = meshgrid(linspace(0,Lx,size(A,2)),linspace(0,Ly,size(A,1)));
%—– Reshape the xy-meshgrid and create a 3xN array of coordinate-data
N = numel(X1);
x1 = reshape(X1,1,N);
y1 = reshape(Y1,1,N);
z1 = zeros(1,N);
P0 = [x1;y1;z1];
%—– Mark the corners of the painting
c1 = [0;0;0]; % lower-left corner of the painting
c2 = [0;Ly;0]; % upper-left corner of the painting
c3 = [Lx;Ly;0]; % upper-right corner of the painting
c4 = [Lx;0;0]; % lower-right corner of the painting
C0 = [c1,c2,c3,c4];
%—– Mark points on the skull
s1 = [661;0;0]; % chin, lower left corner
s2 = [615;58;0]; % left eye-socket, left extremum
s3 = [1083.2;337;0]; % top mid-point
s4 = [1530;472;0]; % right extremum
s5 = [1210;257;0]; % corner of jawbone
S0 = [s1,s2,s3,s4,s5];
%————————————————————————–
%—– Convert the 3xN matrices to 4xN matrices
%————————————————————————–
P1 = cat(1,P0,ones(1,length(P0)));
C1 = cat(1,C0,ones(1,length(C0)));
S1 = cat(1,S0,ones(1,length(S0)));
%————————————————————————–
%—– Translate and rotate painting into viewing aspect
%————————————————————————–
disp(‘Transforming Coordinate Grid’)
%—– 1. Make a translation matrix to set the origin at (x0,y0)
Txy = [1,0,0,-x0;
0,1,0,-y0;
0,0,1,0;
0,0,0,1];
%—– 2. Make a matrix for rotation about the y-axis
Ry = @(theta) [cosd(theta),0,-sind(theta),0;
0,1,0,0;
sind(theta),0,cosd(theta),0;
0,0,0,1];
%—– 3. Make a matrix for translation along the z-axis
Tz = [1,0,0,0;
0,1,0,0;
0,0,1,-R;
0,0,0,1];
%—– 4. Make a matrix for rotation about the x-axis
Rx = @(theta) [1,0,0,0;
0,cosd(theta),-sind(theta),0;
0,sind(theta),cosd(theta),0;
0,0,0,1];
%—– Transform the coordinates
P = Rx(beta)*Tz*Ry(alpha)*Txy*P1;
C = Rx(beta)*Tz*Ry(alpha)*Txy*C1;
S = Rx(beta)*Tz*Ry(alpha)*Txy*S1;
%————————————————————————–
%—– Set the Dimensions of the Projection-Screen
%————————————————————————–
%—– Find the xyz extrema of the painting over an entire 90 deg. rotation
corners = [];
counter = 0;
for dummyAngle = 0:.1:90
counter = counter+1;
corners = cat(2,corners,Rx(beta)*Tz*Ry(dummyAngle)*Txy*C1);
end
zscreen = max(corners(3,:));
xmin = floor(min(zscreen*(corners(1,:)./corners(3,:))));
xmax = ceil(max(zscreen*(corners(1,:)./corners(3,:))));
ymin = floor(min(zscreen*(corners(2,:)./corners(3,:))));
ymax = ceil(max(zscreen*(corners(2,:)./corners(3,:))));
%—– Project painting onto a screen locaed at zscreen
x2 = zscreen*(P(1,:)./P(3,:));
y2 = zscreen*(P(2,:)./P(3,:));
cx2 = zscreen*(C(1,:)./C(3,:));
cy2 = zscreen*(C(2,:)./C(3,:));
sx2 = zscreen*(S(1,:)./S(3,:));
sy2 = zscreen*(S(2,:)./S(3,:));
clear x1 y1 z1 X1 Y1 P0 P1 P2 P3 P corners
%————————————————————————–
%—– Draw the RGB Interpolation Maps
%————————————————————————–
%—– Reshape the red, green and blue values to Nx1 arrays
red1 = double(reshape(A(:,:,1),N,1));
green1 = double(reshape(A(:,:,2),N,1));
blue1 = double(reshape(A(:,:,3),N,1));
%—– Draw the red, green and blue tri-scattered interpolation maps
disp(‘Drawing red map’)
tic
redMap = TriScatteredInterp(x2′,y2′,red1);
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
disp(‘Drawing green map’)
tic
greenMap = TriScatteredInterp(x2′,y2′,green1);
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
disp(‘Drawing blue map’)
tic
blueMap = TriScatteredInterp(x2′,y2′,blue1);
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
%—– Produce an XY mesh grid to contain the transformed image
[X2,Y2] = meshgrid(floor(min(cx2)):resolution:ceil(max(cx2)),floor(min(cy2)):resolution:ceil(max(cy2)));
clear A x2 y2 red1 green1 blue1
%—– Interpolate the RGB data onto the new mesh grid
disp(‘Interpolating red values’)
tic
red2 = uint8(round(redMap(X2,Y2)));
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
clear redMap
disp(‘Interpolating green values’)
tic
green2 = uint8(round(greenMap(X2,Y2)));
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
clear greenMap
disp(‘Interpolating blue values’)
tic
blue2 = uint8(round(blueMap(X2,Y2)));
tau = toc;
disp([‘Elapsed time: ‘,num2str(floor(tau/60)),’:’,num2str(mod(tau,60),’%02.0f’)])
clear blueMap
clear X2 Y2
%————————————————————————–
%—– Form the new RGB image
%————————————————————————–
xx = xmin:resolution:xmax; % the x-coordinates of the screen
yy = ymin:resolution:ymax; % the y-coordinates of the screen
%—– Within the screen, the rows and columns occupied by the projection
rowspan = 1+max([0,(floor(min(cy2))-ymin)/resolution]):max([0,(floor(min(cy2))-ymin)/resolution])+size(red2,1);
colspan = 1+max([0,(floor(min(cx2))-xmin)/resolution]):max([0,(floor(min(cx2))-xmin)/resolution])+size(red2,2);
rowspan = round(rowspan);
colspan = round(colspan);
%—– Construct the RGB image
B = uint8(zeros(length(yy),length(xx),3));
B(rowspan,colspan,1) = red2;
B(rowspan,colspan,2) = green2;
B(rowspan,colspan,3) = blue2;
%————————————————————————–
%—– Display and save the image
%————————————————————————–
%—– Initialize Figure Window and Axes —–
close
scrsz = get(0,’ScreenSize’);
hMainFigure = figure(‘Position’,scrsz,’Color’,[1,1,1]);
%—– Display the image
image(xx,yy,B)
axis xy
axis image
set(gca,’Visible’,’off’)
%—– Display output values
jawAngle = atand((sy2(5)-sy2(1))/(sx2(5)-sx2(1)));
disp(‘ ‘)
disp([‘R: ‘,num2str(R/Lx,’%1.3f’),’*Lx’]);
disp([‘alpha: ‘,num2str(alpha)]);
disp([‘beta: ‘,num2str(beta)]);
disp([‘skull width: ‘,num2str(max(sx2)-min(sx2))])
disp([‘skull height: ‘,num2str(max(sy2)-min(sy2))])
disp([‘jaw angle: ‘,num2str(jawAngle,’%3.4f’)])
%—– Write the image to a file
B(:,:,1) = flipud(B(:,:,1));
B(:,:,2) = flipud(B(:,:,2));
B(:,:,3) = flipud(B(:,:,3));
str = [‘~/Desktop/Y’,num2str(alpha),’.jpg’];
imwrite(B,str,’jpg’)