%OUT = COMPUTE1(OPERATION,IN,OUT_DIP_TYPE)
%    This is the basis function for all the overloaded unary operators.
%    OUT_DIP_TYPE is optional, and can be used to overrule the output
%    type generated by the operation. OUT_DIP_TYPE is either a specific
%    DIPlib data type, or one of: 'forcedouble' or 'keepsame'.

% (C) Copyright 1999-2011, Cris Luengo.
% Centre for Image Analysis, Uppsala, Sweden.
%
% Cris Luengo, June 2011, based on the original DO1INPUT, COMPUTE1, and
% COMPUTE1ARRAY.

function img = compute1(operation,img,force_type)
if nargin < 2
   error('Erroneus input.')
elseif nargin < 3
   force_type = '';
end

for ii=1:prod(imarsize(img))

   switch force_type
      case ''
         out_dip_type = di_findtype(img(ii).dip_type,'sfloat');
      case 'forcedouble'
         out_dip_type = di_findtype(img(ii).dip_type,'dfloat');
      case 'keepsame'
         out_dip_type = img(ii).dip_type;
      otherwise
         out_dip_type = force_type;
   end
   [out_type,complexity] = di_mattype(out_dip_type);
   indata = img(ii).data;

   direct = 0;
   if isa(indata,'double')
      direct = 1;
   else
      dosingle = 0;
      if matlabver_ge([7,0])
         if ~strcmp(out_type,'double')
            dosingle = 1;
         end
      end
      if dosingle & isa(indata,'single')
         direct = 1;
      else
         s = size(indata);
         N = prod(s);
         buffersize = dipgetpref('ComputationLimit');
         if dosingle
            buffersize = floor(buffersize/4);
            if N <= buffersize
               direct = 1;
               indata = single(indata);
            end
         else
            buffersize = floor(buffersize/8);
            if N <= buffersize
               direct = 1;
               indata = double(indata);
            end
         end
      end
   end

   if direct
      outdata = feval(operation,indata);
      %#function int8, uint8, int16, uint16, int32, uint32, int64, uint64, single, double
      outdata = feval(out_type,outdata);
   else
      outdata = di_create(s,out_type);
      if complexity, outdata = complex(outdata); end
      s = 1:buffersize:N;
      %disp(['We''re doing this in ',num2str(length(s)),' parts!']);
      if dosingle
         for ii=1:length(s)-1
            outdata(s(ii):s(ii+1)-1) = feval(operation,single(indata(s(ii):s(ii+1)-1)));
         end
         outdata(s(end):N) = feval(operation,single(indata(s(end):N)));
      else
         for ii=1:length(s)-1
            outdata(s(ii):s(ii+1)-1) = feval(operation,double(indata(s(ii):s(ii+1)-1)));
         end
         outdata(s(end):N) = feval(operation,double(indata(s(end):N)));
      end
   end

   % Make sure things are going fine if we produced a complex output:
   if ~isreal(outdata)
      if strcmp(out_dip_type,'bin') | strcmp(out_dip_type(2:4),'int')
         error('Operation unexpectedly produced complex values.')
      else
         out_dip_type = [out_dip_type(1),'complex'];
      end
   end

   img(ii).data = outdata;
   img(ii).dip_type = out_dip_type;
end
