function imgProj = projectCD(imgColor, frame, desc, imgColorDB, frameDB, descDB, corners)
% imgProj = projectCD(imgColor, frame, desc, imgColorDB, frameDB, descDB, corners)
% Projects a single CD cover into a rectified square.

% -------------------------------------------------------------------------
% Project from query space into database space
% -------------------------------------------------------------------------

% Extract feature correspondences
corners = corners(2:end);
X = corners(1:2:7); % query side (x,y) pairs
Y = corners(2:2:8);
XY = corners(1:8);
xy = corners(9:12); % database side (x,y) pairs
x = [xy(1) xy(3) xy(1) xy(3)];
y = [xy(2) xy(2) xy(4) xy(4)];

% Calculate projective model
c = XY.';
T = zeros(8,8);
for n = 1:4
    T(2*n-1,:) = [x(n) y(n) 1 0 0 0 -X(n)*x(n) -X(n)*y(n)];
    T(2*n,:) = [0 0 0 x(n) y(n) 1 -Y(n)*x(n) -Y(n)*y(n)];
end
a = inv(T) * c;
disp('Projective model:'); disp(a);

% Map query into database space
disp('Projecting query image ...'); tic;
imgColor = double(imgColor);
imgColorDB = double(imgColorDB);
minx = min(x); maxx = max(x);
miny = min(y); maxy = max(y);
[heightDB, widthDB, channels] = size(imgColorDB);
[heightQ, widthQ, channels] = size(imgColor);
imgProj = 128 * ones(heightDB, widthDB, channels);
for x = minx:maxx
    for y = miny:maxy
        
        denom = a(7)*x + a(8)*y + 1;
        X = (a(1)*x + a(2)*y + a(3)) / denom;
        Y = (a(4)*x + a(5)*y + a(6)) / denom;
        X = round(X);
        Y = round(Y);
        if (X > 0) && (X <= widthQ) && (Y > 0) && (Y <= heightQ)
           imgProj(y,x,1:channels) = imgColor(Y,X,1:channels); 
        end
        
    end % y
end % x
toc;
% subplot(2,2,3);
% warning off; imshow(uint8(imgProj)); title('Projected'); warning on;
pause(1);

% -------------------------------------------------------------------------
% Image registration to minimize shift errors
% -------------------------------------------------------------------------

% Register query projection with database image
disp('Registering query and database images ...'); tic;
minError = inf;
for deltax = -10:10
    for deltay = -10:10
        
        % Shift query
        imgProjNew = 128 * ones(heightDB, widthDB, channels);
        rangex = max(1,1+deltax):min(widthDB,widthDB+deltax);
        rangey = max(1,1+deltay):min(heightDB,heightDB+deltay);
        rangex2 = max(1,1-deltax):min(widthDB,widthDB-deltax);
        rangey2 = max(1,1-deltay):min(heightDB,heightDB-deltay);
        imgProjNew(rangey, rangex, 1:channels) = imgProj(rangey2, rangex2, 1:channels);        
%         deltax
%         deltay
%         subplot(2,2,4);
%         warning off; imshow(uint8(imgProjNew)); title('Shifted'); warning on; 
%         pause(5);
        
        % Calculate error
        error = sum( abs(imgProjNew(:) - imgColorDB(:)) );
        if error < minError
           minError = error;
           bestShift = [deltax deltay];
           imgProjBestShift = imgProjNew;
        end
        
    end % deltay
end % deltax
if minError < inf
   imgProj = imgProjBestShift; 
end
toc;
disp('Best shift (dx,dy):'); disp(bestShift);
pause(1);

warning off; imgProj = uint8(imgProj); warning on;