function [FinalLattice,FinalFood,FinalGarbage,Gsums,Rsums,FoodRunning,GarbageRunning] = DoubleBrothLoopF02(InitialLattice,MaskLattice,InitialFoodLattice,InitialGarbageLattice,fc,fd,g,rc,rd,rg,rho,LysingRate,DFchannel,DGchannel,DFres,FoodDiffusionSigma,GarbageDiffusionSigma,MHPCoords,N,FilePrefix)
%   ================================================
%   David Liao
%   2008 February 10
%   DoubleBrothLoopF
%   ================================================
%
%   Requires DoubleBrothF.m and IsingTrioFluorescence.m
%   [FinalLattice,FinalFood,FinalGarbage,Gsums,Rsums,FoodRunning,GarbageRun
%   ning] =
%   DoubleBrothLoopF01(InitialLattice,MaskLattice,InitialFoodLattice,Initia
%   lGarbageLattice,fc,fd,g,rc,rd,rg,rho,LysingRate,DFchannel,DGchannel,DFr
%   es,FoodDiffusionSigma,GarbageDiffusionSigma,MHPCoords,N,FilePrefix)

L = size(InitialLattice);

%   Nearest neighbor filter: do fft once
NeighborFilter = zeros(L(1)+2,L(2)+2,L(3));
NeighborFilter(2,1,1) = 1; NeighborFilter(1,2,1) = 1; NeighborFilter(1,1,2) = 1;
NeighborFilter(L(1)+2,1,1) = 1; NeighborFilter(1,L(2)+2,1) = 1; NeighborFilter(1,1,L(3)) = 1;
NeighborFilter(1,1,1) = 1;
NeighborFilterfft = fftn(NeighborFilter);

GreatestDiffLength = min(max(FoodDiffusionSigma,GarbageDiffusionSigma),MHPCoords(1,4)/2);

CAPMU = [MHPCoords(1,4)+2*GreatestDiffLength MHPCoords(1,5)+2*GreatestDiffLength L(3)];

%   Food diffusion filter: do fft once
FoodGaussian = zeros(CAPMU);
for i = 1:CAPMU(1)
    for j = 1:CAPMU(2)
        for k = 1:CAPMU(3)
            FoodGaussian(i,j,k) = exp(-(i^2 + j^2 + k^2)/FoodDiffusionSigma^2);
        end
    end
end

FoodGaussian = FoodGaussian + FoodGaussian(CAPMU(1):-1:1,:,:) + FoodGaussian(:,CAPMU(2):-1:1,:) + FoodGaussian(:,:,CAPMU(3):-1:1) + FoodGaussian(CAPMU(1):-1:1,CAPMU(2):-1:1,:) + FoodGaussian(:,CAPMU(2):-1:1,CAPMU(3):-1:1) + FoodGaussian(CAPMU(1):-1:1,:,CAPMU(3):-1:1) + FoodGaussian(CAPMU(1):-1:1,CAPMU(2):-1:1,CAPMU(3):-1:1);
FoodGaussian = FoodGaussian/sum(sum(sum(FoodGaussian)));
FoodGaussianfft = fftn(FoodGaussian);

%   Garbage diffusion filter: do fft once
GarbageGaussian = zeros(CAPMU);
for i = 1:CAPMU(1)
    for j = 1:CAPMU(2)
        for k = 1:CAPMU(3)
            GarbageGaussian(i,j,k) = exp(-(i^2 + j^2 + k^2)/GarbageDiffusionSigma^2);
        end
    end
end

GarbageGaussian = GarbageGaussian + GarbageGaussian(CAPMU(1):-1:1,:,:) + GarbageGaussian(:,CAPMU(2):-1:1,:) + GarbageGaussian(:,:,CAPMU(3):-1:1) + GarbageGaussian(CAPMU(1):-1:1,CAPMU(2):-1:1,:) + GarbageGaussian(:,CAPMU(2):-1:1,CAPMU(3):-1:1) + GarbageGaussian(CAPMU(1):-1:1,:,CAPMU(3):-1:1) + GarbageGaussian(CAPMU(1):-1:1,CAPMU(2):-1:1,CAPMU(3):-1:1);
GarbageGaussian = GarbageGaussian/sum(sum(sum(GarbageGaussian)));
GarbageGaussianfft = fftn(GarbageGaussian);

%   Write first picture
[ZerothColorImage,ZerothGsum,ZerothRsum] = Fluorescence(InitialLattice);
imwrite(ZerothColorImage,[FilePrefix '000000.png']);
CurrentLattice = InitialLattice;
CurrentFoodLattice = InitialFoodLattice;
CurrentGarbageLattice = InitialGarbageLattice;

Gsums = ZerothGsum;
Rsums = ZerothRsum;

FoodRunning = sum(sum(sum(CurrentFoodLattice)));
GarbageRunning = sum(sum(sum(CurrentGarbageLattice)));

tic;

for i = 1:N
    [CurrentLattice,CurrentFoodLattice,CurrentGarbageLattice] = DoubleBrothF02(CurrentLattice,MaskLattice,CurrentFoodLattice,CurrentGarbageLattice,fc,fd,g,rc,rd,rg,rho,LysingRate,DFchannel,DGchannel,DFres,FoodGaussianfft,FoodDiffusionSigma,GarbageGaussianfft,GarbageDiffusionSigma,NeighborFilterfft,MHPCoords);
    [CurrentColorImage,Gsum,Rsum] = Fluorescence(CurrentLattice);
    imwrite(CurrentColorImage,[FilePrefix num2str(i,'%06g') '.png']);
    Gsums = [Gsums; Gsum];
    Rsums = [Rsums; Rsum];
    FoodRunning = [FoodRunning; sum(sum(sum(CurrentFoodLattice)))];
    GarbageRunning = [GarbageRunning; sum(sum(sum(CurrentFoodLattice)))];
    ElapsedTime = toc;
    
    RemainingTime = (N - i)*ElapsedTime/i;
    fprintf(['Remaining time: ' ,num2str(RemainingTime), ' s.  Elapsed time: ' ,num2str(ElapsedTime), ' s.  \n']);
end

FinalLattice = CurrentLattice;
FinalFood = CurrentFoodLattice;
FinalGarbage = CurrentGarbageLattice;