function filename = rsetwrite(sourceImage, varargin)
%RSETWRITE Create reduced-resolution dataset (R-Set) from image file.
%   RSETFILE = rsetwrite(FILENAME) creates a reduced resolution dataset (an
%   R-Set) from the specified TIFF or NITF image file.  The generated .rset
%   file is written to the current working directory with a name based on
%   the input FILENAME.  (For example, if FILENAME is
%   'VeryLargeImage.tiff', then RSETFILE will be 'VeryLargeImage.rset'.)
%
%   RSETFILE = rsetwrite(FILENAME, OUTPUT_FILENAME) creates an R-Set from
%   the specified TIFF image file, using OUTPUT_FILENAME as the name of the
%   new file.
%
%   RSETFILE = rsetwrite(ADAPTER, OUTPUT_FILENAME) creates an R-Set from
%   the specified image adapter object, ADAPTER.  Image adapters are
%   user-defined classes that provide RSETWRITE a common API for reading a
%   particular image file format.  See the documentation for ImageAdapter
%   for more details.
%
%   When OUTPUT_FILENAME is specified, RSETFILE and OUTPUT_FILENAME will be
%   the same.
%
%   Because R-Set creation can be time-consuming, a progress bar shows the
%   status of the operation.  If the progress bar is cancelled, processing
%   stops, no file is written, and the RSETFILE variable will be empty.
%
%   Note
%   ----
%   An R-Set file contains multiple views of the same large image which has
%   been divided into spatial tiles and which has been resampled into
%   different resolution levels.  The R-Set file also contains a compressed
%   copy of the original full-resolution data.
%
%   The .rset file generated by this function is an HDF5 file with the
%   multi-resolution data organized for quick access.  Use the IMTOOL
%   function to explore the imagery within the .rset file.
%
%   This function is designed to make working with very large imagery
%   easier.  While it is possible to create an R-Set from an image where
%   the image dimensions are smaller than the size of a single R-Set
%   tile, the resulting .rset file will likely be larger and take longer
%   to load than the original file.
%
%   If a TIFF or NITF file contains multiple images, only the first one
%   is used.
%
%   When given a NITF file, RSETWRITE requires that the file have a
%   version of at least 2.0, contain an image, and be uncompressed.
%   Files with more than three bands are not supported.  Floating point
%   data is not supported, and images with more than one data band must
%   contain unsigned integer data.
%
%   Example
%   -------
%   % Visualize a very large image using a reduced resolution dataset.
%   % (Replace 'MyReallyBigImage.tif' with the filename of your large image.)
%   big_file = 'MyReallyBigImage.tif';
%   rset_file = rsetwrite(big_file);
%   imtool(rset_file)
%
%   % Create .rset files for every TIFF file in a directory containing
%   % very large images.  Put them into a temporary directory.
%   d = dir('*.tif*');
%   image_dir = pwd;
%   cd(tempdir)
%   for p = 1:numel(d)
%       big_file = fullfile(image_dir, d(p).name);
%       rsetwrite(big_file);
%   end
%
%   See also IMREAD, IMTOOL.

%   Copyright 2008-2017 The MathWorks, Inc.

sourceImage = matlab.images.internal.stringToChar(sourceImage);
varargin = matlab.images.internal.stringToChar(varargin);
validateInputs(sourceImage,varargin{:});

adapter = buildAdapter(sourceImage);
adapter.isSupported(sourceImage);

metadata = getBasefileInfo(sourceImage, adapter);
metadata.details = computeStaticImageDetails(sourceImage, adapter);

obj = images.internal.rset.RSetWriter(sourceImage,varargin{:});

obj.setupFile(metadata);



baseGroupID = obj.createLevelGroup(0);
groupID = obj.createLevelGroup(1);

% Create the waitbar.
numTiles = computeNumberOfNonVirtualTiles(obj.metadata, obj.tileSize, obj.maxLevel);
[waitState, cleanup_waitState] = createWaitbar(numTiles);

% Build the L0 and L1 datasets.
numSpanningTilesRows = obj.numSpanningTilesRows / 2;
numSpanningTilesCols = obj.numSpanningTilesCols / 2;

tileSize = obj.tileSize;
ids = obj.ids;
details = obj.metadata.details;
zero_tile = zeros(obj.tileDims, details.mlType);

virtualTileTable = false(numSpanningTilesRows, numSpanningTilesCols);

