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{:});

[filenameRSet, outputDir] = buildFilename(sourceImage, varargin{:});

adapter = buildAdapter(sourceImage);

adapter.isSupported(sourceImage);
verifyOutputDir(outputDir);

tileSize = 512;
info = getBasefileInfo(sourceImage, adapter);
tileDims = getTileDims(tileSize, info.SamplesPerPixel);

[numSpanningTilesRows, numSpanningTilesCols, maxLevel] = computeRSetDetails(info, tileSize);
details = computeStaticImageDetails(sourceImage, adapter);

[ids, zero_tile] = setupFile(filenameRSet, tileDims, details);

baseGroupID = createLevelGroup(ids, 0);
groupID = createLevelGroup(ids, 1);

% Create the waitbar.
numTiles = computeNumberOfNonVirtualTiles(info, tileSize, maxLevel);
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));

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

virtualTileTable = false(numSpanningTilesRows, numSpanningTilesCols);

for row = 0:numSpanningTilesRows-1
    
    for col = 0:numSpanningTilesCols-1
        
        if (waitState.isCancelled)
            cancelAction(ids, groupID, baseGroupID, 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) <= info.Height) && (imageCols(1) <= info.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 = padTile(data, tileSize);
        end
        
        if (~isVirtualTile)
            
            waitState.update();
            
            % Write the L0 data.
            datasetName = sprintf('r%d_c%d', 2*row, 2*col);
            writeTile(data(1:tileSize, 1:tileSize, :), ids, details, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col);
            writeTile(data((tileSize+1):(2*tileSize), 1:tileSize, :), ids, details, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row, 2*col + 1);
            writeTile(data(1:tileSize, (tileSize+1):(2*tileSize), :), ids, details, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col + 1);
            writeTile(data((tileSize+1):(2*tileSize), (tileSize+1):(2*tileSize), :), ids, details, baseGroupID, datasetName)
            
            % Write the L1 data.
            datasetName = sprintf('r%d_c%d', row, col);
            data = resizeImage(data, details.cmap);
            writeTile(data, ids, details, groupID, datasetName)
            
        else
            
            % Write the L0 data.
            datasetName = sprintf('r%d_c%d', 2*row, 2*col);
            writeTileRef(ids, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col);
            writeTileRef(ids, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row, 2*col + 1);
            writeTileRef(ids, baseGroupID, datasetName)
            
            datasetName = sprintf('r%d_c%d', 2*row + 1, 2*col + 1);
            writeTileRef(ids, baseGroupID, datasetName)
            
            % Write the L1 data.
            datasetName = sprintf('r%d_c%d', row, col);
            writeTileRef(ids, groupID, datasetName)
            
        end
        
    end
end

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

% Build the other levels.
for level = 2:maxLevel
    
    groupID = createLevelGroup(ids, level);
    
    maxInboundRow = ceil(info.Height / (tileSize * 2^level)) - 1;
    maxInboundCol = ceil(info.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;
                writeTileRef(ids, groupID, datasetName)
                continue
                
            end
            
            waitState.update();
            
            % Upper Left
            data = readTile(ids, level-1, 2*row, 2*col, virtualTileTable);
            new_tile(1:tileSize/2,1:tileSize/2,:) = resizeImage(data, details.cmap);
            
            % Upper Right
            data = readTile(ids, level-1, 2*row, 2*col+1, virtualTileTable);
            new_tile(1:tileSize/2,tileSize/2+1:end,:) = resizeImage(data, details.cmap);
            
            % Lower Left
            data = readTile(ids, level-1, 2*row+1, 2*col, virtualTileTable);
            new_tile(tileSize/2+1:end,1:tileSize/2,:) = resizeImage(data, details.cmap);
            
            % Lower Right
            data = readTile(ids, level-1, 2*row+1, 2*col+1, virtualTileTable);
            new_tile(tileSize/2+1:end,tileSize/2+1:end,:) = resizeImage(data, details.cmap);
            
            % write new tile
            writeTile(new_tile, ids, details, groupID, datasetName)
            
        end
    end
    
    H5G.close(groupID);
    
    virtualTileTable = virtualTileTable2;
    
end

% clean up wait bar
clear cleanup_waitState;

writeColormap(ids.file, details.cmap);
writeMetadata(ids.file, sourceImage, info, tileSize, maxLevel);
writeVersion(ids.file);

H5G.create(ids.rootID, '/BaseFileMetadata', 65536);

closeIdentifiers(ids);

writeAllBasefileMetadata(filenameRSet, info.meta);

filename = 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 = ~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 info = getBasefileInfo(sourceImage, adapter)

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

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

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


function tileDims = getTileDims(tileSize, samplesPerPixel)

if (samplesPerPixel > 1)
    
    tileDims = [tileSize, tileSize, samplesPerPixel];
    
    if (samplesPerPixel ~= 3)
        warning(message('images:rsetwrite:samplesPerPixel', samplesPerPixel));
    end
    
else
    tileDims = [tileSize, tileSize];
end


function h5Sid = createDataspace(tileDims)
h5Sid = H5S.create_simple(numel(tileDims), fliplr(tileDims), fliplr(tileDims));


function h5Fid = createFile(filenameRSet)

h5Fid = H5F.create(filenameRSet, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');


function h5Pid = createDataspaceProps(tileDims, details)

h5Pid = H5P.create('H5P_DATASET_CREATE');

chunkSide = min(128, tileDims(1));

if (numel(tileDims) == 2)
    chunkDims = [chunkSide, chunkSide];
else
    % Reverse the tile dimensions to compensate for majority change.
    chunkDims = [tileDims(end:-1:3), chunkSide, chunkSide];
end

H5P.set_chunk(h5Pid, chunkDims);

switch (class(details.mlType))
    case {'uint8', 'int8'}
        % Don't shuffle.
    otherwise
        H5P.set_shuffle(h5Pid);
end

H5P.set_deflate(h5Pid, 6);


function [ids, zero_tile] = setupFile(filenameRSet, tileDims, details)

% Identifiers common to all datasets.
ids.file     = createFile(filenameRSet);
ids.dspace   = createDataspace(tileDims);
ids.propList = createDataspaceProps(tileDims, details);

% Build the root for all layer data.
ids.rootID = H5G.create(ids.file, '/RSetData', 512);

% Create a "virtual" dataset that padding tiles can use via reference
% instead of bloating the file by writing the same tile.
ids.virtualGroup = H5G.create(ids.rootID, 'Virtual', 32);
zero_tile = zeros(tileDims, details.mlType);
padDsetID = H5D.create(ids.virtualGroup, 'padding', details.hdfType, ids.dspace, 'H5P_DEFAULT');
H5D.write(padDsetID, 'H5ML_DEFAULT', ids.dspace, ids.dspace, 'H5P_DEFAULT', zero_tile);
H5D.close(padDsetID);

ids.dtypeRef = H5T.copy('H5T_STD_REF_OBJ');
ids.ref = H5R.create(ids.virtualGroup, 'padding', 'H5R_OBJECT', -1);
ids.dspaceRef = H5S.create('H5S_SCALAR');


function closeIdentifiers(ids)

% Close in opposite order of opening.
if (~isequal(ids.propList, 'H5P_DEFAULT'))
    H5P.close(ids.propList);
end

H5T.close(ids.dtypeRef);
H5S.close(ids.dspace);
H5S.close(ids.dspaceRef);
H5G.close(ids.virtualGroup);
H5G.close(ids.rootID);
H5F.close(ids.file);


function data = readTile(ids, level, row, column, virtualTileTable)

dsetName = sprintf('L%d/r%d_c%d', level, row, column);
dset = H5D.open(ids.rootID, dsetName);

%if (H5ML.compare_values(H5T.get_class(H5D.get_type(dset)), ...
%                        H5ML.get_constant_value('H5T_REFERENCE')))
if (virtualTileTable(row+1, column+1))
    
    % Read a "virtual" padding tile that is stored elsewhere in the file.
    dsetDeref = H5R.dereference(dset, 'H5R_OBJECT', ids.ref);
    data = H5D.read(dsetDeref, 'H5ML_DEFAULT', ids.dspace, ids.dspace, 'H5P_DEFAULT');
    H5D.close(dsetDeref);
    
else
    
    data = H5D.read(dset, 'H5ML_DEFAULT', ids.dspace, ids.dspace, 'H5P_DEFAULT');
    
end

H5D.close(dset);


function writeColormap(h5Fid, cmap)

groupID = H5G.create(h5Fid, 'Colormap', 32);

if (~isempty(cmap))
    dspaceID = H5S.create_simple(2, fliplr(size(cmap)), []);
    dsetID  = H5D.create(groupID, 'map', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
    H5D.write(dsetID, 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', cmap);
    H5D.close(dsetID);
    H5S.close(dspaceID);
end

H5G.close(groupID);


function writeMetadata(h5Fid, sourceImage, imageMetadata, tileSize, maxLevel)

% Attach these values to the "/Metadata" group.
groupID = H5G.create(h5Fid, 'Metadata', 32);
dspaceID = H5S.create('H5S_SCALAR');

% The side length of each tile.
attrID = H5A.create(groupID, 'TileSize', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', tileSize);
H5A.close(attrID);

% The dimensions of the full-resolution image.
attrID = H5A.create(groupID, 'FullImageHeight', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', imageMetadata.Height);
H5A.close(attrID);

attrID = H5A.create(groupID, 'FullImageWidth', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', imageMetadata.Width);
H5A.close(attrID);

attrID = H5A.create(groupID, 'MaxLevel', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', maxLevel);
H5A.close(attrID);

% The name of the full resolution image file.
if isa(sourceImage,'ImageAdapter')
    source_name = [class(sourceImage),' object'];
else
    source_name = sourceImage;
end
stringType = H5T.copy('H5T_C_S1');
H5T.set_size(stringType, length(source_name));
H5T.set_strpad(stringType, 'H5T_STR_NULLTERM');
attrID = H5A.create(groupID, 'OriginalBaseFile', stringType, dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, stringType, source_name);
H5T.close(stringType);
H5A.close(attrID);

attrID = H5A.create(groupID, 'FileCreationTime', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', now);
H5A.close(attrID);

attrID = H5A.create(groupID, 'BaseFileCreationTime', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', imageMetadata.Datenum);
H5A.close(attrID);

H5S.close(dspaceID);
H5G.close(groupID);


function writeVersion(fileID)

groupID = H5G.create(fileID, 'FormatInfo', 32);
dspaceID = H5S.create('H5S_SCALAR');

stringType = H5T.copy('H5T_C_S1');
H5T.set_size(stringType, length('RSET'));
attrID = H5A.create(groupID, 'FileType', stringType, dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, stringType, 'RSET');
H5A.close(attrID);
H5T.close(stringType);

formatDescription = 'Multi-resolution, large imagery format';
stringType = H5T.copy('H5T_C_S1');
H5T.set_size(stringType, length(formatDescription));
attrID = H5A.create(groupID, 'Description', stringType, dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, stringType, formatDescription);
H5A.close(attrID);
H5T.close(stringType);

attrID = H5A.create(groupID, 'Version', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', iptui.RSet.currentRSetVersion);
H5A.close(attrID);

attrID = H5A.create(groupID, 'BackwardVersion', 'H5T_NATIVE_DOUBLE', dspaceID, 'H5P_DEFAULT');
H5A.write(attrID, 'H5ML_DEFAULT', iptui.RSet.backwardRSetVersion);
H5A.close(attrID);

H5S.close(dspaceID);
H5G.close(groupID);


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));


function out = resizeImage(in, cmap)

if isempty(cmap)
    out = imresize(in, 0.5);
else
    out = imresize(in, cmap, 0.5, 'Colormap', 'original');
end


function [numSpanningTilesRows, numSpanningTilesCols, maxLevel] = computeRSetDetails(info, tileSize)

numSpanningRows = 2 ^ ceil(log2(info.Height));
numSpanningTilesRows = numSpanningRows / tileSize;

numSpanningCols = 2 ^ ceil(log2(info.Width));
numSpanningTilesCols = numSpanningCols / tileSize;

min_dim = min(info.Width, info.Height);
base_dim = 2 ^ ceil(log2(min_dim));
base_tiles = base_dim / tileSize;
maxLevel = max(ceil(log2(base_tiles)), 1);

% Handle small images.
if ((info.Height < 3) || (info.Width < 3))
    
    error(message('images:rsetwrite:imageTooSmall'))
    
elseif ((numSpanningTilesRows < 2) || (numSpanningTilesCols < 2))
    
    warning(message('images:rsetwrite:smallImage'))
    
    numSpanningTilesRows = max(numSpanningTilesRows, 2);
    numSpanningTilesCols = max(numSpanningTilesCols, 2);
    
end


function [filenameRSet, outputDir] = buildFilename(sourceImage, varargin)

if numel(varargin) > 0
    % the image adapter syntax requires an explicitly specified output
    % filename, so this code path will always be followed for image adapter
    % inputs.
    filenameRSet = varargin{1};
    outputDir = fileparts(filenameRSet);
    
else
    % if a filename has been specified, we can generate an output filename
    % based on the input file.
    outputDir = pwd;
    [~, fname] = fileparts(sourceImage);
    filenameRSet = fullfile(outputDir, [fname '.rset']);
    
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 groupID = createLevelGroup(ids, level)

groupID = H5G.create(ids.rootID, sprintf('L%d', level), 65536);


function writeTile(data, ids, details, groupID, datasetName)

dsetID = H5D.create(groupID, datasetName, details.hdfType, ids.dspace, ids.propList);
H5D.write(dsetID, 'H5ML_DEFAULT', ids.dspace, ids.dspace, 'H5P_DEFAULT', data);
H5D.close(dsetID);


function writeTileRef(ids, groupID, datasetName)

dsetID = H5D.create(groupID, datasetName, ids.dtypeRef, ids.dspaceRef, 'H5P_DEFAULT');
H5D.write(dsetID, 'H5T_STD_REF_OBJ', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT', ids.ref)
H5D.close(dsetID);


function verifyOutputDir(outputDir)

[success, attributes] = fileattrib(outputDir);

if (~success)
    
    error(message('images:rsetwrite:nonexistentOutputDir', outputDir))
    
elseif (~attributes.UserWrite)
    
    error(message('images:rsetwrite:outputDirPerms', outputDir))
    
end


function data = padTile(data, tileSize)

padAmount = [2*tileSize - size(data,1), 2*tileSize - size(data,2), 0];
data = padarray(data, padAmount, 'post', 'replicate');


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 writeAllBasefileMetadata(filenameRSet, info)

% Write metadata found in the top-level of the source file's info
% structure, skipping empty values.  Do not recurse down the hierarchy.

f = fieldnames(info);

for i = 1:numel(f)
    
    attr = info(1).(f{i});
    
    if (~isempty(attr) && ...
            (ischar(attr) || isnumeric(attr)))
        
        % Attributes should be smaller than 16k bytes.
        attrDetails = whos('attr');
        if (attrDetails.bytes <= 16384)
            h5writeatt(filenameRSet, '/BaseFileMetadata', f{i}, attr);
        end
    end
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 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);



