function [out hout]=resize_byN(in,NN,hin)
%
% function [out hout]=resize_byN(in,N,<,hin>)
%
% increases or decreases array size by integer N (N=1 or 0 is no change)
% Optionally converts the matlab .sir header array. Note that EASE2
% projections require that N be a power of 2
%
% increase: (N>0) (2==double size)
%  pixel replicates the input array 'in' to make larger array
%  or (N<0) smaller array where N is an integer 
%
% decrease: (N<0) (-2==half size)
%  averages NxN pixels of in input array 'in' into smaller output array
%

% written 17 Aug 2009 by DGL at BYU
% revised 23 Mar 2020 by DGL at BYU + handle more header variables

N=round(NN);
N1=abs(N);
if N1<2
  return;
end

% create output image and copy image data
if N>0  % increase size of array

  out=zeros(size(in)*N);
  [i j]=ind2sub(size(in),1:prod(size(in)));
  for x=1:N
    for y=1:N
      ind2=sub2ind(size(out),(i-1)*N+x,(j-1)*N+y);
      out(ind2)=in;
    end
  end

else   % decrease size of array

  [m n]=size(in);
  m2=floor(m/N1);
  n2=floor(n/N1);
  out=zeros([m2 n2]);
  [Y X]=meshgrid(1:m2,1:n2);
  ind=sub2ind([m2 n2],Y,X);

  % note: this computation does not handle nodata values
  sub=[1-N1:0];
  for x=1:N1
    for y=1:N1
      ind2=sub2ind([m n],N1*Y+sub(y),N1*X+sub(x));
      out(ind)=out(ind)+in(ind2);
    end
  end
  out=out/N1^2;

end

% if header is input, modify and return
if exist('hin')
  % pull header values for convenience
  nsx    = hin(1);
  nsy    = hin(2);
  xdeg   = hin(3);
  ydeg   = hin(4);
  ascale = hin(6);
  bscale = hin(7);
  a0     = hin(8);
  b0     = hin(9);
  iopt   = hin(17);

  if N<1, N=1/N1; end;
  
  nsx = floor(N * nsx);
  nsy = floor(N * nsy);
  
  % set header values, first line creates output header
  hout=setsirhead('nsx',hin,nsx);
  hout=setsirhead('nsy',hout,nsy);

  switch (iopt)
    case -1  % image only
      xdeg = N * xdeg;
      ydeg = N * ydeg;
      hout=setsirhead('xdeg',hout,xdeg);  
      hout=setsirhead('ydeg',hout,ydeg);  
    case 0   % lat/lon
      xdeg = N * xdeg;
      ydeg = N * ydeg;
      hout=setsirhead('xdeg',hout,xdeg);  
      hout=setsirhead('ydeg',hout,ydeg);  
    case 1   % Lambert
      ascale = N * ascale; 
      bscale = N * bscale; 
      hout=setsirhead('ascale',hout,ascale);  
      hout=setsirhead('bscale',hout,bscale);  
    case 2   % Lambert
      ascale = N * ascale; 
      bscale = N * bscale; 
      hout=setsirhead('ascale',hout,ascale);  
      hout=setsirhead('bscale',hout,bscale);  
    case 5   % polar stereographic
      ascale = ascale / N; 
      bscale = bscale / N; 
      hout=setsirhead('ascale',hout,ascale);  
      hout=setsirhead('bscale',hout,bscale);  
    case 8   % EASE2 grid
      n=floor(log(N)/log(2));
      ascale = ascale + n;
      hout=setsirhead('ascale',hout,ascale);
    case 9   % EASE2 grid
      n=floor(log(N)/log(2));
      ascale = ascale + n;
      hout=setsirhead('ascale',hout,ascale);
    case 10  % EASE2 grid
      n=floor(log(N)/log(2));
      ascale = ascale + n;
      hout=setsirhead('ascale',hout,ascale);
    case 11  % EASE1 grid
      n=floor(log(N)/log(2));
      ascale = ascale / N;
      bscale = N * bscale;
    case 12  % EASE1 grid
      n=floor(log(N)/log(2));
      ascale = ascale / N;
      bscale = N * bscale;
    case 13  % EASE1 grid
      n=floor(log(N)/log(2));
      ascale = ascale / N;
      bscale = N * bscale;
    otherwise% unknown
  end
end

