function dF_saveDataToOMETiff(dFLIMImageData,dFLIMXYTData,dFLIMHistoData)
% dFLIM.ome.tif format
% save a single dFLIM acquisition (contents of a .dFLIM.mat file) to
%   the OME-Tiff format.
% series 0 = dFLIM_Image (XYCTZ, with subchannels of C)
% series 1 = dFLIM_XYT   (XYtCZ)
% series 2 = dFLIM_Histo (xyCTZ, where x=t and y=1)

%% give the user a chance to set the save name
fileName='';
[fname, path]=uiputfile('*.ome.tif','Save data to OME-Tiff as...',fileName);
fname=[path fname];
fname0=fname(1:end-8); % without the extension

%% first figure out the sizes and numbers to be saved
sz=size(dFLIMImageData);

% number of Z slices
if numel(sz)<3
    szZ=1;
else
    szZ=sz(3);
end

% number of dFLIM channels
used = arrayfun(@(x) x.nPixels~=0, dFLIMImageData);
nChan = sum(used(:,1,1));
channelsAvail = find(arrayfun(@(x) ~isempty(x.img), dFLIMImageData(:,1,1))');

% size of images
typicalIndex = find(used,1,'first');
typicalDFI=dFLIMImageData(typicalIndex);
szX = typicalDFI.nPixels;
szY = typicalDFI.nLines;
nFrames = sz(2);

% histo info
nBins = numel(dFLIMHistoData(typicalIndex).Data);

% XYT available?
xytAvail = ~isempty(dFLIMXYTData) && ~isempty(dFLIMXYTData(channelsAvail(1),1).Data);
tic;

%% Create metadata
% global state
toInt = @(x) ome.xml.model.primitives.PositiveInteger(java.lang.Integer(x));
toNNInt = @(x) ome.xml.model.primitives.NonNegativeInteger(java.lang.Integer(x));
OMEXMLService = loci.formats.services.OMEXMLServiceImpl();
metadata = OMEXMLService.createOMEXMLMetadata();
metadata.createRoot();
dimensionOrderEnumHandler = ome.xml.model.enums.handlers.DimensionOrderEnumHandler();
experimentTypeEnumHandler = ome.xml.model.enums.handlers.ExperimentTypeEnumHandler();
experimentType = experimentTypeEnumHandler.getEnumeration('FluorescenceLifetime');
pixelTypeEnumHandler = ome.xml.model.enums.handlers.PixelTypeEnumHandler();

dimOrder = {'XYCTZ', 'XYTCZ', 'XYCTZ'};  % by series
ifd = 0;  % keep count of the current IFD
skipFirst=0;
%% for the moduloAlongT annotation:
coreMetadata = loci.formats.CoreMetadata();

% see the loci.formats.Modulo javadoc for details of what fields can be set
% this will be turned into a Modulo XML annotation later
coreMetadata.moduloT.start = 0;
coreMetadata.moduloT.step = 0.05;
coreMetadata.moduloT.end = 12.75;
coreMetadata.moduloT.type = loci.formats.FormatTools.LIFETIME;
coreMetadata.moduloT.unit = 'ns';

%% create the metadata for each series
for series=0:(2-skipFirst)
    sSeries=num2str(series);
    metadata.setImageID(['Image:' sSeries], series);
    metadata.setPixelsID(['Pixels:' sSeries], series);
    metadata.setPixelsBinDataBigEndian (java.lang.Boolean.TRUE, series, 0); % (java.lang.Boolean.TRUE, 0, 0);  NOT SURE ABOUT THIRD PARAMETER USAGE 
    
    % Set dimension order and experiment type for the series
    dimensionOrder = dimensionOrderEnumHandler.getEnumeration(dimOrder{series+skipFirst+1}); % GY (ip.Results.dimensionOrder);
    metadata.setPixelsDimensionOrder(dimensionOrder, series);
    
    if series==0
        metadata.setExperimentType(experimentType,series);
        metadata.setExperimentID('dFLIM_v2.1', series);
        metadata.setExperimentDescription(dF_makeheaderstr, series);
    end
    
    % Set pixels type 
    if (series+skipFirst)==2
        pixelsType = pixelTypeEnumHandler.getEnumeration('double'); 
    else
        pixelsType = pixelTypeEnumHandler.getEnumeration('uint16'); 
    end
    metadata.setPixelsType(pixelsType, series);
    
    switch series+skipFirst
        case 0  % dFLIMImages
            szXYZCT = [szX szY szZ 3*nChan nFrames];
        case 1  % XYT
            if xytAvail
                szXYZCT = [szX szY szZ nChan nBins];
            else
                szXYZCT = [0 0 0 0 0];
            end
            OMEXMLService.addModuloAlong(metadata, coreMetadata, series);
        case 2 % Histos
            szXYZCT = [nBins 1 szZ nChan nFrames];
    end % switch series (size assignment)
    
    metadata.setPixelsSizeX(toInt(szXYZCT(1)), series);
    metadata.setPixelsSizeY(toInt(szXYZCT(2)), series);
    metadata.setPixelsSizeZ(toInt(szXYZCT(3)), series);
    metadata.setPixelsSizeC(toInt(szXYZCT(4)), series);
    metadata.setPixelsSizeT(toInt(szXYZCT(5)), series);
    
    count = prod(szXYZCT(3:5));
    
    plane=0; % not sure, but think this can be used
    metadata.setTiffDataIFD(toNNInt(ifd), series, plane);
    %metadata.setUUIDFileName(filename, series, plane);
    %metadata.setUUIDValue(uuid, series, plane);
    metadata.setTiffDataPlaneCount(toNNInt(count), series, plane);
    ifd = ifd + count;
    
    frameDeltaT = 0.26;  % TODO
    if (series+skipFirst)==0
        for k=0:count-1
            metadata.setPlaneDeltaT(java.lang.Double(frameDeltaT*floor(k/3)),series,k);
        end
    elseif (series+skipFirst)==2
        for k=0:count-1
            metadata.setPlaneDeltaT(java.lang.Double(frameDeltaT*k),series,k);
        end
    end
    
    % Set channels ID and samples per pixel
    suffix={':img' ':nmg' ':tmg'};
    for i = 1:numel(channelsAvail)
        chID = num2str(channelsAvail(i),'%02i');
        if (series+skipFirst)==0
            for j=1:3
                ch = (i-1)*3+j-1;
                chname = ['dChannel:' chID suffix{j}];
                % chname = ['Channel:' chID(2) ':' num2str(j-1)];
                metadata.setChannelID(chname, series, ch);   % (['Channel:0:' num2str(i-1)], 0, i-1);
                metadata.setChannelFluor(chname, series, ch);
                metadata.setChannelSamplesPerPixel(toInt(1), series, ch);
            end
        else
            ch = (i-1);
            chname = ['dChannel:' chID ];
            metadata.setChannelID(chname, series, ch);   % (['Channel:0:' num2str(i-1)], 0, i-1);
            metadata.setChannelFluor(chname, series, ch);
            metadata.setChannelSamplesPerPixel(toInt(1), series, ch);
        end
    end % for channels avail
    
end % for series

sMetadata=OMEXMLService.getOMEXML(metadata);


%% metadata are created, now write the metadata and the data
% use the MATLAB Tiff facility
t = Tiff(fname,'w8'); % write a BigTiff

% values for first two series
tagstruct.ImageLength = szY;
tagstruct.ImageWidth = szX;
tagstruct.Photometric = Tiff.Photometric.MinIsBlack;
tagstruct.BitsPerSample = 16;
tagstruct.SamplesPerPixel = 1;
tagstruct.RowsPerStrip = szY;
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;
tagstruct.Compression = Tiff.Compression.LZW;
tagstruct.Software = 'MATLAB';


ts0 = tagstruct; % should be a copy
ts0.ImageDescription = char(sMetadata);  % only for first IFD

verbose=0;



%% SERIES 0 [dFLIM_Images]
% define the data conversion
if ~skipFirst,
index=0;
for z=1:szZ
    for frame=1:nFrames
        for ch=channelsAvail
            infostr=['DFI: z=' num2str(z) '/' num2str(szZ) '; ' ...
                'fr=' num2str(frame) '/' num2str(nFrames) '; ' ...
                'dChannel=' num2str(ch) ' [' num2str(find(ch==channelsAvail)) '/' num2str(nChan) ']; '];
            
            if verbose, disp([num2str(index) ' ' infostr]); end
            
            % get the data
            dfi = dFLIMImageData(ch,frame,z);
            
            if index==0,
                t.setTag(ts0);
            else
                t.setTag(tagstruct);
            end
            t.write(uint16(dfi.img));
            t.writeDirectory(); 
            
            t.setTag(tagstruct);
            t.write(uint16(dfi.nmg));
            t.writeDirectory(); 
            
            t.setTag(tagstruct);
            t.write(uint16(dfi.tmg));
            t.writeDirectory(); 
            index = index + 3;
        end
    end
end
end % skipFirst
%% SERIES 1 [dFLIM_XYTs]
index=0;
for z=1:szZ
    for ch=channelsAvail
        % get the xyt data structure (not shifted!!)
        dfxyt = dFLIMXYTData(ch,1,z).Data;
        dfxyt = permute(dfxyt,[2 3 1]); % this structure is ordered y/x/bin
        
        for bin=1:nBins
            infostr=['XYT: z=' num2str(z) '/' num2str(szZ) '; ' ...
                'bin=' num2str(bin) '/' num2str(nBins) '; ' ...
                'dChannel=' num2str(ch) ' [' num2str(find(ch==channelsAvail)) '/' num2str(nChan) ']; '];
            
            if verbose, disp([num2str(index) ' ' infostr]); end
            
            if skipFirst && index==0
                t.setTag(ts0);
            else
                t.setTag(tagstruct);
            end
            t.write(dfxyt(:,:,bin)); index=index+1;
            t.writeDirectory();  
        end
    end
end

%% SERIES 2 [dFLIM_Histos]
% different tag values for histos
tagstruct.ImageLength = 1;
tagstruct.ImageWidth = nBins;
tagstruct.BitsPerSample = 64;
tagstruct.RowsPerStrip = 1;
tagstruct.SampleFormat = Tiff.SampleFormat.IEEEFP;

index=0;
for z=1:szZ
    for frame=1:nFrames
        for ch=channelsAvail
            % get the histo data structure (not shifted!!)
            dfh = dFLIMHistoData(ch,frame,z);
            
            infostr=['Histo: z=' num2str(z) '/' num2str(szZ) '; ' ...
                'fr=' num2str(frame) '/' num2str(nFrames) '; ' ...
                'dChannel=' num2str(ch) ' [' num2str(find(ch==channelsAvail)) '/' num2str(nChan) ']; '];
            
            if verbose, disp([num2str(index) ' ' infostr]); end
            
            % write the "xy plane" for this histo
            t.setTag(tagstruct);
            t.write(dfh.Data(:)'); index=index+1;
            t.writeDirectory();
        end
    end
end

t.close;
toc;

end % function [dF_saveDataToOMETiff]