for row = 0:numSpanningTilesRows-1
    
    for col = 0:numSpanningTilesCols-1
        
        if (waitState.isCancelled)
            cancelAction(obj.ids, groupID, baseGroupID, obj.filenameRSet)
            filename = [];
            return
        end
        
        isVirtualTile = false;
        
        % Read the full resolution data.
        imageRows = [row*2*tileSize + 1, (row+1)*2*tileSize];
        imageCols = [col*2*tileSize + 1, (col+1)*2*tileSize];
        
        if ((imageRows(1) <= obj.metadata.Height) && (imageCols(1) <= obj.metadata.Width))
            data = adapter.read(sourceImage, 'PixelRegion', {imageRows imageCols});
        else
            data = zero_tile;
            isVirtualTile = true;
            virtualTileTable(row+1, col+1) = true;
        end
        
        % Pad the data if the tile is on the edge and a little short.
        if ((size(data,1) ~= 2*tileSize) || (size(data,2) ~= 2*tileSize))
            data = obj.padTile(data);
        end
        
        if (~isVirtualTile)
            
            waitState.update();
            
            % Write the L0 data.
            datasetName = sprintf('r%d_c%d', 2*row, 2*col);
            obj.writeTile(data(1:tileSize, 1:tileSize, :), baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col);
            obj.writeTile(data((tileSize+1):(2*tileSize), 1:tileSize, :), baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row, 2*col + 1);
            obj.writeTile(data(1:tileSize, (tileSize+1):(2*tileSize), :), baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col + 1);
            obj.writeTile(data((tileSize+1):(2*tileSize), (tileSize+1):(2*tileSize), :), baseGroupID, datasetName)
            
            % Write the L1 data.
            datasetName = sprintf('r%d_c%d', row, col);
            data = obj.resizeImage(data);
            obj.writeTile(data, groupID, datasetName)
            
        else
            
            % Write the L0 data.
            datasetName = sprintf('r%d_c%d', 2*row, 2*col);
            
            obj.writeTileRef(baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col);
            obj.writeTileRef(baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row, 2*col + 1);
            obj.writeTileRef(baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col + 1);
            obj.writeTileRef(baseGroupID, datasetName)
            
            % Write the L1 data.
            datasetName = sprintf('r%d_c%d', row, col);
            obj.writeTileRef(groupID, datasetName)
            
        end
        
    end
end

H5G.close(groupID);
H5G.close(baseGroupID);

% Build the other levels.
for level = 2:obj.maxLevel
    
    groupID = obj.createLevelGroup(level);
    
    maxInboundRow = ceil(obj.metadata.Height / (tileSize * 2^level)) - 1;
    maxInboundCol = ceil(obj.metadata.Width  / (tileSize * 2^level)) - 1;
    
    numSpanningTilesRows = numSpanningTilesRows/2;
    numSpanningTilesCols = numSpanningTilesCols/2;
    
    virtualTileTable2 = false(numSpanningTilesRows, numSpanningTilesCols);
    
    for row = 0:numSpanningTilesRows-1
        for col = 0:numSpanningTilesCols-1
            
            if (waitState.isCancelled)
                cancelAction(ids, groupID, [], filenameRSet)
                filename = [];
                return
            end
            
            datasetName = sprintf('r%d_c%d', row, col);
            new_tile = zero_tile;
            
            if ((row > maxInboundRow) || (col > maxInboundCol))
                
                virtualTileTable2(row+1, col+1) = true;
                obj.writeTileRef(groupID, datasetName)
                continue
                
            end
            
            waitState.update();
            
            % Upper Left
            data = obj.readTile(level-1, 2*row, 2*col, virtualTileTable);
            new_tile(1:tileSize/2,1:tileSize/2,:) = obj.resizeImage(data);
            
            % Upper Right
            data = obj.readTile(level-1, 2*row, 2*col+1, virtualTileTable);
            new_tile(1:tileSize/2,tileSize/2+1:end,:) = obj.resizeImage(data);
            
            % Lower Left
            data = obj.readTile(level-1, 2*row+1, 2*col, virtualTileTable);
            new_tile(tileSize/2+1:end,1:tileSize/2,:) = obj.resizeImage(data);
            
            % Lower Right
            data = obj.readTile(level-1, 2*row+1, 2*col+1, virtualTileTable);
            new_tile(tileSize/2+1:end,tileSize/2+1:end,:) = obj.resizeImage(data);
            
            % write new tile
            obj.writeTile(new_tile, groupID, datasetName)
            
        end
    end
    
    H5G.close(groupID);
    
    virtualTileTable = virtualTileTable2;
    
