function PPPM = Func_Wiener_Filter(RRR,SSS,lambda)

[Rsize Csize Hsize] = size(RRR);
Tsize = Rsize;
Thalf = (Tsize-1)/2;
Tcent = Thalf + 1;
PZ = 121; 
Th = 3;

RRR(RRR<0.001*max(max(max(RRR)))) = 0;
figure,imagesc(squeeze(sum(RRR(:,:,PZ-Th:PZ+Th),3))); axis image; colormap(gray); drawnow

III = (abs(fftshift(fftn(RRR)))).^2;

Rin = 40;
BS = ones(Tsize,Tsize,Tsize);
for a = 1:Tsize
    for b = 1:Tsize
        for c = 1:Tsize
            if (sqrt((a-Tcent)^2+(b-Tcent)^2+(c-Tcent)^2)<=Rin)
                BS(a,b,c) = 0;
            end
        end
    end
end

DDD = III;
Inoise = zeros(1,209);
Cnoise = zeros(1,209);
for a = 1:209
    for b = 1:209
        for c = 1:209
            if DDD(a,b,c)~=0
                D = sqrt((a-Tcent)^2+(b-Tcent)^2+(c-Tcent)^2);
                px = round(D)+1;
                Inoise(px) = Inoise(px) + DDD(a,b,c);
                Cnoise(px) = Cnoise(px) + 1;
            end 
        end
    end
end
InoiseAVE = zeros(1,209);
for a = 1:209
    if Cnoise(a)~=0
        InoiseAVE(a) = Inoise(a) / Cnoise(a);
    end
end
InoiseAVE = smooth(InoiseAVE,3);

NNN = zeros(Tsize,Tsize,Tsize);
del = 1e-10;
for a = 1:Tsize
    for b = 1:Tsize
        for c = 1:Tsize
            D = sqrt((a+del-Tcent)^2+(b+del-Tcent)^2+(c+del-Tcent)^2);
            Dd = floor(D);
            Sd = D - Dd;
            Du = ceil(D);
            Su = Du - D;
            if D<=206
                NNN(a,b,c) = 1/(Sd+Su) * (Sd*InoiseAVE(Du+1) + Su*(InoiseAVE(Dd+1)));
            else
                NNN(a,b,c) = InoiseAVE(207);
            end
        end
    end
end

KKK = fftshift(fftn(RRR));
AAA = (abs(fftshift(fftn(RRR)))).^2;
PPP = AAA - lambda*NNN;
PPP(PPP<0) = 0;
FL = PPP ./ (PPP+lambda*NNN);

KKKM = FL .* KKK;
KKKM(Tcent,Tcent,Tcent) = 0;
PPPM = real(ifftn(ifftshift(KKKM)));
PPPM(SSS==0) = 0;

figure,imagesc(squeeze(sum(PPPM(:,:,PZ-Th:PZ+Th),3))); axis image; colormap(gray); drawnow

