% Copyright (c) 2015 The Johns Hopkins University/
% Applied Physics Laboratory. (JHU/APL)
% 
% This material was developed by JHU/APL under U.S. Government contract.
% The U.S. Government has Unlimited Rights in this material per
% DFARS 252.227-7013/7014, as appropriate.
% 
% For any other permission, please contact the Legal Office at JHU/APL.
%
function [ outImage, outRef ] = imwarpGeoTiff(varargin)

    % imwarpGeoTiff( filenameFrom, filenameTo, fillValue )
    % imwarpGeoTiff( imageFrom, infoFrom, infoTo, fillValue )

    fillValue = [];
    if ischar( varargin{1} ) && ischar( varargin{2} )
        filenameFrom =  varargin{1};
        filenameTo   =  varargin{2};
        
        if nargin >= 3
            fillValue    = varargin{3};
        end
        
        imageFrom = geotiffread(filenameFrom);
        infoFrom  = geotiffinfo(filenameFrom);
        infoTo    = geotiffinfo(filenameTo);
    
    elseif  (isnumeric( varargin{1} ) ||  islogical( varargin{1} )  )&&  ...
            isstruct( varargin{2} ) && isstruct( varargin{3} )  
        
        imageFrom = varargin{1};
        infoFrom  = varargin{2};
        infoTo    = varargin{3};
        
        if nargin >= 4
            fillValue    = varargin{4};
        end
        
        
        
    end
        
    refMatFrom = infoFrom.RefMatrix;   
    refMatTo   = infoTo.RefMatrix;  % row/col -> xy
    
    
    
    % Check if nothing to do
    if false %all(all( abs(refMatFrom - refMatTo) < 1e-12 ))  && ...
          ( infoFrom.Height == infoTo.Height ) && ...
          ( infoFrom.Width  == infoTo.Width  )
        
        outImage = imageFrom;
        if ~isempty(fillValue)
            outImage(imageFrom == fillValue) = nan;
        end
        return;
    else

        outSize = [infoTo.Height infoTo.Width];
      
        pts =  [1 1; 1 size(imageFrom,1);  size(imageFrom,2)  1; size(imageFrom,2) size(imageFrom,1)];
        [x, y] = pix2map(refMatFrom,pts(:, [2 1]));
        [r, c] = map2pix(refMatTo, x, y);
        tform = fitgeotrans(pts, [c r], 'Affine');
        
        % This should match variables from [~,ref] = geotiffread(file)
        XIntrinsicLimits = [0 outSize(2)] + 0.5;
        YIntrinsicLimits = [0 outSize(1)] + 0.5;
        
        outRef = imref2d(outSize, XIntrinsicLimits, YIntrinsicLimits);  
        [outImage,outRef] = imwarp(imageFrom, tform, 'OutputView', outRef, 'FillValues', nan);
  
        outImage = double(outImage);
        % Mask out any values influenced by fill values
        if ~isempty(fillValue)
            tmp = single(imageFrom == fillValue);
            tmp = imwarp(tmp, tform, 'OutputView', outRef, 'FillValues', nan);
            outImage(tmp ~= 0) = nan;
        end

    end




end

