【问题标题】:Estimate the amount of the area of the hemisphere’s surface which can be seen by the eye?估计肉眼可以看到的半球表面面积?
【发布时间】:2016-10-22 22:00:16
【问题描述】:

我们假设 3D 空间中有一个半球和三个三角形。半球底部的中心点用 C 表示。半球底部的半径用 R 表示。半球底部的法向量用 n 表示。

第一个三角形由三个点 V1、V2 和 V3 给出。第二个三角形由三个点 V4、V5 和 V6 给出。第三个三角形由三个点 V7、V8 和 V9 给出。点 V1、V2、……、V9 的位置是任意的。现在,我们进一步假设眼睛位于点 E。请注意,三角形可能会阻挡从眼睛到半球表面的视线;因此,肉眼可能看不到半球表面的某些区域。

请开发一种方法来估计肉眼可以看到的半球表面面积? 这是矩形而不是半球的代码:

function r = month_1()
P1 = [0.918559, 0.750000, -0.918559];
P2 = [0.653394, 0.649519, 1.183724];
P3 = [-0.918559, -0.750000, 0.918559];
P4 = [-0.653394, -0.649519, -1.183724];

V1 = [0.989609, -1.125000, 0.071051];
V2 = [1.377838, -0.808013, -0.317178];
V3 = [1.265766, -0.850481, 0.571351];

V4 = [0.989609, -1.125000, 0.071051];
V5 = [1.265766, -0.850481, 0.571351];
V6 = [0.601381, -1.441987, 0.459279];

V7 = [0.989609, -1.125000, 0.071051];
V8 = [1.377838, -0.808013, -0.317178];
V9 = [0.713453, -1.399519, -0.429250];

E = [1.714054, -1.948557, 0.123064];

C = [0,1,0];
Radius = 2;
n = [0,1,0]; 

%hold on
patches.vertices(1,:)= P1;
patches.vertices(2,:)= P2;
patches.vertices(3,:)= P3;
patches.vertices(4,:)= P4;

patches.vertices(5,:)= V1;
patches.vertices(6,:)= V2;
patches.vertices(7,:)= V3;
patches.vertices(8,:)= V4;
patches.vertices(9,:)= V5;
patches.vertices(10,:)= V6;
patches.vertices(11,:)= V7;
patches.vertices(12,:)= V8;
patches.vertices(13,:)= V9;

patches.faces(1,:)= [5,6,7];
patches.faces(2,:)= [8,9,10];
patches.faces(3,:)= [11,12,13];
patches.faces(4,:)= [1,2,3];
patches.faces(5,:)= [1,3,4];

patches.facevertexcdata = 1;
patch(patches);
shading faceted; view (3);
% dispatch([1,1,1])
hold on

Num = 20;
Sum = 0;
%Srec = norm(cross(P1 - P4, P3 - P4));
for i = 1:Num
    x1 = rand;
    x2 = rand;
    v1 = P1-P4;
    v2 = P3-P4;
    Samp = P4+v1*x1+v2*x2;
    ABC = fun_f(E, Samp, V1,V2,V3)*fun_f(E,Samp, V4, V5, V6)*fun_f(E,Samp, V7,V8,V9);
    Sum = Sum + ABC;
   % if ABC ==1
   %      plot3(Samp(1), Samp(2), Samp(3),'r +'), hold on
   % else
   %      plot3(Samp(1), Samp(2), Samp(3),'b +'), hold on
   % end
%............................
[x, y, z]= sphere;
patches = surf2patch(x,y,z,z);
patches.vertices(:,3) = abs(patches.vertices(:,3));
patches.facevertexcdata = 1;
patch(patches)
shading faceted; view(3)
daspect([1 1 1])
%............................
end

%r = Sum/Num;
%view(31, 15)
%end


r = Sum/Num*norm(cross (P1-P4,P3-P4));
disp(sprintf('the integration is: %.5f',r));
disp(sprintf('the accurate result is: %.5f',norm(cross(P1-P4,P3-P4)/4)));



