% Use ICM to remove the noise from the given image. % * covar is the known covariance of the Gaussian noise. % * max_diff is the maximum contribution to the potential % of the difference between two neighbouring pixel values. % * weight_diff is the weighting attached to the component of the potential % due to the difference between two neighbouring pixel values. % * iterations is the number of iterations to perform. function dst = restore_image(src, covar, max_diff, weight_diff, iterations) % Maintain two buffer images. % In alternate iterations, one will be the % source image, the other the destination. buffer = zeros(size(src,1), size(src,2), 2); buffer(:,:,1) = src; s = 2; d = 1; % This value is guaranteed to be larger than the % potential of any configuration of pixel values. V_max = (size(src,1) * size(src,2)) * ... ((256)^2 / (2*covar) + 4 * weight_diff * max_diff); for i = 1 : iterations % Switch source and destination buffers. if s == 1 s = 2; d = 1; else s = 1; d = 2; end % Vary each pixel individually to find the % values that minimise the local potentials. for r = 1 : size(src,1) for c = 1 : size(src,2) V_local = V_max; min_val = -1; for val = 0 : 255 % The component of the potential due to the known data. V_data = (val - src(r,c))^2 / (2 * covar); % The component of the potential due to the % difference between neighbouring pixel values. V_diff = 0; if r > 1 V_diff = V_diff + ... min( (val - buffer(r-1,c,s))^2, max_diff ); end if r < size(src,1) V_diff = V_diff + ... min( (val - buffer(r+1,c,s))^2, max_diff ); end if c > 1 V_diff = V_diff + ... min( (val - buffer(r,c-1,s))^2, max_diff ); end if c < size(src,2) V_diff = V_diff + ... min( (val - buffer(r,c+1,s))^2, max_diff ); end V_current = V_data + weight_diff * V_diff; if V_current < V_local min_val = val; V_local = V_current; end end buffer(r,c,d) = min_val; end end i end dst = buffer(:,:,d);