end

% clean up wait bar
clear cleanup_waitState;

obj.writeColormap();
obj.writeAllMetadata();

filename = obj.filenameRSet;


function validateInputs(sourceImage,varargin)

    valid_file = ischar(sourceImage) && isequal(exist(sourceImage,'file'),2);
    valid_adapter = isa(sourceImage,'ImageAdapter');

    if valid_file

        % Validate file type
        [~, ~, ext] = fileparts(sourceImage);
        tiff_extensions = {'.tiff', '.tif'};
        %is_tiff = validatestring(ext,{'.tiff', '.tif'});
        is_tiff = ~isempty(strmatch(lower(ext),lower(tiff_extensions)));
        is_nitf = isnitf(sourceImage);

        if ~(is_tiff || is_nitf)
            error(message('images:rsetwrite:unsupportedInputFormat'));
        end

    elseif valid_adapter

        % An output file is required when the sourceImage is an ImageAdapter
        if numel(varargin) < 1
            error(message('images:rsetwrite:mustSpecifyOutputFile'));
        end

        % Validate Colormap, must be m-by-3 matrix of numeric
        map = sourceImage.Colormap;
        if ~isempty(map)
            if ~isequal(ndims(map),2) || ~isequal(size(map,2),3) || ~isnumeric(map)
                error(message('images:rsetwrite:invalidColormap'));
            end
        end

    else
        error(message('images:rsetwrite:invalidSourceImage'))
    end
    
function adapter = buildAdapter(sourceImage)
    if isa(sourceImage,'ImageAdapter')
        % these image adapter helper functions mirror the imread/imfinfo
        % syntaxes that the rest of RSETWRITE expects.

        adapter.format             = 'ImageAdapter';
        adapter.read               = @readFromIA;
        adapter.readImageAndCmap   = @readFromIAWithCMap;
        adapter.getInfo            = @infoFromIA;
        adapter.getHeight          = @(meta) meta.Height;
        adapter.getWidth           = @(meta) meta.Width;
        adapter.getSamplesPerPixel = @(meta) meta.SamplesPerPixel;
        adapter.isSupported        = @(filename) true;

    elseif isnitf(sourceImage)

        adapter.format = 'NITF';
        adapter.read = @nitfread;
        adapter.getInfo = @nitfinfo;
        adapter.getHeight = @(meta) ...
            meta.ImageSubheaderMetadata.ImageSubheader001.NumberOfSignificantRowsInImage;
        adapter.getWidth = @(meta) ...
            meta.ImageSubheaderMetadata.ImageSubheader001.NumberOfSignificantColumnsInImage;
        adapter.getSamplesPerPixel = @(meta) computeNumberOfNitfBands(meta);

        % NITFREAD only returns one output argument, so always assign an
        % empty colormap.
        adapter.readImageAndCmap = @(varargin) deal(nitfread(varargin{:}), []);

        % Only certain NITF files are supported.
        adapter.isSupported = @isNitfSupported;

    else

        [~, ~, ext] = fileparts(sourceImage);
        switch (lower(ext))
            case {'.tiff', '.tif'}

                adapter.format             = 'TIFF';
                adapter.read               = @imread;
                adapter.readImageAndCmap   = @imread;
                adapter.getInfo            = @imfinfo;
                adapter.getHeight          = @(meta) meta.Height;
                adapter.getWidth           = @(meta) meta.Width;
                adapter.getSamplesPerPixel = @(meta) meta.SamplesPerPixel;
                adapter.isSupported        = @(filename) true;
        end

    end
    
function metadata = getBasefileInfo(sourceImage, adapter)

    meta = adapter.getInfo(sourceImage);
    % If there are multiple images in the file, use the first.
    metadata.meta = meta(1);

    metadata.SamplesPerPixel = adapter.getSamplesPerPixel(metadata.meta);
    metadata.Height = adapter.getHeight(metadata.meta);
    metadata.Width = adapter.getWidth(metadata.meta);

    % for image adapter objects, set the .Datenum field to now.
    if isa(sourceImage,'ImageAdapter')
        metadata.Datenum = now;
    else
        % Both NITF and TIFF have a .Filename field on the info struct
        d = dir(metadata.meta.Filename);
        metadata.Datenum = d.datenum;
    end
        
