clear all

parent_directory = uigetdir('C:\','Select a directory with spectral z-images.  Images should be named with "imagename_z-point#_channel#.tif');

directory_search_name=strcat(parent_directory,'\*.tif');
file_names=ls(directory_search_name);

unmixed_directory=strcat(parent_directory,'\Unmixed');
mkdir(unmixed_directory);

prompt={'Enter the imagename (without zslice# and channel#):','Enter the number of wavelength bands:','Enter the number of ZSlices:','Is this a FRET measurement (y/n)?','End-Member Number of Donor (in Library):','End-Member Number of Acceptor (in Library):'};
name='Wavelength and z-stack Information';
numlines=1;
def={'imagename','31','76'};
answer=inputdlg(prompt,name,numlines);

image_name=answer{1};
num_wavelengths=str2num(answer{2});
num_zpoints=str2num(answer{3});
isFRET=answer{4};
Donor_Number=str2num(answer{5});
Acceptor_Number=str2num(answer{6});

Z_Point_Digits=numel(num2str(num_zpoints));
Z_Point_Format=strcat('%0',num2str(Z_Point_Digits),'.0f');

[Wavelength_name,Wavelength_path]=uigetfile('*.mat','Select the wavelength information file',parent_directory);
load(strcat(Wavelength_path,Wavelength_name));

[Library_name,Library_path]=uigetfile('*.mat','Select the spectra library file',Wavelength_path);
load(strcat(Library_path,Library_name));
[O P] = size(Library);

temp_image=imread(strcat(parent_directory,'\',image_name,'z',num2str(1,Z_Point_Format),'c1.tif'));
[M N]=size(temp_image);

% Pre-allocate size for large image matrices
Spectral_Stack=zeros(M,N,num_wavelengths);
Unmixed_Images=zeros(M,N,P,num_zpoints);
RMS_Error=zeros(M,N,num_zpoints);

Waitbar_Handle = waitbar(0,'Please wait...','Name','Zpoint #');

for i=1:num_zpoints
    waitbar(i/num_zpoints,Waitbar_Handle,sprintf('Converting Zpoint %5.0f out of %5.0f',i,num_zpoints))
    for j=1:num_wavelengths
        temp_image=imread(strcat(parent_directory,'\',image_name,'z',num2str(i,Z_Point_Format),'c',num2str(j),'.tif'));
        %temp_image=uint16(temp_image); %convert to 16-bit image
        Spectral_Stack(:,:,j)=temp_image;
    end
    
    [A E E_Percent]=linear_unmix_vectorized(Wavelength,Library,Spectral_Stack);
%    Unmixed_Images(:,:,:,i)=A;
%    RMS_Error(:,:,i)=E;
    Mean_RMS_Error(i)=mean(mean(E));
    Mean_Percent_RMS_Error(i)=mean(mean(E_Percent));
   
    % Write unmixed images as separate tiffs
    for k=1:P
%        temp_unmixed=uint16(A(:,:,k));
        imwrite(uint16(A(:,:,k)),strcat(unmixed_directory,'\',image_name,'_',num2str(i,Z_Point_Format),'_Unmixed_',EndMember_Name{k},'.tif'),'tif');
    end
    imwrite(uint16(E),strcat(unmixed_directory,'\',image_name,'_',num2str(i,Z_Point_Format),'_RMS_Error','.tif'),'tif');
    imwrite(uint16(E_Percent*2^16),strcat(unmixed_directory,'\',image_name,'_',num2str(i,Z_Point_Format),'_RMS_Percent_Error','.tif'),'tif');

    % Perform FRET calculation and save FRET image if FRET is indicated
    % Note that scaled FRET efficiency is scaled by 16-bits (2^16) because
    % efficiencies between 0-1 would be rounded to either 0 or 1
    if isFRET=='y'
        FRET_Efficiency(:,:,i)=1-(A(:,:,Donor_Number)./(A(:,:,Donor_Number)+A(:,:,Acceptor_Number)));
        FRET_Efficiency(isnan(FRET_Efficiency))=0;
        imwrite(uint16(FRET_Efficiency(:,:,i)*(2^16)),strcat(unmixed_directory,'\',image_name,'_',num2str(i,Z_Point_Format),'_FRET','.tif'),'tif');
        Donor_Plus_Acceptor_Signal(:,:,i)=A(:,:,Donor_Number)+A(:,:,Acceptor_Number);
        imwrite(uint16(Donor_Plus_Acceptor_Signal(:,:,i)),strcat(unmixed_directory,'\',image_name,'_',num2str(i,Z_Point_Format),'_DonorAcceptorTotal','.tif'),'tif');
    end

end

close(Waitbar_Handle);


% multi_tiff_filename=strcat(parent_directory,'\','multi_tif_',file_names(1,1:(length(file_names(1,:))-7)),'.tif');
% band_sequential_filename=strcat(parent_directory,'\','multi_tif_',file_names(1,1:(length(file_names(1,:))-7)),'.bsq');
% band_sequential_header_filename=strcat(parent_directory,'\','multi_tif_',file_names(1,1:(length(file_names(1,:))-7)),'.hdr');
% 
% for i=1:length(file_names(:,1))
%     temp_image=imread(strcat(parent_directory,'/',file_names(i,:)));
%     wavelength_str=file_names(i,(length(file_names(1,:))-6):(length(file_names(1,:))-4));
%     wavelength(i)=str2num(wavelength_str);
% %    temp_image=uint16(temp_image); %convert to 16-bit image
%     multi_tiff(:,:,i)=temp_image; 
% end
% 
% 
% % Save band-sequential tiff file
% multibandwrite(multi_tiff,band_sequential_filename,'bsq');
% 
% % Generate ENVI header file
% [l,m,n]=size(multi_tiff);
% redband=22;
% greenband=15;
% blueband=8;
% fid = fopen(band_sequential_header_filename, 'wt');
% count=fprintf(fid,'ENVI\ndescription = {\n  File Imported into ENVI.}\nsamples = ');
% count=fprintf(fid,'%d\n',l);
% count=fprintf(fid,'lines   = ');
% count=fprintf(fid,'%d\n',m);
% count=fprintf(fid,'bands   = ');
% count=fprintf(fid,'%d\n',n);
% count=fprintf(fid,'header offset = 0\n');
% count=fprintf(fid,'file type = ENVI Standard\n');
% count=fprintf(fid,'data type = 12\n');
% count=fprintf(fid,'interleave = bsq\n');
% count=fprintf(fid,'sensor type = Unknown\n');
% count=fprintf(fid,'byte order = 0\n');
% count=fprintf(fid,'default bands = {');
% count=fprintf(fid,'%d,',redband);
% count=fprintf(fid,'%d,',greenband);
% count=fprintf(fid,'%d}\n',blueband);
% count=fprintf(fid,'wavelength units = nm\n');
% count=fprintf(fid,'wavelength = {\n');
% count=fprintf(fid,' %d',wavelength(1));
% for i=1:(n-1)
%     count=fprintf(fid,',\n %d',wavelength(i+1));
% end
% count=fprintf(fid,'}\n');
% fclose(fid);