Estimate the amount of the area of the hemisphere’s surface which can be seen by the eye?

Question

We suppose that there are one hemisphere and three triangles in a 3D space. The center point of the hemisphere’s bottom is denoted by C. The radius of the hemisphere’s bottom is represented by R. The normal vector to hemisphere’s bottom is denoted by n.

The first triangle is given by the three points V1, V2 and V3. The second triangle is given by the three points V4, V5 and V6. The third triangle is given by the three points V7, V8 and V9. The positions of the points V1, V2, …, V9 are arbitrary. Now, we further assume that an eye is located at the point E. Note that the triangles may block the line of sight from the eye to the hemisphere’s surface; hence some region of the hemisphere’s surface may not be seen by the eye.

Please develop a method to estimate the amount of the area of the hemisphere’s surface which can be seen by the eye? Here is a code for the rectangle rather than hemisphere:

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

Answer 1

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))

Answer 2

Take a small element of surface area centered on , dimensions . The area element is given by

The idea is to loop over these elements on the sphere; calculate the center point of the element at , and work out if the line segment between this point and the camera intersects a triangle. More here: https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm.

Now we need to find the points; trivially this means incrementing by over all of the hemisphere. But this would make the sampling resolution uneven - the factor would make the elements far larger near the apex of the hemisphere than at its edge.

Instead:

• Set a fixed number N of rings to loop around, i.e. number of iterations of .

• Set the minimum iteration area . The number of iterations of , M is gven by

  

  Where is the total area of the ring at :

  

• And of course from above the increments are given by

  

• Loop over all the rings, being careful that gives the middle line of each ring (So start at ); same concern need not apply to due to the symmetry. For each ring loop over each , and do the line intersection test as mentioned above.

The above method reduces the amount of bias in the area sampling resolution at small .

An even better way would be Fibonacci lattices, but they are more complicated. See this paper: http://geonaut.eu/published/016_Fibonacci_Lattice.pdf