function details = computeStaticImageDetails(sourceImage, adapter)

    % For the sake of efficiency, precompute some "constant" values that we
    % would otherwise only discover every time through the loop.
    [data, details.cmap] = adapter.readImageAndCmap(sourceImage, ...
        'PixelRegion', {[1 1], [1 1]});

    dtypeMap = makeDatatypeMap;
    details.mlType = class(data);

    % It's nontrivial to get 1-bit data back out of the HDF .rset file.
    % Disable the use of RSETWRITE with logical data until we have requests
    % for it.
    if (isequal(details.mlType, 'logical'))

        error(message('images:rsetwrite:logicalImage'))

    end

    details.hdfType = dtypeMap(details.mlType);
    
function [waitState, cleanup_waitState] =  createWaitbar(numTiles)

if images.internal.isFigureAvailable()
    waitState = iptui.cancellableWaitbar('Building R-Set:',...
        'Building R-Set with %d tiles', numTiles);
else
    waitState = iptui.textWaitUpdater('Building R-Set with %d tiles.',...
        'Completed %d of %d tiles.', numTiles);
end
cleanup_waitState = onCleanup(@() destroy(waitState));

function cancelAction(ids, groupID, baseGroupID, filenameRSet)

H5G.close(groupID);

if (~isempty(baseGroupID))
    H5G.close(baseGroupID);
end

closeIdentifiers(ids);

delete(filenameRSet);

function numTiles = computeNumberOfNonVirtualTiles(info, tileSize, maxLevel)

h = info.Height;
w = info.Width;

numTiles = 0;
for level = 1:maxLevel
    
    % The image dimensions at each level are half the size of the preceding
    % level.  Divide by tileSize and round up to compute the number of
    % tiles for a given level.
    numTiles = numTiles + (ceil(h / (tileSize * 2^level)) * ...
        ceil(w / (tileSize * 2^level)));
end

function dtypeMap = makeDatatypeMap

typeMapping = {'int8',    'H5T_NATIVE_INT8'
    'uint8',   'H5T_NATIVE_UINT8'
    'int16',   'H5T_NATIVE_INT16'
    'uint16',  'H5T_NATIVE_UINT16'
    'int32',   'H5T_NATIVE_INT32'
    'uint32',  'H5T_NATIVE_UINT32'
    'single',  'H5T_NATIVE_FLOAT'
    'double',  'H5T_NATIVE_DOUBLE'};
dtypeMap = containers.Map(typeMapping(:,1), typeMapping(:,2));


%......... Adapter functions ......... %
function tf = isNitfSupported(filename)

[tf, eid, msg] = iptui.isNitfSupported(filename);

if (~tf && ~isequal(eid, 'images:isNitfSupported:nitfNumberOfBands'))
    eid = sprintf(eid, 'rsetwrite');
    error(eid, msg)
end


function numBands = computeNumberOfNitfBands(meta)

if (isfield(meta.ImageSubheaderMetadata.ImageSubheader001, 'NumberOfMultiSpectralBands'))
    numBands = meta.ImageSubheaderMetadata.ImageSubheader001.NumberOfMultiSpectralBands;
else
    numBands = meta.ImageSubheaderMetadata.ImageSubheader001.NumberOfBands;
end


function data = readFromIA(sourceImage,~,region)

last_row = min(sourceImage.ImageSize(1),region{1}(2));
last_col = min(sourceImage.ImageSize(2),region{2}(2));
row_start = region{1}(1);
row_size = last_row - region{1}(1) + 1;
col_start = region{2}(1);
col_size = last_col - region{2}(1) + 1;
data = sourceImage.readRegion([row_start col_start],[row_size col_size]);


function [data, map] = readFromIAWithCMap(sourceImage,~,region)
map  = sourceImage.Colormap;
data = readFromIA(sourceImage,'PixelRegion',region);


function info = infoFromIA(sourceImage)
info.SamplesPerPixel = getSamplesfromIA(sourceImage);
info.Height = sourceImage.ImageSize(1);
info.Width = sourceImage.ImageSize(2);


function spp = getSamplesfromIA(sourceImage)
sample = sourceImage.readRegion([1 1],[1 1]);
spp = numel(sample);


function [rows, cols] = getL1SpanningTiles(tileSize, xlims, ylims)
% Compute the spanning tiles, clamping to in-bound values.
col_min = floor(xlims(1) / tileSize);
col_max = floor(xlims(2) / tileSize);
row_min = floor(ylims(1) / tileSize);
row_max = floor(ylims(2) / tileSize);

cols = (col_min:col_max)/2;
rows = (row_min:row_max)/2;



