function geotiffwrite(varargin)
%

% Copyright 2010-2023 The MathWorks, Inc.

% Verify the number of inputs.
narginchk(3,inf);

% Parse the inputs.
[varargin{:}] = convertContainedStringsToChars(varargin{:});
[filename, A, cmap, R, Params] = parseInputs(varargin{:});

% Validate the inputs.
[filename, A, cmap, R, Params] = validateInputs(filename, A, cmap, R, Params);

% Construct TIFF tags structure based on input image characteristics and
% the tag values in the TiffTags parameter.
TiffTags = constructTiffTags(A, cmap, Params.TiffTags);

% Construct GeoTIFF tags structure based on the spatial referencing
% information from R and the coordinate reference information in
% the GeoKeyDirectoryTag parameter.
GeoTiffTags = map.geotiff.internal.constructGeoTiffTags(R, Params.GeoKeyDirectoryTag);

% Write the GeoTIFF file.
rpcTag = Params.RPCCoefficientTag;
isClassicTIFF = Params.IsClassicTIFF;
map.geotiff.internal.writeGeoTiffFile( ...
    filename, A, TiffTags, GeoTiffTags, rpcTag, isClassicTIFF);

%--------------------------------------------------------------------------

function [filename, A, cmap, R, Params] = parseInputs(filename, A, varargin)
% Parse the input arguments. Prior to calling this function, VARARGIN must
% be validated to contain at least one element (or the calling workspace
% varargin must be validated to contain at least three elements.) The
% values or expected class of the arguments are not validated in this
% function.
% 
% FILENAME and A are the first two elements of VARARGIN (from the calling
% workspace) and can be assigned directly using input variable names.
%
% PARAMS is a structure containing field names that match the parameter
% names of the interface:
%  'TiffTags', 'CoordRefSysCode', 'GeoKeyDirectoryTag',
%  'RPCCoefficientTag', 'TiffType'
% The fields are empty if the parameter is not supplied in VARARGIN.

% Get number of data parameters.
numDataArgs = internal.map.getNumberOfDataArgs(varargin{:});

% Obtain required data parameters.
switch numDataArgs
    % (FILENAME, A, 'PARAM1', ...)
    case 0
        error(message('map:geotiff:tooFewDataArgs', 'GEOTIFFWRITE', 'A', 'R'))
        
    case 1
        % (FILENAME, A, R, 'PARAM1', ...)
        R = varargin{1};
        cmap = [];
        varargin(1) = [];
              
    case 2
        % (FILENAME, X, CMAP, R, 'PARAM1' ...)       
        cmap = varargin{1};
        R = varargin{2};
        varargin(1:2) = [];
        
    otherwise
        % Include filename and A in the count.
        error(message('map:geotiff:tooManyDataArgs', 'GEOTIFFWRITE', 'A', 'CMAP', 'R', numDataArgs + 2))
end

% Parse the optional parameters and return a structure.
Params = parseParams(varargin);

%--------------------------------------------------------------------------

function Params = parseParams(inputs)
% Parse the optional input parameters and return them in the structure
% PARAMS which contain field names that match the parameter names of the
% interface:
%  'TiffTags', 'CoordRefSysCode', 'GeoKeyDirectoryTag',
%  'RPCCoefficientTag', 'TiffType'
% The fields are empty [] if the parameter is not supplied in the cell 
% array INPUTS.

map.internal.assert(mod(numel(inputs), 2) == 0, ...
    'map:geotiff:parameterCountNotEven');

% Create output structure with paramNames as the field names of the
% structure. The field values are all [].
paramNames = {'TiffTags', 'CoordRefSysCode', 'GeoKeyDirectoryTag', ...
    'RPCCoefficientTag', 'TiffType'};
Params = cell2struct(cell(size(paramNames)), paramNames, 2);

% Copy the 'PARAM1', VALUE1, ... parameters to the Params structure. Use
% validatestring to match the field names with the parameter names.
if ~isempty(inputs)
    names  = inputs(1:2:end);
    values = inputs(2:2:end);
    
    for k=1:numel(names)
        name = validatestring(names{k}, paramNames);
        Params.(name) = values{k};
    end
end

%--------------------------------------------------------------------------

function [filename, A, cmap, R, Params] = validateInputs( ...
    filename, A, cmap, R, Params)
% Validate the inputs.

