function [bgmx,bgmn,static,mod1,mod2,wmod,mn,mx]=...
  bgcalc4(nn,bb,Nbg,Nx,Ny,Nc,wid,f1,f2,mncout,mxcout);

MINPIXVAL=0;
MAXPIXVAL=255;

[nx ny]=ndgrid(1:Nx,1:Ny);
fprintf(1,'processing.');
cnn=cumsum(nn,3);
[jj mn]=max(cnn>=Nbg*f1,[],3);
fprintf(1,'.');
[jj mx]=max(cnn>=Nbg*(1-f1),[],3);
fprintf(1,'.');
irn=fix(((mn)+(mx))/2);
idn=fix(((mx)-(mn)+1));
[jj mod1]=max(nn,[],3);
fprintf(1,'.');
wmod=(reshape(cnn(sub2ind(size(nn),nx,ny,clip(mod1+wid,1,Nc))),Nx,Ny))-(reshape(cnn(sub2ind(size(nn),nx,ny,clip(mod1-wid,1,Nc))),Nx,Ny));
fprintf(1,'.');
clear('cnn');
bmod1=reshape(repmat(mod1,1,Nc),Nx,Ny,Nc);
zzz=reshape(repmat(bb,Nx*Ny,1),Nx,Ny,Nc);
ii=(zzz<bmod1-wid) | (zzz>bmod1+wid);
clear('bmod1','zzz');
[jj mod2]=max(nn.*ii,[],3);
fprintf(1,'.');
clear('ii');
static=wmod>Nbg*(1-f2);
bgmx=max(mod1,mod2);
bgmn=min(mod1,mod2);
bgmean=(bgmx+bgmn)/2;
if(length(find(~static))==0)
  [junk p1]=max(hist(mod1(:),bb));
  ii=(bb<bb(p1)-wid) | (bb>bb(p1)+wid);
  [junk p2]=max(hist(mod1(:),bb).*ii);
  cut=(bb(p1)+bb(p2))/2;
else
  cut=mean(bgmean(~static));
end
mmmx=(((bgmx-mx)./(bgmx-bgmn).*~static)<mncout);
mmmn=(((bgmx-mn)./(bgmx-bgmn).*~static)>mxcout);
bgmn(static)=mod1(static);
bgmx(static)=mod1(static);
bgmx((static) & (mod1<=cut))=MAXPIXVAL;
bgmn((static) & (mod1>cut))=MINPIXVAL;
bgmx(mmmx)=mx(mmmx);
bgmn(mmmn)=mn(mmmn);
fprintf(1,'Done\n');
  