# L0 Smoothing 笔记（二）

    

本篇主要记录Image Smoothing via L0 Gradient Minimization这篇论文的MATLAB代码以及笔者改写的python代码。笔者初学python，若代码有不合理之处，还望读者指正

（一）原文的MATLAB代码：

function S = L0Smoothing(Im, lambda, kappa)

if ~exist('kappa','var')
kappa = 2.0;
end
if ~exist('lambda','var')
lambda = 2e-2;
end
S = im2double(Im);
betamax = 1e5;
fx = [1, -1];
fy = [1; -1];
[N,M,D] = size(Im);
sizeI2D = [N,M];
otfFx = psf2otf(fx,sizeI2D);
otfFy = psf2otf(fy,sizeI2D);
Normin1 = fft2(S);
Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
if D>1
Denormin2 = repmat(Denormin2,[1,1,D]);
end
beta = 2*lambda;
while beta < betamax
Denormin   = 1 + beta*Denormin2;
% h-v subproblem
h = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
if D==1
t = (h.^2+v.^2)<lambda/beta;
else
t = sum((h.^2+v.^2),3)<lambda/beta;
t = repmat(t,[1,1,D]);
end
h(t)=0; v(t)=0;
% S subproblem
Normin2 = [h(:,end,:) - h(:, 1,:), -diff(h,1,2)];
Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
FS = (Normin1 + beta*fft2(Normin2))./Denormin;
S = real(ifft2(FS));
beta = beta*kappa;
fprintf('.');
end
fprintf('\n');
end

(二)python代码：

def psf2otf(psf,size):
import numpy as np
if not(0 in psf):
#Pad the PSF to outsize
psf=np.double(psf)
psfsize=np.shape(psf)
psfsize=np.array(psfsize)
#Circularly shift otf so that the "center" of the PSF is at the (1,1) element of the array.
psf=np.roll(psf,-np.array(np.floor(psfsize/2),'i'),axis=(0,1))
#Compute the OTF
otf=np.fft.fftn(psf,axes=(0,1))
#Estimate the rough number of operations involved in the computation of the FFT.
nElem=np.prod(psfsize,axis=0)
nOps=0
for k in range(0,np.ndim(psf)):
nffts=nElem/psfsize[k]
nOps=nOps+psfsize[k]*np.log2(psfsize[k])*nffts
mx1=(abs(np.imag(otf[:])).max(0)).max(0)
mx2=(abs(otf[:]).max(0)).max(0)
eps= 2.2204e-16
if mx1/mx2<=nOps*eps:
otf=np.real(otf)
else:
otf=np.zeros(size)
return otf

def L0Smoothing(Im,lamda=2e-2,kappa=2.0):
import numpy as np
S=Im/255
betamax=1e5
fx=np.array([[1,-1]])
fy=np.array([[1],[-1]])
N,M,D=np.shape(Im)
sizeI2D=np.array([N,M])
otfFx=psf2otf(fx,sizeI2D)
otfFy=psf2otf(fy,sizeI2D)
Normin1=np.fft.fft2(S,axes=(0,1))
Denormin2=abs(otfFx)**2+abs(otfFy)**2
if D>1:
D2=np.zeros((N,M,D),dtype=np.double)
for i in range(D):
D2[:,:,i]=Denormin2
Denormin2=D2
beta=lamda*2
while beta<betamax:
Denormin=1+beta*Denormin2
#h-v subproblem
h1=np.diff(S,1,1)
h2=np.reshape(S[:,0],(N,1,3))-np.reshape(S[:,-1],(N,1,3))
h=np.hstack((h1,h2))
v1=np.diff(S,1,0)
v2=np.reshape(S[0,:],(1,M,3))-np.reshape(S[-1,:],(1,M,3))
v=np.vstack((v1,v2))
if D==1:
t=(h**2+v**2)<lamda/beta
else:
t=np.sum((h**2+v**2),2)<lamda/beta
t1=np.zeros((N,M,D),dtype=np.bool)
for i in range(D):
t1[:,:,i]=t
t=t1
h[t]=0
v[t]=0
#S subproblem
Normin2=np.hstack((np.reshape(h[:,-1],(N,1,3))-np.reshape(h[:,0],(N,1,3)),-np.diff(h,1,1)))
Normin2=Normin2+np.vstack((np.reshape(v[-1,:],(1,M,3))-np.reshape(v[0,:],(1,M,3)),-np.diff(v,1,0)))
FS=(Normin1+beta*np.fft.fft2(Normin2,axes=(0,1)))/Denormin
S=np.real(np.fft.ifft2(FS,axes=(0,1)))
beta*=kappa
print('.')
print('\n')
return S

def main():
import PIL
import numpy as np
import pylab
Im=np.array(PIL.Image.open("文件名"),'d')
S=L0Smoothing(Im,0.01)
pylab.imshow(S)

main()		

1、http://blog.csdn.net/panda1234lee/article/details/52825113