try
    % Validate filename.
    filename = validateFilename(filename);
    
    % Validate image (or data grid) parameters.
    [A, cmap] = validateImage(A, cmap);
    
    % Validate Params.
    Params = validateParams(Params);
    
    % Obtain the coordinate system type: 'projected' or 'geographic'
    type = getGeoKeyCRSType(Params.GeoKeyDirectoryTag);
    
    % Validate R and if required, convert it to a referencing object.
    hasColorMap = ~isempty(cmap);
    R = validateR(R, size(A), hasColorMap, type);
    
    % Validate GeoKeyDirectoryTag raster type code ('PixelIsArea' or
    % 'PixelIsPoint') and update the GeoKeyDirectoryTag based on the raster
    % type, if needed.
    Params.GeoKeyDirectoryTag = validateRasterType(R, Params.GeoKeyDirectoryTag);
    
    % Validate image shape (rows/columns start from position).
    [A, R] = validateImageShape(A, R);
    
    % Validate the photometric interpretation value, if set, with the image.
    validatePhotometricWithImage(Params.TiffTags, A);
catch e
    throwAsCaller(e)
end

%--------------------------------------------------------------------------

function filename = validateFilename(filename)
% Validate the filename input.

validateattributes(filename, {'char','string'}, {'scalartext','nonempty'}, ...
    mfilename, 'FILENAME', 1);

% Add .tif extension to filename if necessary.
tifExt = '.tif';
[pathname, basename, ext] = fileparts(filename);
map.internal.assert(strcmpi(ext, tifExt) || isempty(ext), ...
   'map:geotiff:invalidExtension', upper( mfilename ), filename, tifExt);
filename = fullfile(pathname,[basename, tifExt]);

%--------------------------------------------------------------------------

function [A, cmap] = validateImage(A, cmap)
% Validate the input image A (or data grid) and its associated colormap
% (which is empty if the image is not indexed).

if isempty(cmap)
    % M-by-N, M-by-N-by-P (data grid or image)
    classes =  {'numeric', 'logical'};
    validateattributes(A, classes, {'nonempty'}, mfilename, 'A', 2);
else
    % M-by-N indexed image
    classes = {'uint8','uint16'};
    validateattributes(A, classes, {'nonempty'}, mfilename, 'X', 2);
    internal.map.checkcmap(cmap, mfilename, 'CMAP', 3)
    
    % Pad the colormap if needed.
    classType = class(A);
    maxNumColors = double(intmax(classType)) + 1;
    numRows = size(cmap,1);
    map.internal.assert(numRows <= maxNumColors, ...
        'map:geotiff:tooManyColormapValues', maxNumColors, numRows);
    cmap(end+1:maxNumColors,:) = 0;
end

%--------------------------------------------------------------------------

function Params = validateParams(S)
% Validate the optional parameter/value pair structure. A new structure is 
% returned containing the field names:
%    'TiffTags', 'GeoKeyDirectoryTag', 'RPCCoefficientTag', 'IsClassicTIFF'
%
% The 'CoordRefSysCode' field of S is merged into the 'GeoKeyDirectoryTag'
% field.

% Validate TiffTags field.
Params.TiffTags = validateTiffTags(S.TiffTags);

% Validate CoordRefSysCode field.
coordRefSysCode = validateCoordRefSysCode(S.CoordRefSysCode);

% Validate GeoKeyDirectoryTag field.
GeoKeyDirectoryTag = validateGeoKeyDirectoryTag(S.GeoKeyDirectoryTag);

% Merge CoordRefSysCode into the GeoKeyDirectoryTag.
Params.GeoKeyDirectoryTag = mergeCoordRefSysCode( ...
    GeoKeyDirectoryTag, coordRefSysCode);

% Validate RPCCoefficientTag, if present.
if ~isempty(S.RPCCoefficientTag)
    validateattributes(S.RPCCoefficientTag, {'map.geotiff.RPCCoefficientTag'}, ...
        {'scalar'}, mfilename, 'RPCCoefficientTag')
end
Params.RPCCoefficientTag = S.RPCCoefficientTag;

% Validate TiffType, if present.
paramName = 'TiffType';
validOptions = {'classictiff', 'bigtiff'};
defaultOption = validOptions{1};
tiffType = S.TiffType;
if ~isempty(tiffType)
    tiffType = validatestring(tiffType, validOptions, mfilename, paramName);
else
    tiffType = defaultOption;
end

 Params.IsClassicTIFF = strcmp(tiffType, defaultOption);

%--------------------------------------------------------------------------

function TiffTags = validateTiffTags(TiffTags)
% Validate the 'TiffTags' parameter value.

if isempty(TiffTags)
    TiffTags = struct();
