matlab-(含教程)基于OBNLM分块贝叶斯非局部优化图像去噪算法matlab仿真

function processedImg = BayesianNLM(img, blockSize, windowSize, gapBwnBlock, h) % Function: Bayesian-based nonlocal means algorithm % Input: % img: the origin image [data type: double or integer] % blockSize: size of the block % windowSize: size of the search window % gapBwnBlock: gap between the search blocks (in order to solve computational burden) % h: filtering parameter controlling the decay of the exponential function % % % _ _ _ _ _ _ _ imgSize_ _ _ _ _ _ _ __ % |* * * * * * * * * * | % |* * | % |* + + + * | % |* + + * | % |* + + + * | % |* blockSize * | % |* * | % |* * | % |* * windowSize * * | % | | % | | % |_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ | % % % Output: the despeckled image % Author: Shuwei Xing % Date: 2019-09-03 % Reference: Coup谷, Pierrick, et al. "Nonlocal means-based speckle filtering for ultrasound images." IEEE transactions on image processing 18.10 (2009): 2221-2229. %% % Define variables for algorithm originImg = double(img); % the raw image with noises imgSize = size(originImg,1); % Assumption: the shape of image is suqare processingImg = zeros(size(originImg)); % the value of processing indexImg = zeros(size(originImg)); % show the number of each pixel denoisedImg = zeros(size(originImg)); % the resotred(denoised) image epsilon = 10^(-13); % handle 0/0 case sideBlock = fix(blockSize/2); a = ceil(blockSize/2); b = imgSize - a + 1; A = ceil(windowSize/2); B = imgSize - A + 1; numCompletedPixel = 0; %% % Optimizaing the search space rowTotal = ceil(a / gapBwnBlock)* gapBwnBlock : gapBwnBlock : b ; colTotal = ceil(a / gapBwnBlock) * gapBwnBlock : gapBwnBlock : b; indexImgTotal = zeros(imgSize); indexImgTotal(rowTotal, colTotal) = 1; mask = logical(img); searchIndexMatrix =double(mask).* indexImgTotal; [rows, cols] = find(searchIndexMatrix); % Bayesian-based NonLocal Means Algorithm % tic % for row = ceil(a / gapBwnBlock)* gapBwnBlock : gapBwnBlock : b % for col = ceil(a / gapBwnBlock) * gapBwnBlock : gapBwnBlock : b for i = 1:size(cols, 1) col = cols(i); row = rows(i); % row = 12; % col = 200; % Get the block matrix of origin image % To do list: if the alogrithm can run fastly, write the specific % function to get the block. originBlock = originImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]); % % Discuss this solution within different conditions % % CONDITION 1 if row >= A && row <= B if col >= A && col <= B weightList1 = []; NLBlockWithoutNormalization1 = zeros(blockSize, blockSize); for row_NB = ceil((row - A + a )/ gapBwnBlock)* gapBwnBlock : gapBwnBlock : (row + A -a) for col_NB = ceil((col - A + a)/ gapBwnBlock) * gapBwnBlock : gapBwnBlock : (col + A - a) neighborBlock = originImg([(row_NB - sideBlock) : (row_NB + sideBlock) ], [(col_NB - sideBlock) : (col_NB + sideBlock)]); % Compute Pearson Distance distance = getPearsonDistance(originBlock, neighborBlock); % Compute the weight without normalization , just for single % neighborBlock weightWithoutNormalization = exp(-distance/(h^2)); weightList1 = [weightList1, weightWithoutNormalization]; % Compute the restored neighbor block NL_u_B without normalization subNLBlock1 = neighborBlock*weightWithoutNormalization; NLBlockWithoutNormalization1 = NLBlockWithoutNormalization1 + subNLBlock1; end end % Sum these weightes and get the final restored neighbor % block NL_u_B NLBlock = NLBlockWithoutNormalization1/sum(weightList1); % pixels'value of the NL_u_B is put in the processingImg processingImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) = processingImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) + NLBlock; % pixels' logic value of the NL_u_B is put in the indexImg indexImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) = indexImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) + ones(blockSize); % disp('run to this position') end % CONDITION 2 if col < A weightList2 = []; NLBlockWithoutNormalization2 = zeros(blockSize, blockSize); for row_NB = ceil((row - A + a )/ gapBwnBlock)* gapBwnBlock : gapBwnBlock : (row + A -a) for col_NB = ceil(a / gapBwnBlock)* gapBwnBlock : gapBwnBlock : (windowSize + 1 -a) neighborBlock = originImg([(row_NB - sideBlock) : (row_NB + sideBlock) ], [(col_NB - sideBlock) : (col_NB + sideBlock)]); % Compute Pearson Distance distance = getPearsonDistance(originBlock, neighborBlock); % Compute the weight without normalization , just for single % neighborBlock weightWithoutNormalization = exp(-distance/(h^2)); weightList2 = [weightList2, weightWithoutNormalization]; % Compute the restored neighbor block NL_u_B without normalization subNLBlock2 = neighborBlock*weightWithoutNormalization; NLBlockWithoutNormalization2 = NLBlockWithoutNormalization2 + subNLBlock2; end end % Sum these weights and get the final restored neighbor % block NL_u_B NLBlock = NLBlockWithoutNormalization2/sum(weightList2); % pixels'value of the NL_u_B is put in the processingImg processingImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) = processingImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) + NLBlock; % pixels' logic value of the NL_u_B is put in the indexImg indexImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) = indexImg([(row - sideBlock) : (row + sideBlock) ], [(col - sideBlock) : (col + sideBlock)]) + ones(blockSize);