注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

东南隅

wantnon的blog

 
 
 

日志

 
 
 
 

正方形覆盖(改进)  

2009-08-18 17:16:55|  分类: 数模 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

在前面帖子“正方形覆盖”的基础上,作如下两点改进:
一,改用中心定位:
前帖中用于决定正方形方位的量(m,n),a中的(m,n)是正方形一角的坐标,
这样不对称,首先修改为(m,n)作为正方形中心的坐标,而正方形边长仍为r,a的意义也不变,
那么此时正方形的四个顶点坐标变为:
p1:(m-r/2*sin(a)+r/2*cos(a),n+r/2*cos(a)+r/2*sin(a))
p2:(m+r/2*cos(a)+r/2*sin(a),n+r/2*sin(a)-r/2*cos(a))
p3:(m+r/2*sin(a)-r/2*cos(a),n-r/2*cos(a)-r/2*sin(a))
p4:(m-r/2*cos(a)-r/2*sin(a),n-r/2*sin(a)+r/2*cos(a))
此时,要判断一点q:(x,y)是否在正方形内(包括边界),同样作与前面相同的旋转,得到q’(x’,y’),
则当且仅当满足:
-r/2<=x’<=r/2
-r/2<=y’<=r/2
时q在正方形内。
二,预先统计密度:
点在区域(0,0)-(1,1)上的分布是不均匀的,各处的疏密程度不同,而前帖中随机覆盖的方法,
将区域内各处等同地对待,这样是盲目的。
可改进为:事先统计一下各处的“密度”,在随机覆盖试验的过程中参考这些信息--即
在随机生成正方形中心点(m,n)时让它落在“密度”大处的机会多些,落在“密度”小处的机会少些。
以上所说“密度”为形象说法,具体如何表达仍是问题,下面采取一种最简单的方法:
将(0,0)-(1,1)等分成N*M个格子,统计每格内点数,然后当生成(m,n)是,让(m,n)落在各格内的概率
与其所含点数成正比。
各分区的编号规则如下图所示(其中M=N=5):
image 
可得出下面关系:
编号为k的分区对应的格子的下标为(I,J)=(k%M==0?floor(k/m):floor(k/m)+1,k%M==0?M:k%M),
坐标为((I-1)*w,(J-1)*h)到(I*w,J*h)。
另外一个问题就是如何使各格被选中的概率与其所含点数成正比:
设1*25的向量nper中已保存了各格所含点数(nper(k)为编号为k的格所含点数),
则显然p=nper/sum(nper)就是各格被选中的概率,于是我们有一个分布列:
1       2      3     …   25
p(1)   p(2)  p(3) …   p(25)
下面自定义函数mrand完成按指定分布列生成随机数:

function rs=mrand(p)
%按如下分布列生成随机数:
%1    2    3    ... length(p)
%p(1) p(2) p(3) ... p(end)
%注:p必须是分布,即必须满足sum(p)==1
v=rand();
s=0;
for i=1:length(p)
    s=s+p(i);
    if s>=v
        break;
    end
end
rs=i;
end

即每次首先用mrand生成随机数k,此即为本次正方形中心要击中的格子,
然后再在k号格子内去随机生成一个点(m,n):

m=rand()*w+(I(k)-1)*w;
n=rand()*h+(J(k)-1)*h;
这样就实现了“正方形中心落在各格的机会与其所含点数成正比”。
下面是完整程序:

%随机生成(0,0)-(1,1)的np个点
np=100;%点数
x=rand(1,np);
y=rand(1,np);
%将有效区域分成N*M个区,统计各区点数
M=5;
N=5;
w=1/M;
h=1/N;
nper=[];
I=[];
J=[];
for k=1:M*N
    %求第k区点数
    %求第k区脚标
    tv=rem(k,M);
    if tv==0
        I(end+1)=M;
        J(end+1)=floor(k/M);
    else
        I(end+1)=tv;
        J(end+1)=floor(k/M)+1;
    end
    %求第k区点数
    tcount=0;
    for i=1:length(x)
        if x(i)>=(I(end)-1)*w&&x(i)<=I(end)*w&&y(i)>=(J(end)-1)*h&&y(i)<=J(end)*h
            tcount=tcount+1;
        end
    end%得到tcount
    nper(end+1)=tcount; 
    %画第k区
    plot([(I(end)-1)*w,I(end)*w,I(end)*w,(I(end)-1)*w,(I(end)-1)*w],...
        [J(end)*h,J(end)*h,(J(end)-1)*h,(J(end)-1)*h,J(end)*h],'-g');
    hold on;
    text((I(end)-1)*w+w/6,J(end)*h,num2str(tcount));
end
plot(x,y,'o','Markersize',3);
r=0.1;%正方形边长
%mcount,mm,mn,ma记录最好的一次实验的count,m,n,a值
mcount=-1;
mm=[];
mn=[];
ma=[];
ntry=10000;%实验次数
p=nper/sum(nper);%被选中分区号的概率分布列
for turn=1:ntry
    k=mrand(p);%本次正方形中心落在第k区
    %生成第k区的随机点(m,n)
    m=rand()*w+(I(k)-1)*w;
    n=rand()*h+(J(k)-1)*h;
    a=rand()*pi/2;%随机转角
    count=0;%本次实验覆盖的点数
    for i=1:length(x)
        %判断x(i),y(i)是否被覆盖
        x_(i)=(x(i)-m)*sin(a)-(y(i)-n)*cos(a);
        y_(i)=(x(i)-m)*cos(a)+(y(i)-n)*sin(a);
        if x_(i)>=-r/2&&x_(i)<=r/2&&y_(i)>=-r/2&&y_(i)<=r/2
            count=count+1;
        end
    end
    %记录最优值
    if count>mcount
        mcount=count;
        mm=m;
        mn=n;
        ma=a;
    end
end
%设置标题
title({'',...%第一行空出为了美观
    ['随机生成(0,0)-(1,1)上',num2str(np),'个点,用边长',num2str(r),'的正方形去覆盖'],...
    ['区域分成',num2str(N),'*',num2str(M),'个格子,各格上数字为其所含点数'],...
    ['所含点数越多,正方形中心落在上面的几率越大'],...
    ['共进行了',num2str(ntry),'次实验'],...
    ['最好的一次覆盖了',num2str(mcount),'个点']},'fontsize',8);
%最优一次实验中矩形的四个顶点坐标p1,p2,p3,p4
p1=[mm-r/2*sin(ma)+r/2*cos(ma),mn+r/2*cos(ma)+r/2*sin(ma)];
p2=[mm+r/2*cos(ma)+r/2*sin(ma),mn+r/2*sin(ma)-r/2*cos(ma)];
p3=[mm+r/2*sin(ma)-r/2*cos(ma),mn-r/2*cos(ma)-r/2*sin(ma)];
p4=[mm-r/2*cos(ma)-r/2*sin(ma),mn-r/2*sin(ma)+r/2*cos(ma)];
%画正方形
plot([p1(1),p2(1),p3(1),p4(1),p1(1)],[p1(2),p2(2),p3(2),p4(2),p1(2)],'-r');
axis equal;
axis([-0.1,1.1,-0.1,1.1]);

运行结果:

image

  评论这张
 
阅读(99)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017