else
    % Validate TiffTags as a scalar structure.
    validateattributes( ...
        TiffTags, {'struct'}, {'scalar'}, mfilename, 'TiffTags')
    
    % Validate the field names of TiffTags.
    TiffTags = validateTiffTagNames(TiffTags);
    
    % Validate Compression field and change string value to numeric value.
    if isfield(TiffTags, 'Compression')
        TiffTags.Compression = validateCompression(TiffTags.Compression);
    end
    
    % Validate Photometric field and change string value to numeric value.
    if isfield(TiffTags, 'Photometric')
        TiffTags.Photometric = validatePhotometric(TiffTags.Photometric);
    end
    
    % Validate RowsPerStrip with TileLength or TileWidth.
    if isfield(TiffTags, 'RowsPerStrip') && ...
            (isfield(TiffTags, 'TileLength') || ...
             isfield(TiffTags, 'TileWidth'))
         error(message('map:geotiff:invalidImageOrganization', 'TiffTags', 'RowsPerStrip', 'TileLength', 'TileWidth'));
    end
end

%--------------------------------------------------------------------------

function TiffTags = validateTiffTagNames(TiffTags)
% Validate the field names of TiffTags to be valid Tiff tag names excluding
% the reserved names. The reserved names are either:
%    1) set by the function,
%    2) GeoTIFF tag names, or
%    3) read-only tag names.

% Obtain all the tag names from the Tiff class.
tiffTagNames = Tiff.getTagNames;

% Replace ASCII with Ascii in order to be consistent in any generated error
% messages (GEOTIFFINFO uses GeoAsciiParamsTag).
tiffTagNames = strrep(tiffTagNames,'GeoASCIIParamsTag','GeoAsciiParamsTag');

% Remove 'Photometric' and replace with 'PhotometricInterpretation'
% validatestring will still allow 'Photometric' to be used with this
% substitution. Replace the field after using validatestring.
tiffTagNames = strrep(tiffTagNames,'Photometric','PhotometricInterpretation');

% Assign the names for the reserved tags.
reservedTagNames = { ...
    'BitsPerSample', 'SampleFormat', 'SamplesPerPixel', ...
    'StripByteCounts', 'StripOffsets', 'SubFileType', 'SubIFD', ...
    'TileByteCounts', 'TileOffsets', ...
    'ColorMap', 'ImageLength', 'ImageWidth', ...
    'GeoAsciiParamsTag', 'GeoDoubleParamsTag', 'GeoKeyDirectoryTag', ...
    'ModelPixelScaleTag', 'ModelTiepointTag', 'ModelTransformationTag'};

% Validate the field names, allowing only the permitted names to be set.
permittedNames = setdiff(tiffTagNames, reservedTagNames);
tagNames = fieldnames(TiffTags);
for k = 1:numel(tagNames)
    inputName = tagNames{k};
    tagName = validatestring( ...
        inputName, permittedNames, mfilename, ['TiffTags.' inputName]);
    if ~isfield(TiffTags, tagName)
        % tagName and inputName may not match even though a match is found
        % since the string may be a partial match or not match case.
        % Switch the case and/or change partial/lower case to mixed
        % case to match the field names of the Tiff class.
        TiffTags.(tagName) = TiffTags.(inputName);
        TiffTags = rmfield(TiffTags, inputName);
    end
end

% Replace 'PhotometricInterpretation' with 'Photometric' (for Tiff class)
photometricInterpretation = 'PhotometricInterpretation';
if isfield(TiffTags, photometricInterpretation)
    TiffTags.Photometric = TiffTags.(photometricInterpretation);
    TiffTags = rmfield(TiffTags, photometricInterpretation);
end

%--------------------------------------------------------------------------

function compression = validateCompression(compression)
% Validate the input COMPRESSION. If it is a string, convert the expected
% value to the value expected by the Tiff class.

% Modify 'uncompressed' to 'none' for compatibility with imfinfo.
if strcmpi(compression, 'uncompressed')
    compression = 'none';
end

permittedValues = {'LZW', 'PackBits', 'Deflate', 'None'};
compression = tiffNameToValue(compression, permittedValues, 'Compression');

%--------------------------------------------------------------------------

function photometric = validatePhotometric(photometric)
% Validate the input PHOTOMETRIC. If it is as a string, convert the
% expected string value to the value expected by the Tiff class.

permittedValues = {'MinIsBlack', 'RGB', 'Palette', 'Separated'};
photometric = tiffNameToValue(photometric, permittedValues, 'Photometric');

%--------------------------------------------------------------------------

