% fdscale (Fermi-Dirac Rescaling) by Dan Piehl (2006)
%
% This is an image rescaler than can accept 2 or 3 dimensional
% arrays.  The third dimension will always be interpreted as
% a color subscript.  The function requires a rescaling factor
% (eg. pxlscale=4) and a "temperature" correction (eg. tempscale=1).
% A mean value correction factor can also be supplied.

function b = fdscale(a,pxlscale,tempscale,meancorrect)

  % Examine parameters and return error messages, if needed
  numDimsA = ndims(a);          % Get the number of input dimensions
  S1=size(a,1);                 % Get the number of rows
  S2=size(a,2);                 % Get the number of columns
  S3=size(a,3);                 % Get the number of color planes
  if (nargin < 3 )
      error('MATLAB:wavelet:NoInputs',['No input arguments specified.'])
  end
  if (numDimsA<2 || numDimsA>3) % Only 2 or 3 dimensional arrays are allowed
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.']) 
  end
  if not (tempscale>0 && pxlscale>0 && pxlscale-floor(pxlscale)==0)
      error('MATLAB:wavelet:NoInputs',['Invalid scaling parameters.'])
  end
  if not ( (S1==1 && S2>1) || (S1>1 && S2==1) || (S1>1 && S2>1) )
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.'])
  end

  % Determine min/max values (this is "linearized" for color images)
  R=2;                  % Reset largest search radius
  CB1=a;                % Start with local color basis
  CB2=a;
  for iteration=1:3
    CB2=1.1*lowpass(CB2);              % Identify local change
  end  
  c=sum((CB2-CB1).^2,3);    % Compute basis square distance
  c=max(cat(3,c,0.0001*ones(S1,S2)),[],3);  % Keep distance bounded
  MinVal=zeros(S1,S2);  % Reset minimum and maximum arrays
  MaxVal=zeros(S1,S2);
  R=2;                  % Reset largest search radius
  for C1=-R:R           % Loop through potential circular shift values
    for C2=-R:R         
      DS=C1^2+C2^2;     % Compute pixel distance squared
      if (DS>0 && DS<R^2)                  % For pixels within accepted radius
        t=circshift(a,[C1 C2 0]);          % Compute shifted matrix
        t=0.5*(c+sum((t-CB1).^2,3)-sum((t-CB2).^2,3))./c;
        MinVal=min(cat(3,MinVal,t),[],3);  % Compute minimum
        MaxVal=max(cat(3,MaxVal,t),[],3);  % Compute maximum
      end
    end
  end
  EpsVal=0.0001;
  MinVal=MinVal-EpsVal; % Add a small epsilon to avoid log(0) problems
  MaxVal=MaxVal+EpsVal;

  % Determine coefficients  (Array K1-K6 corresponds to values A-F)
  % NOTE:  P(X,Y) = AX2 + BXY + CY2 + DX + EY + F = ln(G/(array-H)-1)
  R=1;                               % Work only with a radius of 1 pixel
  K=zeros(S1,S2,6);                  % Reset coefficient array
  for C1=-R:R                        % Loop through shift values
    for C2=-R:R         
      DS=C1^2+C2^2;                  % Compute pixel distance squared
      if (DS<=2)                     % If within accepted radius squared
        t=circshift(a,[C1 C2 0]);    % Compute shifted matrix
        t=0.5*(c+sum((t-CB1).^2,3)-sum((t-CB2).^2,3))./c;
        t=real(log(((MaxVal-MinVal)./(t-MinVal))-1)); % Compute P(X,Y)
        if (C1==0 && C2==0)
          K(:,:,6)=t;               % Store F
        elseif ( C1==1  && C2==0 )
          K(:,:,1)=t;               % Store A+D+F              
        elseif ( C1==-1 && C2==0 )    
          K(:,:,4)=t;               % Store A-D+F
        elseif ( C1==0  && C2==1 )
          K(:,:,3)=t;               % Store C+E+F
        elseif ( C1==0  && C2==-1)
          K(:,:,5)=t;               % Store C-E+F
        elseif ( C1==C2)
          K(:,:,2)=K(:,:,2)+t;      % Modify B
        else
          K(:,:,2)=K(:,:,2)-t;      % Modify B
        end
      end
    end
  end
  K(:,:,2)=K(:,:,2)/4;              % Compute B
  t=K(:,:,1);                       % Get A+D+F
  K(:,:,1)=(t+K(:,:,4))/2-K(:,:,6); % Compute A
  K(:,:,4)=(t-K(:,:,4))/2;          % Compute D
  t=K(:,:,3);                       % Get C+E+F
  K(:,:,3)=(t+K(:,:,5))/2-K(:,:,6); % Compute C
  K(:,:,5)=(t-K(:,:,5))/2;          % Compute E

  % Build output array 
  R=pxlscale;                       % Rescale according to supplied parameter
  b=zeros(S1*R,S2*R,S3);            % Reset return array
  m=zeros(S1,S2,S3);                % Reset mean array
  for C1=1:R                        % Loop through shift values
    for C2=1:R
      F1=0.5-(C1-0.5)/R;            % Compute fractional pixel position
      F2=0.5-(C2-0.5)/R;            % ... then compute P(X+FX,Y+FY) array
      t=K(:,:,1)*(F1^2)+K(:,:,2)*F1*F2+K(:,:,3)*(F2^2)+K(:,:,4)*F1+K(:,:,5)*F2+K(:,:,6);
      t=MinVal+(MaxVal-MinVal)./(1+exp(t/tempscale));
      for color=1:S3
        c=(1-t).*CB1(:,:,color)+t.*CB2(:,:,color);      % Compute color value
        b(C1:R:(S1*R),C2:R:(S2*R),color)=c;             % Store color value
        m(:,:,color)=m(:,:,color)+c;                    % Increment mean value
      end
    end
  end
 
  % Normalize the mean values and check bounds
  if (nargin==4)                    % If a mean correction was supplied
    m=meancorrect*((m/(R^2))-a);    % Normalize the mean and compute adjustment
  else                              % Otherwise, if no mean correction...
    m=zeros(S1,S2,S3);              % Clear mean correction value
  end
  for C1=1:R                        % Loop through shift values
    for C2=1:R
      t=b(C1:R:(S1*R),C2:R:(S2*R),:)-m;     % Normalize the array
      t=max(cat(4,t,zeros(S1,S2,S3)),[],4); % Enforce minimum of zero
      t=min(cat(4,t,ones(S1,S2,S3)),[],4);  % Enforce maximum of one
      b(C1:R:(S1*R),C2:R:(S2*R),:)=t;       % Store corrected array
    end
  end