function [Unmixed_Image RMS_Error_Image RMS_Percent_Error_Image] = linear_unmix(Wavelength,Library,Image)
% clear all
[L M N]=size(Image);
[O P]=size(Library);

Image_Reshaped=(reshape(Image,(L*M),N))';
Library_Reshaped=Library';

%Spectral_Stack_Reshaped_Back=reshape(Spectral_Stack_Reshaped,M,N,O);

Unmixed_Array=lsqnonnegvect(Library,Image_Reshaped);

Unmixed_Image=reshape(Unmixed_Array',L,M,P);

% Calculate Fit Spectra
Library_z=permute(Library,[3 2 1]);
Fit_Spectral_Image=zeros(L,M,N);
Library_matrix_temp=zeros(L,M,P);
coeff_matrix_temp=zeros(L,M,P);
for i=1:P
    Library_matrix_temp=repmat(Library_z(1,i,:),[L,M]);
    coeff_matrix_temp=repmat(Unmixed_Image(:,:,i),[1,1,O]);
    Fit_Spectral_Image=Fit_Spectral_Image+(coeff_matrix_temp.*Library_matrix_temp);
end

% Calculate Error
Residual_Spectral_Image=Image-Fit_Spectral_Image;
RMS_Error_Image=((sum((Residual_Spectral_Image.^2),3))./P).^0.5;
RMS_Signal_Image=((sum((Image.^2),3))./P).^0.5;
RMS_Percent_Error_Image=RMS_Error_Image./RMS_Signal_Image;