function value = tiffNameToValue(value, permittedValues, fieldValue)
% Validate the input, VALUE, to be an element in the cell array,
% permittedValues or a numeric value that is defined by the property,
% Tiff.(fieldValue). If VALUE is a string, convert the valid string to the
% numeric value that is expected by the Tiff class (Tiff.(fieldValue))

% Obtain the values assigned by the Tiff class for the permittedValues.
values = zeros(numel(permittedValues), 1);
for k=1:numel(values)
    values(k) = Tiff.(fieldValue).(permittedValues{k});
end

% Obtain the numerical value.
valueName = ['TiffTags.' fieldValue];
if ischar(value)
    strvalue = validatestring(value, permittedValues, mfilename, valueName);
    index = strcmp(permittedValues, strvalue);
    value = values(index);
else
    validateattributes(value, ...
        {'numeric'}, {'scalar', 'finite'},  mfilename, valueName) 
    index = value == values;
    if ~any(index)
        cvalues = sprintf('%d, ', values);
        cvalues(end-1:end) = '. ';
        error(message('map:geotiff:invalidTiffTagValue', valueName, cvalues, value))
    end
end
      
%--------------------------------------------------------------------------

function GeoKeyTag = validateGeoKeyDirectoryTag(GeoKeyTag)
% Validate the 'GeoKeyDirectoryTag' parameter value.

if isempty(GeoKeyTag)
    GeoKeyTag = struct();
else
    % Validate as a scalar structure.
    validateattributes(GeoKeyTag, ...
        {'struct'}, {'scalar'}, mfilename, 'GeoKeyDirectoryTag')
    
    % Validate field names.
    GeoKeyTag = validateGeoKeyDirectoryTagNames(GeoKeyTag);
    
    % Validate data type of field values.
    validateGeoKeyDirectoryTagValues(GeoKeyTag);
    
    % Validate the coordinate reference system if defined in the
    % GeoKeyDirectoryTag value.
    GeoKeyTag = validateGeoKeyCRS(GeoKeyTag);
end

%--------------------------------------------------------------------------

function GeoKeyTag = validateGeoKeyDirectoryTagNames(GeoKeyTag)
% Validate the field names of the GeoKeyDirectoryTag.

% Construct the key ID to name map.
idToNameMap = map.geotiff.internal.constructGeoKeyDirectoryMap;
permittedNames = idToNameMap.values;

% Remove 'UserDefined'
permittedNames = strrep(permittedNames, 'UserDefined', '');

keyNames = fieldnames(GeoKeyTag);
for k = 1:numel(keyNames)
    inputName = keyNames{k};
    keyName = validatestring(inputName, permittedNames, mfilename,  ...
        ['GeoKeyDirectoryTag.' inputName]);
    if ~isfield(GeoKeyTag, keyName)
        % Switch case or partial name.
       GeoKeyTag.(keyName) = GeoKeyTag.(inputName);
       GeoKeyTag = rmfield(GeoKeyTag, inputName);
    end
end

%--------------------------------------------------------------------------

function validateGeoKeyDirectoryTagValues(GeoKeyTag)
% Validate the values in the 'GeoKeyDirectoryTag' structure, GeoKeyTag. 

% Any *Citation* key name must be a char and not empty.
names = fieldnames(GeoKeyTag);
cellIndex = regexpi(names, regexptranslate('wildcard', 'Citation'));
stringFields = names(~cellfun(@isempty, cellIndex));
numericFields = setdiff(names, stringFields);
for k=1:numel(stringFields)
    field = stringFields{k};
    validateattributes(GeoKeyTag.(field), ...
        {'char'}, {'nonempty'}, mfilename, ['GeoKeyDirectoryTag.' field])
end

% Validate other fields as numeric, non-empty, and finite.
for k=1:numel(numericFields)
    field = numericFields{k};
    validateattributes(GeoKeyTag.(field), {'numeric'},  ...
        {'nonempty', 'finite'}, mfilename, ['GeoKeyDirectoryTag.' field])
end

%--------------------------------------------------------------------------

function GeoKeyTag = validateGeoKeyCRS(GeoKeyTag)
% Validate consistency of coordinate reference system information in the
% 'GeoKeyDirectoryTag' structure, GeoKeyTag.

modelField = 'GTModelTypeGeoKey';
projField  = 'ProjectedCSTypeGeoKey';
geoField   = 'GeographicTypeGeoKey';

tfModelField = isfield(GeoKeyTag, modelField);
tfProjField  = isfield(GeoKeyTag, projField);
tfGeoField   = isfield(GeoKeyTag, geoField);

