function [finalout]=dbl(FFTrv,FFTiv,IFFTiv,size1,size2,isreal1,isreal2,finetuning)
% [finalout]=dbl(FFTrv,FFTiv,IFFTiv,size1,size2,isreal1,isreal2)
% Determine the best parameters for Overlap-Add FFT-based convolution.
%
% INPUT
% FFTrv:      vector with costs of FFT for real 1d vectors
% FFTiv:      vector with costs of FFT for complex 1d vectors
% IFFTiv:     vector with costs of IFFT for complex 1d vectors
% size1:      size(first_image)
% size2:      size(second_image)
% isreal1:    1 if first image is real, 0 otherwise (complex)
% isreal2:    1 if second image is real, 0 otherwise (complex)
% finetuning: once theoretical parameters have been found, in order to
%             optimize the final output, a finer search is performed.
%             This step will is computationally expensive.
%
% OUTPUT
% out:    the optimized parameters:
%         out.inverse:     if 1 the two input have to be inverted
%         out.fftxfirst:   if one the image has to be fft first along
%                          x-dimension
%         out.ifftxfirst:  if one the product of spectra has to be ifft
%                          first along x-dimension. NOTE: Always set to 1.
%         out.nfftx:       the best length for fft transform along
%                          x-dimension
%         out.nffty:       the best length for fft transform along
%                          y-dimension
%         out.filterxfirst if 1 the filter has to be fft fisrt alng
%                          x-dimension
%
%
if isreal1 && isreal2
    % A real image, B real filter
    out{1} = test(FFTrv,FFTiv,FFTrv,FFTiv,IFFTiv,size1,size2);
    % B real image, A real filter
    out{2} = test(FFTrv,FFTiv,FFTrv,FFTiv,IFFTiv,size2,size1);
end
if isreal1 && ~isreal2
    % A real image, B complex filter
    out{1} = test(FFTrv,FFTiv,FFTiv,FFTiv,IFFTiv,size1,size2);
    % B complex image, A real filter
    out{2} = test(FFTiv,FFTiv,FFTrv,FFTiv,IFFTiv,size2,size1);
end
if ~isreal1 && isreal2
    % A complex image, B real filter
    out{1} = test(FFTiv,FFTiv,FFTrv,FFTiv,IFFTiv,size1,size2);
    % B real image, A complex filter
    out{2} = test(FFTrv,FFTiv,FFTiv,FFTiv,IFFTiv,size2,size1);
end
if ~isreal1 && ~isreal2
    % A complex image, B complex filter
    out{1} = test(FFTiv,FFTiv,FFTiv,FFTiv,IFFTiv,size1,size2);
    % B complex image, A complex filter
    out{2} = test(FFTiv,FFTiv,FFTiv,FFTiv,IFFTiv,size2,size1);
end
if nargin<8 || ((nargin == 8) && (finetuning == 0))
    [finalout] = selectout(out);
    return;
else
    disp('Optimization of parameters in progress... Please wait.');
    if out{1}.flops<out{2}.flops
        pos     = 1;
        inverse = 0;
        a       = rand(size1);
        b       = rand(size2);
        if ~isreal1
            a = a + i*rand(size1);
        end
        if ~isreal2
            b = b + i*rand(size2);
        end
    else
        pos     = 2;
        inverse = 1;
        a       = rand(size2);
        b       = rand(size1);
        if ~isreal2
            a = a + i*rand(size2);
        end
        if ~isreal1
            b = b + i*rand(size1);
        end
    end
    coeff = 1.2;
    m     = out{pos}.flopmatrix <= out{pos}.flops*coeff;
    pos   = find(m);
    [x,y] = ind2sub(size(m),pos);
    L = length(pos);
    time_req = zeros(L,1);
    h = waitbar(0,'Please wait...');
    for ii=1:L
        tic;
        fftolam(a,b,x(ii),y(ii));
        tr           = toc;
        time_req(ii) = tr;
        waitbar(ii/L)
    end
    close(h);
    [valmin,posval] = min(time_req);
    xok             = x(posval);
    yok             = y(posval);
    
    [fftxfirst,filterxfirst] = findpm(FFTrv,FFTiv,IFFTiv,size(a),size(b),isreal(a),isreal(b),xok,yok);
end
finalout.fftxfirst    = fftxfirst;
finalout.ifftxfirst   = 1;%out{pos}.ifftxfirst;
finalout.nfftx        = xok;
finalout.nffty        = yok;
finalout.filterxfirst = filterxfirst;
finalout.inverse      = inverse;
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [finalout] = selectout(out)
% If no further fine tuning is selected, this function is used to find the
% best parameters that have been just calculated.
if out{1}.flops < out{2}.flops
    pos = 1;
else
    pos = 2;
end
if pos == 1
    finalout.inverse = 0;
else
    finalout.inverse = 1;
end
finalout.fftxfirst    = out{pos}.fftxfirst;
finalout.ifftxfirst   = 1;%out{pos}.ifftxfirst;
finalout.nfftx        = out{pos}.nfftx;
finalout.nffty        = out{pos}.nffty;
finalout.filterxfirst = out{pos}.filterxfirst;
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [out] = test(fft_A_1,fft_A_2,fft_B_1,fft_B_2,ifft_1,sizeA,sizeB)
% This function detremines the number of flops for overlap-add method
% varying the length of vectors that have to be FFT- and IFFT- transformed.
% For example fft2 of a matrix NxM is equivalent to do N FFTs with vectors
% of length M and then M FFTs of vectors of length N. See fft2 in Matlab
% help.
ax = sizeA(1);
ay = sizeA(2);
bx = sizeB(1);
by = sizeB(2);

