如何通过已知参数产生椭圆矩阵

ll710915

如何产生一椭圆矩阵，在椭圆内部值为1，在外部为0，椭圆参数知道（长轴,短轴，,角度），谢谢！

[ 本帖最后由 eight 于 2007-10-21 23:23 编辑 ]

花如月

如果没有好办法，就直接赋值吧

花如月

试了下，最笨的方法，也是最直接的方法

1. clear,clc;
   
2. N=111;A=ones(N);o=floor((N+1)/2);
   
3. for i=1:N
   
4.   for j=1:N
   
5.     if ((i-o)^2/9+(j-o)^2/16)>10^2
   
6.         A(i,j)=0;
   
7.     end
   
8.    end
   
9. end
   
10. imshow(A)

复制代码

[ 本帖最后由 花如月 于 2007-10-21 18:09 编辑 ]

ll710915

可是我还有个参数，带角度的，希望不用循环，实现的话万分感谢

ll710915

如果采用花如月 的方法代码可以如下修改：
function A= EliMatrix(N,R1,R2,r,alpha)

A=ones(N);o=floor((N+1)/2);
for i=1:N
  for j=1:N
x=i-o;
y=j-o;
c1=cos(alpha)^2/R1^2+sin(alpha)^2/R2^2;
c2=sin(2*alpha)/R1^2-sin(2*alpha)/R2^2;
c3=sin(alpha)^2/R1^2+cos(alpha)^2/R2^2;
if (c1*x^2+c2*x*y+c3*y^2)>r^2
        A(i,j)=0;
    end
   end
end
imshow(A)
在双核处理下时间为0.35秒，时间太长，这只是一个mask而已，所以寻求不是循环的，另外还要解决矩阵与椭圆外切的关系，进行中。。。

花如月

一步一步来吧，欢迎分享你的解决办法

eight

原帖由 ll710915 于 2007-10-22 14:09 发表

                               
登录/注册后可看大图

如果采用花如月 的方法代码可以如下修改：
function A= EliMatrix(N,R1,R2,r,alpha)

A=ones(N);o=floor((N+1)/2);
for i=1:N
  for j=1:N
x=i-o;
y=j-o;
c1=cos(alpha)^2/R1^2+sin(alpha)^2/R2^2;
c2 ...

自行按照“矢量化编程技巧”的帖子（置顶帖可以找到）修改一下即可去掉所有循环

花如月

看来我太愚钝了，怎么我看来了2遍矢量化编程。
还是没有找到去掉循环的办法

eight

原帖由 花如月 于 2007-10-23 15:26 发表

                               
登录/注册后可看大图

看来我太愚钝了，怎么我看来了2遍矢量化编程。
还是没有找到去掉循环的办法

c1,c2,c3值是固定的，条件判断 if 中 r 也是固定的，只有 x、y 是可变的，A 随着 x、y 的改变有不同的取值，这就可以摆脱循环

花如月

依然没有看出去掉所有循环的办法:@L

eight

原帖由 花如月 于 2007-10-23 15:53 发表

                               
登录/注册后可看大图

依然没有看出去掉所有循环的办法:@L

1. 
   
2. clear,clc;
   
3. N=1024;C=ones(N);o=floor((N+1)/2);
   
4. X=repmat([1:N]',1,N);
   
5. Y=repmat(1:N,N,1);
   
6. C(((X-o).^2/9+(Y-o).^2/16)>10^2) = 0;
   
7. imshow(C)

复制代码

两个 repmat 语句可以用一个 meshgrid 语句代替


速度比较：


[ 本帖最后由 eight 于 2007-10-23 16:58 编辑 ]

花如月

忽然想起来，曾经见到过这种用法。当时不太理解，所以不会用。eight总是能给人惊喜，再次感谢！

ll710915

谢谢！

ll710915

function A= EliMatrix(R1,R2,r,alpha)
N=2*max(R1,R2)*r;
A=ones(N);o=floor((N+1)/2);
c1=cos(alpha)^2/R1^2+sin(alpha)^2/R2^2;
c2=sin(2*alpha)/R1^2-sin(2*alpha)/R2^2;
c3=sin(alpha)^2/R1^2+cos(alpha)^2/R2^2;
X=repmat([1:N]',1,N);
Y=repmat(1:N,N,1);
A((c1*(X-o).^2+c2.*(X-o).*(Y-o)+c3*(Y-o).^2)>r^2)=0;
imshow(A)