function res = fun_f(LineB, LineE, V1, V2, V3)
O = LineB;
Len = norm(LineE-LineB);
v = (LineE-LineB)/Len;
N = cross(V2-V1, V3-V1)/norm(cross(V2-V1, V3-V1));
if dot(N,v)~=0
    tp = dot(N, V1-O)/ dot(N,v);
   % if tp >=0 && tp <= (1:3);
    P = LineB+tp*v(1:3);

    A = V1 - P;
    B = V2 - P;
    Stri1 = norm(cross(A,B))/2;

    A = V1 - P;
    B = V3 - P;
    Stri2 = norm(cross(A,B))/2;

    A = V3 - P;
    B = V2 - P;
    Stri3 = norm(cross(A,B))/2;

    A = V1 - V2;
    B = V3 - V2;
    Stotal = norm(cross(A,B))/2;

    if (Stri1 + Stri2 + Stri3)> (Stotal + 1e-8)
        res = 1;
    else
        res = 0;
    end
 else
    res =1;
end
end
end

【问题讨论】:

  • 我有一个想法,但我只能给你数学,而不是 MatLab。你想要吗?
  • @willywonkadailyblah 是的,我想要。请在数学中给出想法,然后我将其转换为matlab。

标签: matlab simulation matlab-figure matlab-cvst geometry-surface


【解决方案1】:

选取以 为中心的一小部分表面积,尺寸为。面积元素由下式给出

这个想法是在球体上循环这些元素;计算处元素的中心点,并计算该点与相机之间的线段是否与三角形相交。更多信息:https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm

现在我们需要找到点;这意味着在整个半球上将 增加但是这会使采样分辨率不均匀 - 因子会使半球顶点附近的元素比其边缘大得多。

改为:

  • 设置一个固定的环数N循环,即的迭代次数。

  • 设置最小迭代区域M的迭代次数由

    给出

    其中 是环在 处的总面积:

  • 当然,上面的增量由下式给出

  • 遍历所有环,注意 给出每个环的中间 行(所以从 开始);由于对称性,同样的问题不必适用于。对于每个 上的每个环,并按照上述方法进行线交叉测试。

上述方法在小处减少了区域采样分辨率的偏差量。

更好的方法是斐波那契格,但它们更复杂。看这篇论文:http://geonaut.eu/published/016_Fibonacci_Lattice.pdf