if tfProjField && tfGeoField
    % User has specified both geographic and projected systems, which is
    % valid. However, additional information must be specified.
    % Expect the GTModelTypeGeoKey to defined and to indicate projected
    % coordinates. (The coordinate reference system code is validated in the
    % next if block).
    map.internal.assert(tfModelField && GeoKeyTag.(modelField) == 1, ...
       'map:geotiff:geoAndMapProjections', projField, geoField, modelField);
end

if tfProjField
    % projected
    key = GeoKeyTag.(projField);
    geoKeyType = internal.map.getCoordRefSysCodeType(key);
    % The geoKeyType must be 'projected' or user-defined.
    map.internal.assert(any(strcmp(geoKeyType, {'projected','user-defined'})), ...
       'map:geotiff:invalidProjectedCSTypeGeoKey', [ 'GeoKeyDirectoryTag.', projField ])
    
    if tfModelField
        % GTModelTypeGeoKey = 1 indicates projected coordinates
        map.internal.assert(GeoKeyTag.(modelField) == 1, ...
           'map:geotiff:invalidProjectedModelType', [ 'GeoKeyDirectoryTag.', modelField ])
    else
        % The GTModelTypeGeoKey is not required to be set by the user but
        % is required to be set in the file. Since it is known now, 
        % set the value to indicate projected coordinates.
        GeoKeyTag.(modelField) = 1;
    end
    
elseif tfGeoField
    % geographic
    key = GeoKeyTag.(geoField);
    geoKeyType = internal.map.getCoordRefSysCodeType(key);
    % The geoKeyType must be 'geographic' or user-defined.
    map.internal.assert(any(strcmp(geoKeyType, {'geographic', 'user-defined'})), ...
      'map:geotiff:invalidGeographicTypeGeoKey', [ 'GeoKeyDirectoryTag.', geoField ])
    
    if tfModelField
        % GTModelTypeGeoKey = 2 indicates geographic coordinates
        map.internal.assert(GeoKeyTag.(modelField) == 2, ...
           'map:geotiff:invalidGeographicModelType', [ 'GeoKeyDirectoryTag.', modelField ])
    else
        % The GTModelTypeGeoKey is not required to be set by the user but
        % is required to be set in the file. Since it is known now, 
        % set the value to indicate geographic coordinates.
        GeoKeyTag.(modelField) = 2;
    end
end

%--------------------------------------------------------------------------

function coordRefSysCode = validateCoordRefSysCode(coordRefSysCode)
% Validate the 'CoordRefSysCode' parameter value. Return a string valued
% 'EPSG:' code to its equivalent double value.

if ~isempty(coordRefSysCode)   
    strName = 'CoordRefSysCode';
    if ischar(coordRefSysCode)
        value = coordRefSysCode;
        epsg = 'EPSG:';
        map.internal.assert(strncmpi(value, epsg, numel(epsg)), ...
           'map:geotiff:notEPSG', strName, epsg)
        coordRefSysCode = str2double(value(numel(epsg)+1:end));
    end
    
    % Validate class type.
    validateattributes(coordRefSysCode, ...
        {'numeric'}, {'scalar', 'finite'},  mfilename, strName)
  
    % Validate coordinate type: 'projected' or 'geographic'.
    coordRefSysCodeType = internal.map.getCoordRefSysCodeType(coordRefSysCode);
    map.internal.assert(strcmp(coordRefSysCodeType, 'projected') ...
        || strcmp(coordRefSysCodeType, 'geographic'), ...
        'map:geotiff:invalidCoordRefSysCodeType', 'CoordRefSysCode');
end

%--------------------------------------------------------------------------

function GeoKeyTag = mergeCoordRefSysCode(GeoKeyTag, coordRefSysCode)
% Merge the coordRefSysCode into the GeoKeyTag structure.

geoField  = 'GeographicTypeGeoKey';
projField = 'ProjectedCSTypeGeoKey';
tfProjField = isfield(GeoKeyTag, projField);
tfGeoField  = isfield(GeoKeyTag, geoField);

if ~tfProjField && ~tfGeoField && isempty(coordRefSysCode)
    % Set default value for coordinate reference system.
    coordRefSysCode = 4326;
    GeoKeyTag.(geoField) = coordRefSysCode;
end

% 'GeoKeyDirectoryTag' coordinate type.
geoKeyType = getGeoKeyCRSType(GeoKeyTag);

