% textmapping
% 11/8/00
% Dorian Miller


function mainTM(angleHor, angleVer, checker)
% usage: mainTM();
% Computes the text map over an arbitray region.

cla
if (checker ==1)
	% construct a checkerboard texture
	tw = 64;
	th = 64;
	SW = 6;
	[ty, tx] = meshgrid(1:th, 1:tw);
	bitmap = mod(fix((tx-1)/SW),2) ~= mod(fix((ty-1)/SW),2);
	% This is the background color
	%bitmap(1,1)=0;
else
   % Create a texture map from a grayscale picture -- picture has to be in 
   % working directory
	bitmap = imread('well2.jpg');
	bitmap = double(bitmap);
	[th, tw] = size(bitmap);
end

% the area to texture map
uv = [0 0 1; 0 th   1; tw   th   1; tw   0 1]';
xy = [0 0 1; 0 2*th 1; 2*tw 2*th 1; 2*tw 0 1]';
xy = transform(xy, [-round(tw/2) -round(th/2) 0]);
% Perform projection of camera view !!!!
xy = projection(xy,angleVer,angleHor,200)
xy = transform(xy, [150 250 0]);
xy = xy.*(1./repmat(xy(3,:),3,1))

%xy = [-0.5 0 1; 0 1 1; 1 1 1; 1.5 0 1]';

% find the translation between bitmap and skewed projection. !!!!
Map = texturemap(uv,xy);

% Map all pointes of the image plane to those of the 
% texture map.
W = 256;
H = 256;
% image plane coordinates 
[imx,imy] = meshgrid(1:H,1:W);
imcolor = zeros(1,H*W);
imCoord = [imx(:)';imy(:)';ones(1,H*W)]';
 size(imCoord)
bitmapCoord = inv(Map)*imCoord';

%scale results -- all z-coords = 1
bitmapCoord = bitmapCoord.*(1./repmat(bitmapCoord(3,:),3,1));
bitDDDD  =max(bitmapCoord(3,:))
bitDDDD  =min(bitmapCoord(3,:))
%find the valid indcies that lye in the bitmap
xCoord = bitmapCoord(1,:);
yCoord = bitmapCoord(2,:);
validColor = ones(1,H*W);
minY = min(yCoord)
minX = min(xCoord)

xCoordFL = floor(xCoord);
yCoordFL = floor(yCoord);
xCoordCL = ceil(xCoord);
yCoordCL = ceil(yCoord);
xCoordFR = xCoord - floor(xCoord);
yCoordFR = yCoord - floor(yCoord);


validX = find(0>=xCoordFL | xCoordFL>tw); % WIDTH BITMAP???
xCoordFL(validX) = 1;
yCoordFL(validX) = 1;
validColor(validX)= 0; 
validY = find(0>=yCoordFL | yCoordFL>th); % HEIGHT BITMAP???
xCoordFL(validY) = 1;
yCoordFL(validY) = 1;
validColor(validY)= 0;

validX = find(0>=xCoordCL | xCoordCL>tw); % WIDTH BITMAP???
xCoordCL(validX) = 1;
yCoordCL(validX) = 1;
validColor(validX)= 0; 
validY = find(0>=yCoordCL | yCoordCL>th); % HEIGHT BITMAP???
xCoordCL(validY) = 1;
yCoordCL(validY) = 1;
validColor(validY)= 0;

minY = min(yCoord)
minX = min(xCoord)

% bitmap index (bI) defines the four corners 
bIFLFL = sub2Ind([th,tw],yCoordFL,xCoordFL);
bICLCL = sub2Ind([th,tw],yCoordCL,xCoordCL);
bIFLCL = sub2Ind([th,tw],yCoordCL,xCoordFL);
bICLFL = sub2Ind([th,tw],yCoordFL,xCoordCL);

% resphae bitmap to index into
bitmap1D = reshape(bitmap,1,tw*th);

% valid color index (ci) is the same length as bitmapIndex
validCI = find(validColor==1);

% x00(1-u)(1-v)+x01(u)(1-v)+x10(1-u)v+x11uv
% u=x v=y
imcolor(validCI) = ...
   (1-xCoordFR(bIFLFL(validCI))).*(1-yCoordFR(bIFLFL(validCI))).*bitmap1D(bIFLFL(validCI))...
   +(xCoordFR(bIFLFL(validCI))).*(1-yCoordFR(bIFLFL(validCI))).*bitmap1D(bICLFL(validCI))...
   +(1-xCoordFR(bIFLFL(validCI))).*(yCoordFR(bIFLFL(validCI))).*bitmap1D(bIFLCL(validCI))...
   +(xCoordFR(bIFLFL(validCI))).*(yCoordFR(bIFLFL(validCI))).*bitmap1D(bICLCL(validCI));

% No filtering
%imcolor(validCI) = bitmap1D(bIFLFL(validCI));



% place the color in the final image.
% reshape in final step
finalImage = zeros(1,W*H);
finalImage(validCI) = imcolor(validCI);
finalImage = reshape(finalImage,H,W);

colormap('gray');
imagesc(finalImage);


function newM = rotation(m, angleHor,angleVer)

newM = rotationHor(m,    angleHor);
newM = rotationVer(newM, angleVer);

% rotation about y,z-plane pointing straight up
% angle given in degrees
% rotation in horizontal plane
function newM = rotationHor(m, angle)
	m = [m; ones(1,size(m,2))];
	cosA = cos(angle*pi/180);
	sinA = sin(angle*pi/180);
	
	rotate = [cosA  sinA 0 0;...
	          -sinA cosA 0 0;...
	         0      0    1 0;...
	       	0      0    0 1];
	newM = rotate*m;
	newM= newM(1:3,:);	
   
function newM = rotationVer(m, angle)
m = [m; ones(1,size(m,2))];
cosA = cos(angle*pi/180);
sinA = sin(angle*pi/180);

rotate = [1 0    0     0;...
          0 cosA -sinA 0;...
          0 sinA cosA  0;...
          0 0    0     1];
newM = rotate*m;
newM= newM(1:3,:);   
   
% transform points by shift vector
function newM = transform(m, shift)

m = [m; ones(1,size(m,2))];
trans = [1 0 0 shift(1);...
         0 1 0 shift(2);...
         0 0 1 shift(3);...
         0 0 0 1]
   
newM = trans*m; 
newM= newM(1:3,:);


function newPoints = projection(data,angleHor,angleVer,fl)
% The one true way
% The Axis of the viewing plane
X = [1  0 0 ];
Y = [0  -1 0];
% Vector from the eye to enter of viewing plane
eye = [0 0 -fl];
P = [X;Y;eye]';
P = rotation(P, angleHor, angleVer);
% Light vector points to light (opposite direction of the eye.)
% Light is fixed at user's head to always ensure a good view

%Projection in Image Plane = inv(P)*(Data-COP)
dataXYZ = [data(1,:);data(2,:);data(3,:)];

% Image plane coordinates, camera 400 pixels from object
copV = rotation([0 0 300]',angleHor, angleVer); %P(:,3)./(-1)
numPoints = size(data,2)
copRep = repmat(copV,1,numPoints);

% NEWLY CALCULATED IMAGE POINTS
% datacop also helps determine visibility
datacop = dataXYZ-copRep;
newPoints = inv(P)*(datacop);