计算元素区域邻居元素块的个数
问问大虾们这样一个问题，比如一个矩阵里面有不同的元素，但是相同的元素都连通在一起（上下左右连通），怎么计算矩阵里面元素1的周围有几个不同的元素区（相同元素的块）？
比如附图中的元素1周围就有2,3,4,5,6共5个元素区，需要计算的结果就是5个，实际问题中矩阵元素块很多，数不清楚，如何用程序实现呢？
希望路过的虾米们指点一下！ 本帖最后由 scott198510 于 2010-12-7 11:36 编辑
clear all;clc;
S=[2 2 2 3 3 4 4
2 2 1 1 1 4 4
6 1 1 1 1 1 4
6 1 1 1 1 1 5
6 6 1 1 1 5 5;
6 6 6 6 5 5 5;
6 6 6 6 5 5 5];
N=7;
Sd=zeros(N+2); % 把状态矩阵扩大至(N+2)×(N+2)
Sd(2:N+1,1)=S(:,N); % 把Sd处理为周期性边界条件
Sd(2:N+1,N+2)=S(:,1); % 把Sd处理为周期性边界条件
Sd(1,2:N+1)=S(N,:); % 把Sd处理为周期性边界条件
Sd(N+2,2:N+1)=S(1,:); % 把Sd处理为周期性边界条件
Sd(1,1)=S(N,N); % 把Sd处理为周期性边界条件
Sd(N+2,N+2)=S(1,1); % 把Sd处理为周期性边界条件
Sd(1,N+2)=S(N,1); % 把Sd处理为周期性边界条件
Sd(N+2,1)=S(1,N); % 把Sd处理为周期性边界条件
Sd(2:end-1, 2:end-1)=S% 把Sd处理为周期性边界条件
y=Sd(2:end-1,2:end-1);
x=unique(y);
s=cell(length(x),2);
for i=1:length(x)
=find(y==x(i));
s{i,1}=x(i);
a=unique(Sd(m,n+1));
b=unique(Sd(m+2,n+1));
c=unique(Sd(m+1,n));
d=unique(Sd(m+1,n+2));
f=horzcat(a',b',c',d');
f(f==x(i))=[ ];
s{i,2}=unique(f);
end%块的个数
C=s
用这种方法得到的是上面矩阵里面每个元素块的邻居元素块的个数，完全准确
但是为什么计算下面这个矩阵的时候，元素8组成的块，周围明明没有邻居1，统计却有1呢？
clear all;clc;
S=[2 2 2 3 3 4 499
2 2 1 1 1 4 499
6 1 1 1 1 1 499
6 1 1 1 1 1 598
6 6 1 1 1 5 588
6 6 6 6 5 5 588
6 6 6 6 5 5 588
7 7 7 7 8 8 888
7 7 7 7 8 8 88 8 ]
N=9;
Sd=zeros(N+2); % 把状态矩阵扩大至(N+2)×(N+2)
Sd(2:N+1,1)=S(:,N); % 把Sd处理为周期性边界条件
Sd(2:N+1,N+2)=S(:,1); % 把Sd处理为周期性边界条件
Sd(1,2:N+1)=S(N,:); % 把Sd处理为周期性边界条件
Sd(N+2,2:N+1)=S(1,:); % 把Sd处理为周期性边界条件
Sd(1,1)=S(N,N); % 把Sd处理为周期性边界条件
Sd(N+2,N+2)=S(1,1); % 把Sd处理为周期性边界条件
Sd(1,N+2)=S(N,1); % 把Sd处理为周期性边界条件
Sd(N+2,1)=S(1,N); % 把Sd处理为周期性边界条件
Sd(2:end-1, 2:end-1)=S% 把Sd处理为周期性边界条件
y=Sd(2:end-1,2:end-1);
x=unique(y);
s=cell(length(x),2);
for i=1:length(x)
=find(y==x(i));
s{i,1}=x(i);
a=unique(Sd(m,n+1));
b=unique(Sd(m+2,n+1));
c=unique(Sd(m+1,n));
d=unique(Sd(m+1,n+2));
f=horzcat(a',b',c',d');
f(f==x(i))=[];
s{i,2}=unique(f);
end%块的个数
C=s
本帖最后由 feynmand 于 2010-12-7 15:14 编辑
你的程序没有兴趣找错误了，这里给出我的解答。希望你能够看懂为什么我乘以了个足够大的数100，这是解决问题的关键。
%%
clear
clc
S=[2 2 2 3 3 4 499
2 2 1 1 1 4 499
6 1 1 1 1 1 499
6 1 1 1 1 1 598
6 6 1 1 1 5 588
6 6 6 6 5 5 588
6 6 6 6 5 5 588
7 7 7 7 8 8 888
7 7 7 7 8 8 88 8 ];
N=length(S);
Sd=zeros(N+2); % 把状态矩阵扩大至(N+2)×(N+2)
Sd(2:N+1,1)=S(:,N); % 把Sd处理为周期性边界条件
Sd(2:N+1,N+2)=S(:,1); % 把Sd处理为周期性边界条件
Sd(1,2:N+1)=S(N,:); % 把Sd处理为周期性边界条件
Sd(N+2,2:N+1)=S(1,:); % 把Sd处理为周期性边界条件
Sd(1,1)=S(N,N); % 把Sd处理为周期性边界条件
Sd(N+2,N+2)=S(1,1); % 把Sd处理为周期性边界条件
Sd(1,N+2)=S(N,1); % 把Sd处理为周期性边界条件
Sd(N+2,1)=S(1,N); % 把Sd处理为周期性边界条件
Sd(2:end-1, 2:end-1)=S;% 把Sd处理为周期性边界条件
x=unique(S);
s=cell(length(x),2);
for n=1:length(x)
s{n,1}=x(n);
S0=Sd;
S0(S0==x(n))=x(n)*100;
S1=S0(2:end,:); % 横向
S2=S0(:,2:end); % 纵向
S3=(S0(1:end-1,:)+S1)/2;
S4=(S0(:,1:end-1)+S2)/2;
S5=;
S6=S5(S5>x(n)*100/2);
xx=unique(S6);
s{n,2}=xx(1:end-1)*2-xx(end);
end
版主的算法确实很高明，我刚才用斑竹的算法计算了一个400X400的矩阵，结果有问题，里面统计的邻居居然有负数，不过，我把版主强调的那个100改成1000，问题就没有了，结果很准确，我好好研究这个100的意义，貌似，计算的矩阵维数越大，那个100要改大才行 足够大的数100
注意这个“足够大”是关键，为了避免出错，这个数尽量大，最好比输入矩阵的最大值大一个数量级。 分享一下：
function C=demo_cal
clear all;
S=[2 2 2 3 3 4 4 9 9
2 2 1 1 1 4 4 9 9
6 1 1 1 1 1 4 9 9
6 1 1 1 1 1 5 9 8
6 6 1 1 1 5 5 8 8
6 6 6 6 5 5 5 8 8
6 6 6 6 5 5 5 8 8
7 7 7 7 8 8 8 8 8
7 7 7 7 8 8 8 8 8 ]
N=9;
Sd=zeros(N+2); % 把状态矩阵扩大至(N+2)×(N+2)
Sd(2:N+1,1)=S(:,N); % 把Sd处理为周期性边界条件
Sd(2:N+1,N+2)=S(:,1); % 把Sd处理为周期性边界条件
Sd(1,2:N+1)=S(N,:); % 把Sd处理为周期性边界条件
Sd(N+2,2:N+1)=S(1,:); % 把Sd处理为周期性边界条件
Sd(1,1)=S(N,N); % 把Sd处理为周期性边界条件
Sd(N+2,N+2)=S(1,1); % 把Sd处理为周期性边界条件
Sd(1,N+2)=S(N,1); % 把Sd处理为周期性边界条件
Sd(N+2,1)=S(1,N); % 把Sd处理为周期性边界条件
Sd(2:end-1, 2:end-1)=S % 把Sd处理为周期性边界条件
y=Sd(2:end-1,2:end-1);
x=unique(y);
s=cell(length(x),2);
for i=1:length(x)
=find(y==x(i));
s{i,1}=x(i);
a=calculate(Sd,m,n+1);
b=calculate(Sd,m+2,n+1);
c=calculate(Sd,m+1,n);
d=calculate(Sd,m+1,n+2);
f=horzcat(a,b,c,d);
i
f(f==x(i))=[ ];
s{i,2}=unique(f);
end % 块的个数
C=s
end
function ff=calculate(A,m,n)
for ii=1:length(m)
ff(ii)=A(m(ii),n(ii));
end
ff=unique(ff);
end 5# feynmand
版主给 liuyalong008加分吧，他这个算法就是我最上面说的那个计算8的邻居有问题的改进版，我试过了，两种算法一样好，都是妙哉的。。完全达到要求 看大家都在写，我也写一个，供参考clear
clc
S=[2 2 2 3 3 4 499
2 2 1 1 1 4 499
6 1 1 1 1 1 499
6 1 1 1 1 1 598
6 6 1 1 1 5 588
6 6 6 6 5 5 588
6 6 6 6 5 5 588
7 7 7 7 8 8 888
7 7 7 7 8 8 88 8 ];
S0=feval(@(x)x(size(S,1):2*size(S,1)+1,size(S,2):2*size(S,2)+1),repmat(S,3,3));
Sn=cell2mat(cellfun(@(x)x(:),{S,S0(1:end-2,2:end-1),S0(3:end,2:end-1),...
S0(2:end-1,1:end-2),S0(2:end-1,3:end)},'UniformOutput',false));
f=arrayfun(@(x)setdiff((unique(Sn(S==x,:))),x),1:max(S(:)),'UniformOutput',false)'; 本帖最后由 bainhome 于 2010-12-7 22:28 编辑
这好像是8邻域算法吧？都是猛人！qibbxxt最近arrayfun用上瘾了，看见别人代码里的“for”，眼睛里就冒着arrayfun的火星儿，呵呵。 8# qibbxxt
高手！！！！ 本帖最后由 scott198510 于 2011-10-29 19:37 编辑
5# feynmand
参照斑竹的求法 求了一下三维空间的邻居问题：
即各个元素块在三维空间的邻居元素块个数（三维空间的面数），二维求的是边数：
clear ;clc ;close all;
loadS
N=length(S);
% ===============把状态矩阵扩大至(N+2)×(N+2)===============================
Sd=zeros(N+2,N+2,N+2); % 把状态矩阵扩大至(N+2)×(N+2)
Sd(2:N+1,1,2:N+1)=S(1:N,N,1:N);% 把Sd处理为周期性边界条件 右移左
Sd(2:N+1,N+2,2:N+1)=S(1:N,1,1:N); % 把Sd处理为周期性边界条件 左移右
Sd(2:N+1,2:N+1,1)=S(1:N,1:N,N);% 把Sd处理为周期性边界条件上移下
Sd(2:N+1,2:N+1,N+2)=S(1:N,1:N,1);% 把Sd处理为周期性边界条件 下移上
Sd(1,2:N+1,2:N+1)=S(N,1:N,1:N); % 把Sd处理为周期性边界条件 前移后
Sd(N+2,2:N+1,2:N+1)=S(1,1:N,1:N); % 把Sd处理为周期性边界条件 后移前
Sd(2:end-1, 2:end-1,2:end-1)= S; % 把Sd处理为周期性边界条件 前移后
% ===============把状态矩阵扩大至(N+2)×(N+2)===============================
x=unique(S);
s=cell(length(x),2);
for n=1:length(x)
s{n,1}=x(n);
S0=Sd;
S0(S0==x(n))=x(n)*50000;
S1=S0(2:end,:,:); % 水平方向
S2=S0(:,2:end,:); % 垂直方向
S3=S0(:,:,2:end); % 竖直方向
S4=(S0(1:end-1,:,:)+S1)/2;
S5=(S0(:,1:end-1,:)+S2)/2;
S6=(S0(:,:,1:end-1)+S3)/2;
S7=;
S8=S7(S7>x(n)*50000/2);
xx=unique(S8);
s{n,2}=xx(1:end-1)*2-xx(end);
end
C=cell(length(x),3);
C(:,1:2)=s;
for i=1:length(x)
C{i,3}=length(C{i,2});
end
AB=(cell2mat(C(:,3)))';
table=tabulate (AB);
bar(table(:,2))
set(gca,'xlim',)
set(gcf,'color','w');
可是按道理求出来的应该是14面最多，结果却是17面的最多？ 11# scott198510 12# scott198510
你把你的数据矩阵贴出来 本帖最后由 scott198510 于 2011-3-3 10:33 编辑
13# qibbxxt
这个矩阵数据很多 530KB 怎么上传不上去 把矩阵维数改了 14# scott198510
在我原来的程序上面稍加修改，时间关系，没有用cell进行优化，你试一试吧
clear
clc
load S;
siz=arrayfun(@(x)x:2*x+1,size(S),'UniformOutput',0);
S0=feval(@(x)x(siz{:}),repmat(S,3+zeros(1,3)));
Sn=cell2mat(cellfun(@(x)x(:),{S,S0(1:end-2,2:end-1,2:end-1),S0(3:end,2:end-1,2:end-1),...
S0(2:end-1,1:end-2,2:end-1),S0(2:end-1,3:end,2:end-1),...
S0(2:end-1,2:end-1,1:end-2),S0(2:end-1,2:end-1,3:end)},...
'UniformOutput',false));
f=arrayfun(@(x)setdiff((unique(Sn(S==x,:))),x),1:max(S(:)),'UniformOutput',false)';
二维矩阵代码。
三维矩阵情况的代码在此基础上拓展即可。（15#附件下不了）
clear all;
clc;
S = [
2 2 2 3 3 4 4 9 9
2 2 1 1 1 4 4 9 9
6 1 1 1 1 1 4 9 9
6 1 1 1 1 1 5 9 8
6 6 1 1 1 5 5 8 8
6 6 6 6 5 5 5 8 8
6 6 6 6 5 5 5 8 8
7 7 7 7 8 8 8 8 8
7 7 7 7 8 8 8 8 8
];
[ mz, nz ] = size(S);
[ dm, dn ] = meshgrid([ -1, 0, 1 ], [ -1, 0, 1 ]);
TheNo = 5; % take 5 for example.
[ m0, n0 ] = find(S == TheNo);
m = kron(m0,ones(3)) + repmat(dm,);
n = kron(n0,ones(3)) + repmat(dn,);
m(m<1) = 1;
m(m>mz) = mz;
n(n<1) = 1;
n(n>nz) = nz;
ind = sub2ind([ mz, nz ], m, n);
mySet = unique(S(ind));
setdiff(mySet( : ), TheNo) % take 5 for example.
本帖最后由 scott198510 于 2011-10-29 19:36 编辑
15# qibbxxt
采用斑竹的方法和我上面套用feynmand 虾米的方法都求出来了，斑竹的方法求出来的速度要快一点，套用feynmand 斑竹方法求的，不知道为什么要运行差不多半小时才出来结果。我对比了一下，发现统计的一样，说明是我的矩阵问题，面数最多的不是14面，都是17面最多
clear ;clc ;close all;
loadS
N=length(S);
siz=arrayfun(@(x)x:2*x+1,size(S),'UniformOutput',0);
S0=feval(@(x)x(siz{:}),repmat(S,3+zeros(1,3)));
Sn=cell2mat(cellfun(@(x)x(:),{S,S0(1:end-2,2:end-1,2:end-1),S0(3:end,2:end-1,2:end-1),S0(2:end-1,1:end-2,2:end-1),S0(2:end-1,3:end,2:end-1),S0(2:end-1,2:end-1,1:end-2),S0(2:end-1,2:end-1,3:end)},'UniformOutput',false));
f=arrayfun(@(x)setdiff((unique(Sn(S==x,:))),x),1:max(S(:)),'UniformOutput',false)';
for i=1:length(unique(S));
M(i) = length(cell2mat(f((i))));
end
AB=M ;
table=tabulate (AB);
bar(table(:,2))
set(gca,'xlim',)
set(gcf,'color','w');
三维问题代码
clear all;
clc;
load S
tic
% S = [ %改为三位矩阵数据
% ...
% ];
[ Xind_max, Yind_max, Zind_max ] = size(S);
[ dXind, dYind, dZind ] = meshgrid([ -1, 0, 1 ], [ -1, 0, 1 ], [ -1, 0, 1 ]);
S_elements = unique(S);
Result = zeros(1,2); % 元素及个数(面数)
for kk = 1 : length(S_elements)
TheNo = S_elements(kk); % take 5 for example.
ind = find(S == TheNo); % take 5 for example.
[ Xind0, Yind0, Zind0 ] = ind2sub([ Xind_max, Yind_max, Zind_max ], ind);
mySet = [ ];
for k = 1 : length(Xind0)
Xind = dXind + Xind0(k);
Xind(Xind<1) = 1;
Xind(Xind>Xind_max) = Xind_max;
Yind = dYind + Yind0(k);
Yind(Yind<1) = 1;
Yind(Yind>Yind_max) = Yind_max;
Zind = dZind + Zind0(k);
Zind(Zind<1) = 1;
Zind(Zind>Zind_max) = Zind_max;
ind = unique(sub2ind([ Xind_max, Yind_max, Zind_max ], Xind, Yind, Zind));
mySet = union(mySet, unique(S(ind)));
end
Result(kk,:) = [ S_elements(kk), length(setdiff(mySet( : ), TheNo)) ]; % take 5 for example.
end
Result
toc
...
7991 11
7992 23
7993 24
7994 25
7995 20
7996 18
7997 16
7998 22
7999 28
8000 20
Elapsed time is 520.034481 seconds.
，
本帖最后由 scott198510 于 2011-10-31 09:59 编辑qibbxxt 发表于 2011-3-3 14:19 http://forum.simwe.com/static/image/common/back.gif
14# scott198510
在我原来的程序上面稍加修改，时间关系，没有用cell进行优化，你试一试吧
...
套用斑竹的方法求三维的面数的时候，比如某个元素50在三维空间内的邻居里面元素9只出现了2次，这时候按照原来的方法求的面数也算进去了9，
实际上当他的某个邻居出现的次数很少的时候，需要把这个邻居剔除，也就是说计算的面数少计算一次，上述只是一个例子，
现在 比如设置一个阀值Fa=8, 当某个元素出现的次数少于8时，就不计算它的面数，
且计算凡是个数少于8的邻居的时候的面数也少计算一次；相应的程序改怎么改呢；
下面附上矩阵 :S750 ,如果按照原来的方法计算面数，得到的分布图如下图，一面的占了很大比例，这显然是不正确的：
clear ;clc ;close all;
loadS
N=length(S);
siz=arrayfun(@(x)x:2*x+1,size(S),'UniformOutput',0);
S0=feval(@(x)x(siz{:}),repmat(S,3+zeros(1,3)));
Sn=cell2mat(cellfun(@(x)x(:),{S,S0(1:end-2,2:end-1,2:end-1),S0(3:end,2:end-1,2:end-1),S0(2:end-1,1:end-2,2:end-1),S0(2:end-1,3:end,2:end-1),S0(2:end-1,2:end-1,1:end-2),S0(2:end-1,2:end-1,3:end)},'UniformOutput',false));
f=arrayfun(@(x)setdiff((unique(Sn(S==x,:))),x),1:max(S(:)),'UniformOutput',false)';
for i=1:length(unique(S));
M(i) = length(cell2mat(f((i))));
end
AB=M ;
table=tabulate (AB);
bar(table(:,2))
set(gca,'xlim',)
set(gcf,'color','w');
xlabel('三维晶粒的面数','Fontsize',18);
ylabel('出现次数','Fontsize',18);
title('统计分布图','Fontsize',18);
页:
[1]