% 'CoordRefSysCode'
if ~isempty(coordRefSysCode)
    coordRefSysCodeType = internal.map.getCoordRefSysCodeType(coordRefSysCode);
    if ~isempty(geoKeyType)
        % Validate that the codes represent the same type of
        % coordinate system.
        map.internal.assert(strcmp(coordRefSysCodeType, geoKeyType), ...
           'map:geotiff:multipleCoordinateReferenceSystemTypes', 'CoordRefSysCode', geoKeyType);
    end
    % Merge the CoordRefSysCode into the GeoKeyTag structure.
    GeoKeyTag = mergeCoordRefSysCodeType( ...
        GeoKeyTag, coordRefSysCode, coordRefSysCodeType);
end
   
%--------------------------------------------------------------------------

function geoKeyType = getGeoKeyCRSType(GeoKeyTag)
% Return the coordinate system type of the input GeoKeyTag structure. The
% output string, geoKeyType will either be:
% 'geographic', 'projected', or '' if undetermined.

projField  = 'ProjectedCSTypeGeoKey';
geoField   = 'GeographicTypeGeoKey';

tfProjField  = isfield(GeoKeyTag, projField);
tfGeoField   = isfield(GeoKeyTag, geoField);
    
if tfProjField
    % projected
    key = GeoKeyTag.(projField);
    geoKeyType = internal.map.getCoordRefSysCodeType(key);
    if strcmp(geoKeyType, 'user-defined')
        geoKeyType = 'projected';
    end
    
elseif tfGeoField
    % geographic
    key = GeoKeyTag.(geoField);
    geoKeyType = internal.map.getCoordRefSysCodeType(key);
    if strcmp(geoKeyType, 'user-defined')
        geoKeyType = 'geographic';
    end
    
else
    geoKeyType = '';
end

%--------------------------------------------------------------------------

function GeoKeyTag = mergeCoordRefSysCodeType(GeoKeyTag, code, type)
% Merge the 'CoordRefSysCode', code into the GeoKeyTag structure based on
% the value of the string, type, indicating either 'projected' or
% 'geographic'. Prior to calling this function, type must be set to one
% value or the other (and is not validated in this function). The GeoKeyTag
% fields are set to indicate either a projected or geographic coordinate
% system and the projection code is set.

modelField = 'GTModelTypeGeoKey';
tfModelField = isfield(GeoKeyTag, modelField);

if strcmp(type, 'projected')
    if tfModelField
        % GTModelTypeGeoKey = 1 indicates projected coordinates
        map.internal.assert(GeoKeyTag.(modelField) == 1, ...
          'map:geotiff:invalidProjectedModelType', [ 'GeoKeyDirectoryTag.', modelField ])
    else
        GeoKeyTag.(modelField) = 1;
    end
    GeoKeyTag.ProjectedCSTypeGeoKey = code;
else
    if tfModelField
        % GTModelTypeGeoKey = 2 indicates geographic coordinates
        map.internal.assert(GeoKeyTag.(modelField) == 2, ...
           'map:geotiff:invalidGeographicModelType', [ 'GeoKeyDirectoryTag.', modelField ])
    else
        GeoKeyTag.GTModelTypeGeoKey = 2;
    end
    GeoKeyTag.GeographicTypeGeoKey = code;
end

%--------------------------------------------------------------------------

function R = validateR(R, rasterSize, hasColorMap, coordRefSysType)
% Validate R and return a referencing object.

if hasColorMap
    % (FILENAME, A, CMAP, R)
    arg_pos = 4;
else
    % (FILENAME, A, R)
    arg_pos = 3;
end

if strcmp(coordRefSysType, 'geographic')
    if isobject(R) && isscalar(R) && isprop(R,'CoordinateSystemType') ...
            && strcmp(R.CoordinateSystemType,'planar')
        error(message('map:geotiff:expectedGeoRasterReference', 'R', ...
            class(R), 'CoordRefSysCode', 'GeoKeyDirectoryTag'));
    else
        angle_units = 'degrees';

        if isequal(size(R),[3 2]) || numel(R) == 3
            warning(message("map:geotiff:refvecmat","GEOTIFFWRITE"))
        end
        R = internal.map.convertToGeoRasterRef(R, rasterSize, angle_units, ...
            upper(mfilename), 'R', arg_pos);
    end   
elseif strcmp(coordRefSysType, 'projected')
    if isequal(size(R),[3 2]) % Check for refmat only since refvec will produce error
        warning(message("map:geotiff:refvecmat","GEOTIFFWRITE"))
    end
    R = convertToMapRasterRef(R, rasterSize, upper(mfilename), 'R', arg_pos);
