用matlab写的径向分布函数RDF

RDF是径向分布函数Radical distribution function的缩写，指的是给定一个空间，在此空间以一个对象为中心，去寻找周围对象的的概率。对于分子模拟的径向分布函数实则也是求解粒子在周期性边界盒子的区域密度和全局密度的比值。
其中：
区域密度 = 每一个球壳的数密度 / 球壳体积。
全局密度 = 盒子里面的粒子总数 / 盒子的体积
径向分布 = 区域密度 / 全局密度

但是， 本文用于计算的RDF是基于二维的。目前很多仿真软件的到结果是：一张图片上有很多个原子。本文是用于计算图像上原子的RDF.
则：
区域密度 = 每个圆内包含的原子个数 / 圆的面积。
全局密度 = 图片上原子的总数 / 图片面积
径向分布 = 区域密度 / 全局密度

matlab代码演示：
Circle_detect.m（检测输入图片上的圆）

function [centers,particle_num] = Circle_detect(image)
%function：检测输入图像的圆
%input:image:用于检测圆的图片彩色图片
%output:centers:检测到的圆的圆心坐标；particle_num:检测到图像上圆的总数

gray = rgb2gray(image);  %灰度化
bw = imbinarize(gray);   %二值化
bw_inverse = imcomplement(bw); %二值图像反转，即0变1,1变0
[centers,radius] = imfindcircles(bw_inverse,[5,20]); %检测图像上半径在5~20之间的圆，也可以是其它，看实际情况需要
particle_num = length(radius(:)); %统计圆的个数
%显示图像
figure;
imshow(bw_inverse);
hold on;
viscircles(centers,radius,'EdgeColor','b');

NumPerRadius.m（计算每个原子的每个单位圆内的原子个数）

function gr = NumPerRadius(centers,rc,dr)

[row,col] = size(centers);
num = round(rc/dr);
gr=zeros(num,1);

for i=1:(row-1)
    for j = i+1:row
        distance = sqrt((centers(i,1)-centers(j,1))^2 + (centers(i,2)-centers(j,2))^2);
        if distance <= rc
            lane = round(distance/dr);
            gr(lane) = gr(lane)+1;
        end
    end
end

RDF.m（计算径向分布）

function [index,percent] = rdf(height,width,gr,particle_num,dr)

[row,col] = size(gr);
percent = zeros(row,1);
gobal_rho = particle_num / (height * width);
for i=1:row
    temp = gr(i);
    temp = temp / particle_num;
    temp = temp / (pi*((i*dr)^2-((i-1)*dr)^2));
    percent(i) = temp / gobal_rho;
end

index = (1:row)/3.4;
percent = reshape(percent,1,row);
figure;
% plot(index,percent,'r');
xlim([0, 20]);
ylim([0, 20]);
values = spcrv([[index(1) index index(end)];[percent(1) percent percent(end)]],3);
plot(values(1,:),values(2,:), 'r');
% hold on;
% plot(index,percent,'go');

main.m（主函数）

close all;
clc;

image = imread('input2.png');
[height,width,channel] = size(image);
rc = width/2; %搜索圆的最大半径
n=15;  
dr = rc/n;
[centers,particle_num] = Circle_detect(image);
gr = NumPerRadius(centers,rc,dr);
[index,percent] = rdf(height,width,gr,particle_num,dr);

输入图片：

检测图像上原子的效果图：

径向分布函数效果图：

需要源代码的可从以下链接下载：
径向分布函数(RDF).zip