图像各向异性平滑滤波

• 介绍
• • 原理解析
  • 实现

介绍

图像处理中先验知识较少的情况下，需要在保留图像边缘信息下同时对图像进行平滑，可以采用将越明显的边缘减少平滑率的方式进行滤波，就是今天要讲的PM各向异性扩散滤波。
           

https://max.book118.com/html/2018/0913/6034150213001215.shtm

https://wenku.baidu.com/view/9c2897cbf78a6529657d5391.html

原理解析

通过迭代计算最佳扩散平滑尺度，

一般迭代1到100次，可以通过获取的最佳效果来决定。

实现

function [filtmap2d] = anisotropicfiltering(map2d, iterations, lambda, contrastK, smoothsigma, Dosharpenfirst, sharpensigma)
%ANISOTROPICFILTER 此处显示有关此函数的摘要
%   此处显示详细说明

	if iterations <= 0
	    iterations = 100;
	end
	if smoothsigma <= 0
	    smoothsigma = 2;
	end
	smoothfilt = fspecial('gaussian',ceil(6*smoothsigma),smoothsigma);
	contrastK_square = contrastK^2;
	if Dosharpenfirst == true
	    result_map2d = tools_sharpen(grapmap3d, sharpensigma);
	    result_map2d(result_map2d < 0) = 0;
	    result_map2d = result_map2d / 256;
	else
	    result_map2d = map2d / 256;
	end
	for iter=1:iterations
	    Delta_N = result_map2d(1:end-2,2:end-1) - result_map2d(2:end-1,2:end-1);
	    Delta_S = result_map2d(3:end,2:end-1) - result_map2d(2:end-1,2:end-1);
	    Delta_E = result_map2d(2:end-1,3:end) - result_map2d(2:end-1,2:end-1);
	    Delta_W = result_map2d(2:end-1,1:end-2) - result_map2d(2:end-1,2:end-1);
	
	% % %   PM filter
	    g_N = exp(-(Delta_N.^2)/contrastK_square);
	    g_S = exp(-(Delta_S.^2)/contrastK_square);
	    g_E = exp(-(Delta_E.^2)/contrastK_square);
	    g_W = exp(-(Delta_W.^2)/contrastK_square);
	
	% % % %   catte et al filter (gaussian smoothing)
	%     temp_map2d = imfilter(result_map2d, smoothfilt);
	%     temp_Delta_N = temp_map2d(1:end-2,2:end-1) - temp_map2d(2:end-1,2:end-1);
	%     temp_Delta_S = temp_map2d(3:end,2:end-1) - temp_map2d(2:end-1,2:end-1);
	%     temp_Delta_E = temp_map2d(2:end-1,3:end) - temp_map2d(2:end-1,2:end-1);
	%     temp_Delta_W = temp_map2d(2:end-1,1:end-2) - temp_map2d(2:end-1,2:end-1);
	%     g_N = 1 ./ (1+(temp_Delta_N.^2)/contrastK_square);
	%     g_S = 1 ./ (1+(temp_Delta_S.^2)/contrastK_square);
	%     g_E = 1 ./ (1+(temp_Delta_E.^2)/contrastK_square);
	%     g_W = 1 ./ (1+(temp_Delta_W.^2)/contrastK_square);
	    
	    result_map2d(2:end-1,2:end-1) = result_map2d(2:end-1,2:end-1) + lambda .* (g_N.*Delta_N + g_S.*Delta_S + g_E.*Delta_E + g_W.*Delta_W);
	    
	    if mod(iter,5) == 0
	        figure;imshow(result_map2d(2:end-1,2:end-1),[]);
	    end
	end
	
	filtmap2d = result_map2d;
	filtmap2d(1,:) = result_map2d(2,:);
	filtmap2d(end,:) = result_map2d(end-1,:);
	filtmap2d(:,1) = result_map2d(:,2);
	filtmap2d(:,end) = result_map2d(:,end-1);
	filtmap2d = (filtmap2d - min(filtmap2d(:))) / (max(filtmap2d(:)) - min(filtmap2d(:))) * 256;
end