end
    
%--------------------------------------------------------------------------

function R = convertToMapRasterRef( ...
    R, rasterSize, func_name, var_name, arg_pos)
% Validate and convert R to a map raster reference object.

if numel(R) == 3
    error(message('map:geotiff:refvecNotAllowed','R','map raster reference'));
end

R = refmatToMapRasterReference(R, rasterSize, func_name, var_name, arg_pos);

map.internal.assert(strcmp(R.TransformationType, 'rectilinear'), ...
   'map:geotiff:skewOrRotation', func_name, func_name)

%--------------------------------------------------------------------------

function GeoKeyTag = validateRasterType(R, GeoKeyTag)
% Validate GeoKeyDirectoryTag raster type code to be consistent with the
% referencing object.

if strcmp(R.RasterInterpretation, 'cells')
    % RasterPixelIsArea
    rasterType = 1;
else
    % RasterPixelIsPoint
    rasterType = 2;
end

field = 'GTRasterTypeGeoKey';
if isfield(GeoKeyTag, field)
    map.internal.assert(rasterType == GeoKeyTag.(field), ...
       'map:geotiff:inconsistentRasterType', [ 'GeoKeyDirectoryTag.', field ],  ...
          rasterType, R.RasterInterpretation, GeoKeyTag.(field));
else
    GeoKeyTag.(field) = rasterType;
end

%--------------------------------------------------------------------------

function [A, R] = validateImageShape(A, R)
% Validate the image shape and flip the rows or columns if needed to be
% consistent with the GeoTIFF specification.

% Assuming rectilinear.
if strcmp(R.ColumnsStartFrom, 'south')
    R.ColumnsStartFrom = 'north';
    A = A(end:-1:1,:,:);
end

if strcmp(R.RowsStartFrom, 'east')
    R.RowsStartFrom = 'west';
    A = A(:,end:-1:1,:);
end

%--------------------------------------------------------------------------

function validatePhotometricWithImage(TiffTags, A)
% Validate the photometric interpretation value with the image.

if isfield(TiffTags, 'Photometric') ...
        && TiffTags.Photometric == Tiff.Photometric.RGB
    map.internal.assert(ndims(A) == 3 && size(A,3) == 3, ...
       'map:geotiff:imageSizeNotRGB', 'PhotometricInterpretation')
end
        

%--------------------------------------------------------------------------

function TiffTags = constructTiffTags(A, cmap, UserTags)
% Construct TIFF tags from R and the GeoKeyTag structure.

% Create TiffTags structure based on characteristics of input image.
TiffTags = assignTiffTags(A, cmap);

% Modify TiffTags based on user inputs.
TiffTags = updateTagsFromInputs(TiffTags, UserTags);

%--------------------------------------------------------------------------

function TiffTags = assignTiffTags(A, cmap)
% Construct the TiffTags structure based on image characteristics.

% Create structure of TIFF tags.
TiffTags = struct( ...
    'Photometric', [], ...
    'Compression', [], ...
    'ImageLength', [], ...
    'ImageWidth',  [], ...
    'SamplesPerPixel', [], ...
    'BitsPerSample',   [], ...
    'SampleFormat', [], ...
    'PlanarConfiguration', [], ...
    'Software', []); 

% Assign the tag value for Photometric.
if ~isempty(cmap)
    photometric = Tiff.Photometric.Palette;
    TiffTags.ColorMap = cmap;
else
    numBands = size(A, 3);
    isRGB = ndims(A) == 3 && numBands == 3 ...
        && any(strcmp(class(A), {'uint8', 'uint16', 'uint32'}));
    if isRGB
        photometric = Tiff.Photometric.RGB;
    else
        photometric = Tiff.Photometric.MinIsBlack;
        if numBands > 1 && ~islogical(A)
            % Signed datatypes require MinIsBlack Photometric
            % Interpolation. SamplesPerPixel is numBands, but MinIsBlack
            % corresponds to only one sample. Create a [1xnumBands-1] array
            % of EmptySamples to represent the remaining samples.
            extraSamples = Tiff.ExtraSamples.Unspecified;
            n = max(1, numBands-1);
            extraSamples(1:n) = extraSamples(1);
            TiffTags.ExtraSamples = extraSamples;
        end
    end
end
TiffTags.Photometric = photometric;

% Assign the default tag value for Compression.
TiffTags.Compression = Tiff.Compression.PackBits;

