function [Output] = MF(Image,Library)

% Matched Filter
% Silas J. Leavesley, PhD
% Сòòò½APP
% Last Updated: 2/7/2017

% To visualize output data:
% figure(1)
% for i = 1:length(Output.data(1,1,:))
%   subplot(1,length(Output.data(1,1,:)),i)
%   imagesc(Output.data(:,:,i))
% end
% colormap(gray) % or any other colormap, as desired

tic

Image=double(Image);

[O P] = size(Library);
[L M N]=size(Image);
Image_Reshaped=(reshape(Image,(L*M),N))'; %L M are x y N=# wavelengths
K=L*M; % Total # pixels in image

% Initialize library variable - library gets shifted by 1 endmember each
% time around the loop
Library_shifted = Library;

% Run matched filter for each endmember
for i = 1:P
    % Specify desired spectral signature (d)
    d = Library_shifted(:,1);
    
    % Remaining endmembers are undersired signatures (U)
    U = Library_shifted(:,2:P); 

    % Define the pseudo-inverse of U
    U_pseudo = (inv(transpose(U)*U))*transpose(U);
    
    % Define the rejection operator
    I = eye(N);
    rejection_operator = (I-(U*U_pseudo));

    % Define the filter function
    q = transpose(d)*rejection_operator;

    % Multiply image by filter function (pixel-by-pixel operation, in
    % vectorized form)
    for j = 1:K
        MF(j,i) = q*Image_Reshaped(:,j);
    end
    % Shift library to next endmember
    Library_shifted = circshift(Library_shifted,[0 -1]);
end
Output.data=reshape(MF,L,M,P);
toc