infinitevalue = 99*10^99;
val0          = infinitevalue;

costcomplexsum  = 2;
costcomplexprod = 6;

L = size(fft_A_1,1);

out.flopmatrix = 1.5270e+050*ones(L,L);


for ii=1:L
    for jj=1:L
        Lx    = ii-bx+1;
        Ly    = jj-by+1;
        if Lx>0 && Ly>0
            nx    = ceil(ax/Lx);
            ny    = ceil(ay/Ly);

            cv1 = nx*ny*(Lx*fft_A_1(jj,2) + jj*fft_A_2(ii,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod)) + (bx*fft_B_1(jj,2) + jj*fft_B_2(ii,2));
            cv3 = nx*ny*(Ly*fft_A_1(ii,2) + ii*fft_A_2(jj,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod)) + (bx*fft_B_1(jj,2) + jj*fft_B_2(ii,2));
            cv5 = nx*ny*(Lx*fft_A_1(jj,2) + jj*fft_A_2(ii,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod)) + (by*fft_B_1(ii,2) + ii*fft_B_2(jj,2));
            cv7 = nx*ny*(Ly*fft_A_1(ii,2) + ii*fft_A_2(jj,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod)) + (by*fft_B_1(ii,2) + ii*fft_B_2(jj,2));

            out.flopmatrix(ii,jj) = (cv1+cv3+cv5+cv7)/4;

            if cv1<val0
                val0         = cv1;
                fftxfirst    = 0;
                filterxfirst = 0;
                nfftx        = ii;
                nffty        = jj;
            end
            if cv3<val0
                val0         = cv3;
                fftxfirst    = 1;
                filterxfirst = 0;
                nfftx        = ii;
                nffty        = jj;
            end
            if cv5<val0
                val0         = cv5;
                fftxfirst    = 0;
                filterxfirst = 1;
                nfftx        = ii;
                nffty        = jj;
            end
            if cv7<val0
                val0         = cv7;
                fftxfirst    = 1;
                filterxfirst = 1;
                nfftx        = ii;
                nffty        = jj;
            end
        end
    end
end

out.flops        = val0;
out.fftxfirst    = fftxfirst;
out.nfftx        = nfftx;
out.nffty        = nffty;
out.filterxfirst = filterxfirst;
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [out]=fftolam(a,b,siz1,siz2)
% [out]=fftolam(a,b,siz1,siz2)
% Simple Overlap-add method FFT-based 2D convolution
% INPUT
% a:     first image (2D double matrix)
% b:     second image (2D double matrix)
% siz1:  the specified x-length for FFT along x-dimension (siz1 > size(b,1))
% siz2:  the specified x-length for FFT along y-dimension (siz2 > size(b,2))
% OUTPUT
% out:   2D convolution of a and b matrices: out = conv2(a,b);

[ax,ay]       = size(a);
[bx,by]       = size(b);
dimx          = ax+bx-1;
dimy          = ay+by-1;
nfftx         = siz1;
nffty         = siz2;
Lx            = nfftx-bx+1;
Ly            = nffty-by+1;
B             = fft2(b,nfftx,nffty);
out           = zeros(dimx,dimy);
x0 = 1;
while x0 <= ax
    x1   = min(x0+Lx-1,ax);
    y0   = 1;
    endx = min(dimx,x0+nfftx-1);
    while y0 <= ay
        y1                   = min(y0+Ly-1,ay);
        endy                 = min(dimy,y0+nffty-1);
        X                    = fft2(a(x0:x1,y0:y1),nfftx,nffty);
        Y                    = ifft2(X.*B);
        out(x0:endx,y0:endy) = out(x0:endx,y0:endy)+Y(1:(endx-x0+1),1:(endy-y0+1));
        y0                   = y0+Ly;
    end
    x0 = x0+Lx;
end
if ~(any(any(imag(a)))||any(any(imag(b))))
    out=real(out);
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [fftxfirst,filterxfirst]=findpm(FFTrv,FFTiv,IFFTiv,size1,size2,isreal1,isreal2,xok,yok)
% If the fine tuning option is selectd, this function determines the best parameters
% for overlap-add method.
ii = xok;
jj = yok;

ax = size1(1);
ay = size1(2);
bx = size2(1);
by = size2(2);
Lx = ii-bx+1;
Ly = jj-by+1;
nx = ceil(ax/Lx);
ny = ceil(ay/Ly);

costcomplexsum  = 2;
costcomplexprod = 6;

ifft_1          = IFFTiv;
if isreal2
    fft_B_1 = FFTrv;
    fft_B_2 = FFTiv;
else
    fft_B_1 = FFTiv;
    fft_B_2 = FFTiv;
end
if isreal1
    fft_A_1 = FFTrv;
    fft_A_2 = FFTiv;
else
    fft_A_1 = FFTiv;
    fft_A_2 = FFTiv;
end
if (bx*fft_B_1(jj,2) + jj*fft_B_2(ii,2))<(by*fft_B_1(ii,2) + ii*fft_B_2(jj,2))
    filterxfirst = 0;
else
    filterxfirst = 1;
end 

p1 = nx*ny*(Lx*fft_A_1(jj,2) + jj*fft_A_2(ii,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod));
p2 = nx*ny*(Ly*fft_A_1(ii,2) + ii*fft_A_2(jj,2) + jj*ifft_1(ii,2) + ii*ifft_1(jj,2) + ii*jj*(costcomplexsum+costcomplexprod));
if p1<p2
    fftxfirst = 0;
else
    fftxfirst = 1;
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------