% Assign the tag values for the image size.
TiffTags.ImageLength = size(A, 1);
TiffTags.ImageWidth  = size(A, 2);

% Assign the tag value for SamplesPerPixel.
map.internal.assert(ndims(A) <= 3, 'map:geotiff:invalidNumberOfDimensions');
if ~islogical(A)
    TiffTags.SamplesPerPixel = size(A, 3);
else
    if ~ismatrix(A)
        % Only one logical channel is supported by Tiff.
        % Use same message as imwrite.
        error(message('MATLAB:imagesci:writetif:unsupportedData'))
    end
    TiffTags.SamplesPerPixel = 1;
end

% Assign the tag values for BitsPerSample and SampleFormat.
classType = class(A);
TiffTags.BitsPerSample = assignBitsPerSample(classType);
TiffTags.SampleFormat  = assignSampleFormat(classType);

% Assign the tag value for PlanarConfiguration.
TiffTags.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky;

% Assign file organization tags.
TiffTags = assignFileOrganizationTags(TiffTags, A);

% Assign the tag value for Software.
v1 = matlabRelease;
v2 = ver('map');
software = ['MATLAB' ' ' char(v1.Release) ', ' v2.Name, ' ' v2.Version];
TiffTags.Software = software;

%--------------------------------------------------------------------------

function bitsPerSample = assignBitsPerSample(classType)
% Assign the value for BitsPerSample based on class type.

switch classType
    case {'logical'}       
        bitsPerSample = 1;
        
    case {'uint8', 'int8'}
        bitsPerSample = 8;
        
    case {'uint16', 'int16'}
        bitsPerSample = 16;
        
    case {'uint32', 'int32', 'single'}
        bitsPerSample = 32;

    case {'double'}
        bitsPerSample = 64;

    otherwise
        error(message('map:geotiff:unknownClass', classType));
end

%--------------------------------------------------------------------------

function sampleFormat = assignSampleFormat(classType)
% Assign the value for SampleFormat based on class type.

switch classType
    case {'logical', 'uint8', 'uint16', 'uint32'}       
        sampleFormat = Tiff.SampleFormat.UInt;
                
    case {'int8', 'int16', 'int32'}
        sampleFormat = Tiff.SampleFormat.Int;
                    
    case {'single', 'double'}
        sampleFormat = Tiff.SampleFormat.IEEEFP;

    otherwise
        error(message('map:geotiff:unknownClass', classType));
end

%--------------------------------------------------------------------------

function TiffTags = assignFileOrganizationTags(TiffTags, A)
% Assign the file organization tags, TileLength and TileWidth or
% RowsPerStrip.

% Set the tile size if the image is at least 160-by-160. Note the minimum
% tile dimension is 16-by-16. TileLength and TileWidth must be multiples of
% 16. Create tiles such that the maximum number of tiles is 100 (which can
% be obtained from the Tiff class method, numberOfTiles).
if size(A,1) >= 160 && size(A,2) >= 160
    maxTileWidth  = size(A,2)/10;
    maxTileLength = size(A,1)/10;
    tileWidth  = fix(maxTileWidth/16)*16;
    tileLength = fix(maxTileLength/16)*16;
    
    if tileWidth ~= maxTileWidth
        % Increase tileWidth to next multiple of 16.
        tileWidth = tileWidth + 16;
    end
    if tileLength ~= maxTileLength
        % Increase tileLength to next multiple of 16.
        tileLength = tileLength + 16;
    end
   
    TiffTags.TileWidth  = tileWidth;
    TiffTags.TileLength = tileLength;
else
    TiffTags.RowsPerStrip = 1;
end

%--------------------------------------------------------------------------

function TiffTags = updateTagsFromInputs(TiffTags, UserTags)
% Copy the fields from UserTags to TiffTags.

userNames = fieldnames(UserTags);

% Remove TileLength and TileWidth if RowsPerStrip has been specified.
tileNames = {'TileLength', 'TileWidth'};
if any(strcmpi('RowsPerStrip', userNames)) ...
        && all(isfield(TiffTags, tileNames))
   TiffTags = rmfield(TiffTags, tileNames);
end

% Remove Software if the user requested the empty string.
field = 'Software';
index = strcmpi(field, userNames);
if any(index) && isempty(UserTags.(field))
    TiffTags = rmfield(TiffTags, field);
    UserTags = rmfield(UserTags, field);
    userNames(index) = [];
end

% Copy or create the fields from UserTags to TiffTags.
for k=1:numel(userNames)
    field = userNames{k};
    TiffTags.(field) = UserTags.(field);
end
