% estimateSensorSpectralResponse: estimate sensor spectral sensitivity 
%
% response = estimateSensorSpectralResponse(input_spectrum,rgb,'method')
%
% input_spectrum: spectrum of input light
% rgb: imaging device response to input light
% method: 
%    Implemented: 
%
%       'gsvd': Generalized SVD method [Default]
%               Reference: Bob Dyas, Robust color sensor response characterization, CIC 2000
%
%       'quad': Quadratic programming
%    
%    Non-implemented yet:
%
%       'weiner': Wiener estimation method 
%               Reference: W.K.Pratt and C.E.Mancill, Spectral estimation
%               techniques for the spectral calibration of a color image
%               scanner, Applied Optics, 1976 p73-75
%
%       'rdpi': Rank-deficient pseudo inverse method 
%
%       'pocs': Projections Onto Convex Sets (POCS)  
%               Reference: Gaurav Sharma and H.J. Trussell, Characterization of scanner 
%                       sensitivity, CIC 1993
%
% General reference: 
%       Paul Hubel et. al. A Comparison of Methods of Sensor Spectral
%       Sensitivity Estimation, 2nd Color Imaging Conference, 1994
%
% Written by  : Feng Xiao
% Last Updated: 03-25-03

function response = estimateSensorSpectralResponse(input_spectrum,rgb,method)

if (nargin < 2 )
   error('estimateSensorSpectralResponse requires at least two parameters');
end

if ~exist('method','var')
    method='gsvd';
end

M = length(rgb);
channels = prod(size(rgb))/M;
N = prod(size(input_spectrum))/M;

if size(input_spectrum,1) ~= M
    input_spectrum = input_spectrum';
end

if size(rgb,1) ~= M
    rgb = rgb';
end

switch lower(method)
    case 'gsvd'
        [U,S,V] = svd(input_spectrum);
        weight = 0.1;
        S_recip = 1./diag(S);
        % figure;
        for r=1:length(U)
            response{r} = V(:,1:r)*diag(S_recip(1:r))*U(:,1:r)'*rgb;
            tmp = input_spectrum*response{r}-rgb;
            err(r) = sum(abs(tmp(:)));
            derive = diff(diff(response{r}));
            smooth(r) = sum(abs(derive(:)));
            % plot(response{r}); title(num2str(r)); drawnow;
            % pause(0.5);
        end
        cost = smooth/max(smooth)*weight + err/max(err);
        [tmp,index]=min(cost);
        response = response{index};
        % figure; semilogy(cost,'*-'); grid on; 
        % figure; plot(response{index}); title(num2str(index));
        
    case 'quad'
        for color=1:channels
            H = input_spectrum'*input_spectrum;
            f = -input_spectrum'*rgb(:,color);
            A = -eye(N); b=zeros(N,1);
            [tmp,lamda] = qp(H,f,A,b); %quadprog(H,f); %A,b);
            lamda %= 1; %rgb(:,color)'*input_spectrum*tmp/(tmp'*H*tmp);
            response(:,color) = tmp*lamda; %*;
            %response(:,color) = tmp*sum(rgb(:,color))/sum(input_spectrum*tmp);
        end    

    otherwise 
        error(['Method ' method '  is not supported']);
end
