本文实现Matlab函数psf2otf()重写为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)
padsize=size-psfsize
psf=np.lib.pad(psf,((0,padsize[0]),(0,padsize[1])),'constant')
#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
下面为使用样例（图像平滑）：
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)
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)
padsize=size-psfsize
psf=np.lib.pad(psf,((0,padsize[0]),(0,padsize[1])),'constant')
#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
该代码的处理效果：