【讨论】:

    【解决方案2】:
    clc
    %% Declaration of Initial Variables
    C = [0.918559, 0.750000, -0.918559];
    R = 10;
    n = [1, 2, 1.5];
    V1 = [0.989609, -1.125000, 0.071051];
    V2 = [1.377838, -0.808013, -0.317178];
    V3 = [1.265766, -0.850481, 0.571351];
    V4 = [0.989609, -1.125000, 0.071051];
    V5 = [1.265766, -0.850481, 0.571351];
    V6 = [0.601381, -1.441987, 0.459279];
    V7 = [0.989609, -1.125000, 0.071051];
    V8 = [1.377838, -0.808013, -0.317178];
    V9 = [0.713453, -1.399519, -0.429250];
    E = [1.714054, -1.948557, 0.123064];
    Num = 10000;
    count1 = 0; count2 = 0; count3 = 0;
    
    %% Calculating the triangles Normal and Area
    N1 = cross((V2-V1),(V3-V1));
    N2 = cross((V5-V4),(V6-V4));
    N3 = cross((V8-V7),(V9-V7));
    Area1 = norm(N1)/2;
    Area2 = norm(N2)/2;
    Area3 = norm(N3)/2;
    
    %% Plotting the triangle
    patch([V1(1),V2(1),V3(1),V1(1)],[V1(2),V2(2),V3(2),V1(2)],[V1(3),V2(3),V3(3),V1(3)], 'green');
    hold on
    patch([V4(1),V5(1),V6(1),V4(1)],[V4(2),V5(2),V6(2),V4(2)],[V4(3),V5(3),V6(3),V4(3)],'green');
    hold on
    patch([V7(1),V8(1),V9(1),V7(1)],[V7(2),V8(2),V9(2),V7(2)],[V7(3),V8(3),V9(3),V7(3)], 'green');
    plot3(E(1),E(2),E(3),'ro')
    hold on
    %% The Loop Section
    for i=1:Num
        x1 = rand;
        x2 = rand;
        Px = R*sqrt(x1*(2-x1))*cos(2*pi*x2)+C(1);
        Py = R*sqrt(x1*(2-x1))*sin(2*pi*x2)+C(2);
        Pz = R*(1 - x1)+C(3);
        z = [0,0,1];
        if norm(cross(n,z)) ~= 0
            v = cross(n,z)/norm(cross(n,z));
            u = cross(v,n)/norm(cross(v,n));
            w = n/norm(n);
        else
            u = (dot(n,z)/abs(dot(n,z))).*[1,0,0];
            v = (dot(n,z)/abs(dot(n,z))).*[0,1,0];
            w = (dot(n,z)/abs(dot(n,z))).*[0,0,1];
        end
        M = [u(1),u(2),u(3),0;v(1),v(2),v(3),0;w(1),w(2),w(3),0;0,0,0,1]*...
            [1,0,0,-C(1);0,1,0,-C(2);0,0,1,-C(3);0,0,0,1];
        L = [Px,Py,Pz,1]*M;
        Points = [L(1),L(2),L(3)];
        Len = norm(E - Points);
        plot3(L(1),L(2),L(3),'b.'),hold on
        dv = (E - Points)/norm(E - Points);
    
        %% Triangle 1
        tp1 = dot(N1,(V1-Points))/dot(N1,dv);
        if tp1>=0 && tp1<=Len
            R1 = Points + tp1*dv;
            A1 = norm(cross((V1-R1),(V2-R1)))/2;
            A2 = norm(cross((V1-R1),(V3-R1)))/2;
            A3 = norm(cross((V2-R1),(V3-R1)))/2;
            if (A1+A2+A3) <= Area1
                count1 = count1 + 1;
                plot3(L(1),L(2),L(3),'r*')
    
            end
        end
        %% Triangle 2
        tp2 = dot(N2,(V4-Points))/dot(N2,dv);
        if tp2>=0 && tp2<=Len
            R2 = Points + tp2*dv;
            A4 = norm(cross((V4-R2),(V5-R2)))/2;
            A5 = norm(cross((V4-R2),(V6-R2)))/2;
            A6 = norm(cross((V5-R2),(V6-R2)))/2;
            if (A4+A5+A6) <= Area2
                count2 = count2 + 1;
                plot3(L(1),L(2),L(3),'r*')
            end
        end
    
        %% Triangle 3
        tp3 = dot(N3,(V7-Points))/dot(N3,dv);
        if tp3>=0 && tp3<=Len
            R3 = Points + tp3*dv;
            A7 = norm(cross((V7-R3),(V8-R3)))/2;
            A8 = norm(cross((V7-R3),(V9-R3)))/2;
            A9 = norm(cross((V8-R3),(V9-R3)))/2;
            if  (A7+A8+A9) <= Area3
                count3 = count3 + 1;
                plot3(L(1),L(2),L(3),'r*')
            end
        end
    
    end
    %% Final Solution
    AreaofHemi = 2*pi*R^2;
    Totalcount = count1 + count2 + count3;
    Areaseen=((Num-Totalcount)/Num)*AreaofHemi;
    disp(fprintf('AreaofHemi %f, AreaSeen %f ',AreaofHemi,Areaseen))
    

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 2013-05-06
      • 2010-11-23
      • 2014-06-15
      • 2011-06-08
      • 2019-04-03
      • 2020-02-24
      相关资源
      最近更新 更多