%OUT = COMPUTE2(OPERATION,IN1,IN2)
%    This is the basis function for all the overloaded operators and
%    functions.
%OUT = COMPUTE2(OPERATION,IN1,IN2,DIMS,OUT_DIP_TYPE,OUT_PHYS)
%    This is the form if you already called DO2INPUTS(). Use this if
%    the format returned is not correct.

% (C) Copyright 1999-2008               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Cris Luengo, April 2000.
% 1 June 2000:     Using MATTYPE created this week.
% 2 June 2000:     Using FEVAL instead of BUILTIN.
% 18 January 2001: Added DIMS parameter. Changed all calls to this
%                  function accordingly.
% 14 August 2006:  Fixed output image being reshaped if input is 1xN.
% February 2008:   Adding pixel dimensions and units to dip_image. (BR)
% 11 March 2008:   Since MATLAB 7.0, computations can be done with singles.

function out = compute2(operation,in1,in2,dims,out_dip_type,out_phys)
if nargin == 3
   [in1,in2,dims,out_dip_type,out_phys] = do2inputs(in1,in2);
elseif nargin ~= 6
   error('Erroneus input.')
end
[out_type,complexity] = di_mattype(out_dip_type);
direct = 0;
if isa(in1,'double') & isa(in2,'double')
   direct = 1;
else
   dosingle = 0;
   if matlabversion>=7.0
      if ~strcmp(out_type,'double')
         dosingle = 1;
      end
   end
   if dosingle & (isa(in1,'single')|isa(in1,'double')) & (isa(in2,'single')|isa(in2,'double'))
      % single+double -> single
      % double+single -> single
      % single+single -> single
      direct = 1;
   else
      isscalar = 0;
      s1 = size(in1); if prod(s1) == 1, isscalar = 1; end % empty arrays are not scalar
      s2 = size(in2); if prod(s2) == 1, isscalar = 2; end
      if isscalar==1
         N = prod(s2);
         s = s2;
      else % isscalar==2 or 0
         N = prod(s1);
         s = s1;
      end
      buffersize = dipgetpref('ComputationLimit');
      if dosingle
         buffersize = floor(buffersize/4);
         if N <= buffersize
            direct = 1;
            in1 = single(in1);
            in2 = single(in2);
         end
      else
         buffersize = floor(buffersize/8);
         if N <= buffersize
            direct = 1;
            in1 = double(in1);
            in2 = double(in2);
         end
      end
   end
end
if direct
   out = feval(operation,in1,in2);
   %#function int8, uint8, int16, uint16, int32, uint32, int64, uint64, single, double
   out = feval(out_type,out);
else
   out = di_create(s,out_type);
   if complexity, out = complex(out); end
   s = 1:buffersize:N;
   %disp(['We''re doing this in ',num2str(length(s)),' parts!']);
   if isscalar==1
      if dosingle
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               single(in1),single(in2(s(ii):s(ii+1)-1)));
         end
         out(s(end):N) = feval(operation,single(in1),single(in2(s(end):N)));
      else
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               double(in1),double(in2(s(ii):s(ii+1)-1)));
         end
         out(s(end):N) = feval(operation,double(in1),double(in2(s(end):N)));
      end
   elseif isscalar==2
      if dosingle
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               single(in1(s(ii):s(ii+1)-1)),single(in2));
         end
         out(s(end):N) = feval(operation,single(in1(s(end):N)),single(in2));
      else
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               double(in1(s(ii):s(ii+1)-1)),double(in2));
         end
         out(s(end):N) = feval(operation,double(in1(s(end):N)),double(in2));
      end
   else
      if dosingle
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               single(in1(s(ii):s(ii+1)-1)),single(in2(s(ii):s(ii+1)-1)));
         end
         out(s(end):N) = feval(operation,single(in1(s(end):N)),single(in2(s(end):N)));
      else
         for ii=1:length(s)-1
            out(s(ii):s(ii+1)-1) = feval(operation,...
               double(in1(s(ii):s(ii+1)-1)),double(in2(s(ii):s(ii+1)-1)));
         end
         out(s(end):N) = feval(operation,double(in1(s(end):N)),double(in2(s(end):N)));
      end
   end
end
% The next step not only creates the dip_image object, but also sets the
% datatype correctly. That is, we expect no conversion is done anymore,
% but a binary image will be labelled as such.
out = dip_image('trust_me',out,out_dip_type,dims,out_phys);
