Code covered by the BSD License  

Highlights from
circtriolap.m

image thumbnail
from circtriolap.m by Chris
Computes overlap area of generally placed triangle and circle.

circtriolap(rc,v1,v2,v3)
%%%  circtriolap : Calculates the overlap area for a generally located
%%%                plane triangle and a circle centered at 0,0
%%%		inputs: rc - radius of circle; v1, v2, v3 - vertex vectors [1 x 2] of plane triangle.
%%%		outputs: A - area of overlap with circle, vout = # vertices outside circle
%%%  Ex. [A,vout] = circtriolap(2,[1 2],[-3 2],[1 1]);  => A = 1.365, vout = 2 (two vertices ([1 2],[-3 2]) are outside)
%%%  The routine has been tested against CAD measurements for several
%%%  cases, but no guarantees. I have a feeling that there must be
%%%  a more elegant way to compute this, so if you have one, I'd love to see it.
%%%  (Have found it useful for computing irradiance from facted (or facet- sampled, surfaces]
%%%  CB, www.oesrochester.com
function [A, vout] = circtriolap(rc,v1,v2,v3)
switch nargin 
	case{2}
   	vo = v1;  % option for 3 vertex vectors input as 3 x 2 matrix
   case{4}      
      vo(1,:) = v1; vo(2,:) = v2; vo(3,:) = v3;
   otherwise
      A = NaN; vout = NaN;
      disp('Enter 3 x 2 matrix or (3) 1 x 2 vectors for vertices after circle radius argument.');
      return;
end
if (size(vo,1) ~= 3) ||  (size(vo,2) ~= 2)
      A = NaN; vout = NaN;
      disp('Enter 3 x 2 matrix or (3) 1 x 2 vectors for vertices after circle radius argument.');
      return;
end  
rc2 = rc*rc; Ac = pi*rc2;
%%%%%%%%%%%%%%%%%%%%%% sort vertices by distance from origin
[ds,vi] = sort([norm(vo(1,:)) norm(vo(2,:)) norm(vo(3,:))]);
v1p = vo(vi(1),:); v2p = vo(vi(2),:); v3p = vo(vi(3),:);
vminhat = v1p/ds(1);
%% Quick, non-exclusive filter for zero overlap area (if not satisfied, A can still be 0)
if (ds(1) >= rc) && (v2p*vminhat' >= rc) && (v3p*vminhat' >= rc) 
   A = 0; 
   vout = 3;
else   
%%%%%%%%%%%%%%%%
%%%%%%%		START COMPUTATION FOR OVERLAP AREA OF TRIANGLE AND CIRCLE
%%%%%%%%%%%%%%%%
%% Check for number of vertices outside circle
%%% vf is an outside-circle index vector, vf(n) = 1 if vertex is outside circle
vf = (ds > rc);
%%%
switch sum(vf) 
    case {0}
   % All vertices in circle (Overlap Area = A_triangle), A_tri = 1/2 cross prod. of any 2 leg vectors          
   	  vout = 0;
      A = 0.5*abs((v1p(1)-v3p(1))*(v2p(2)-v3p(2)) - (v1p(2)-v3p(2))*(v2p(1)-v3p(1))); 
   case {1}  % Nearest 2 vertices (v1p, v2p) in circle	
        vout = 1;
        vl = v3p - v1p; % triangle side vector
        vu = v3p - v2p; % triangle side vector    
      	bl = (sqrt((v3p*vl')^2 + (rc2 - v3p*v3p')*(vl*vl')) - vl*v3p')/(vl*vl'); 
      	al = v3p + bl*vl;
      	bu = (sqrt((v3p*vu')^2 + (rc2 - v3p*v3p')*(vu*vu')) - vu*v3p')/(vu*vu'); 
         au = v3p + bu*vu;
         alv = al - v3p; auv = au - v3p;
         At1 = 0.5*abs(alv(1)*auv(2)-alv(2)*auv(1));   					
         At2 = 0.5*abs(al(1)*au(2)-al(2)*au(1));
         At3 = 0.5*abs(vl(1)*vu(2)-vl(2)*vu(1));
      	 Aw = 0.5*acos((al*au')/rc2)*rc2;
         A =  Aw - At2 + At3 - At1;	% chord + difference in triangles          
   case {2}   % Only nearest-vertex (v1p) in circle
      vout = 2;
         vl = v2p - v1p; vu = v3p - v1p;     										
         bl = (sqrt((v1p*vl')^2 + (rc2 - (v1p*v1p'))*(vl*vl')) - vl*v1p')/(vl*vl');         	
         bu = (sqrt((v1p*vu')^2 + (rc2 - (v1p*v1p'))*(vu*vu')) - vu*v1p')/(vu*vu');
         al = v1p + bl*vl; au = v1p + bu*vu; % al, au - radial vectors to intersection of triangle w/ circle
         alv = al - v1p; auv = au - v1p; 
         Aw = 0.5*acos((al*au')/rc2)*rc2;
         At1 = 0.5*abs(alv(1)*auv(2)-alv(2)*auv(1));   					
         At2 = 0.5*abs(al(1)*au(2)-al(2)*au(1));
         A = Aw - At2 + At1;   % chord + triangle
% Check for chord inside circle from outside vertices
        b32 = (v3p-v2p)/norm(v3p-v2p); k3 = v3p*b32';
         if sqrt(v3p*v3p' - k3*k3) < rc  % perpendicular to chord inside circle < radius
             k2 = v2p*b32';
             gll = k3 - sign(k3)*sqrt(k3*k3 - (v3p*v3p' - rc2));
             guu = -k2 + sign(k2)*sqrt(k2*k2 - (v2p*v2p' - rc2));
             vll = v3p - gll*b32;
             vuu = v2p + guu*b32;   % vll, vuu are intercepts of chord on circle
             
             Aww = 0.5*acos((vll*vuu')/rc2)*rc2;
             Att = 0.5*abs(vll(1)*vuu(2) - vuu(1)*vll(2));
             A = A - Aww + Att;  % subtract chord area
         end           
  case {3}   % sum(vf) = 3; all vertices outside circle, most conditional case
     vout = 3; 
%        Form perp. vectors (qc) from origin to triangle sides 
%        Check if triangle sides make chords in the circle
         vd = (v3p - v2p); % side vector
         b(1) = (v3p*vd')/(vd*vd');  % rel. distance from furthest vertex along triangle side to perp (|b| > 1 => no chord)
         cperp(1,:) = v3p - b(1)*vd; % perp vector to triangle side
         nq(1) = norm(cperp(1,:));
%         
         vd = (v3p - v1p); % side vector;
         b(2) = (v3p*vd')/(vd*vd'); % (|b| > 1 => no chord)
         cperp(2,:) = v3p - b(2)*vd;
         nq(2) = norm(cperp(2,:));
%     
         vd = (v2p - v1p); % side vector
         b(3) = (v2p*vd')/(vd*vd'); % (|b| > 1 => no chord)
         cperp(3,:) = v2p - b(3)*vd;
         nq(3) = norm(cperp(3,:)); 
%                  
         [nqs,nqi] = sort(nq); % Sort acc. to shortest perp vector
         disp(nqs);
         nch = sum((nq < rc)); % number of POTENTIAL chords thru circle
         A = 0;
         switch nch
             case {0}
                 if (cperp(nqi(1),:)*(cperp(nqi(2),:))' <= 0)
                     A = Ac; % obtuse angle betw. perps => circle completely inside triangle
                 end                     
             case {1}
                  if abs(b(nqi(1))) > 1 % not a chord if this is true
                      A = 0;
                  else
                      if (cperp(nqi(1),:)*(cperp(nqi(2),:))' < 0)
                         A = Ac - chordarea(nqs(1),rc);
                      else
                          A = chordarea(nqs(1),rc);
                      end
                 end
             case {2}
                 A1 = 0; A2 = 0;
                 if abs(b(nqi(1))) <= 1
                    A1 = chordarea(nqs(1),rc);
                 end
                 if abs(b(nqi(2))) <= 1
                     A2 = chordarea(nqs(2),rc);
                 end
                 if (cperp(nqi(1),:)*(cperp(nqi(2),:))' < 0)
                     A = Ac - A1 - A2;
                 else
                     A = A1 - A2;
                 end
             case {3}
                 A1 = 0; A2 = 0; A3 = 0;
                 if abs(b(nqi(1))) <= 1
                    A1 = chordarea(nqs(1),rc);
                 end
                 if abs(b(nqi(2))) <= 1
                     A2 = chordarea(nqs(2),rc);
                 end
                 if abs(b(nqi(3))) <= 1
                     A3 = chordarea(nqs(3),rc);
                 end 
                 if (cperp(nqi(1),:)*(cperp(nqi(2),:))' < 0)
                     A = Ac - A1 - A2 - A3;
                     disp([A A1 A2 A3 -1]);
                 else
                     A = A1 - A2 - A3;
                     disp([A A1 A2 A3 1]);
                 end                                  
         end %% switch nch
end   %% switch sum(vf)
end  %% Simple check 
 %
function ch = chordarea(d,r)  % d is perp dist. to chord, r is circle radius
if d >= r
   ch = 0;
else
   u = d/r;
   ch = r*r*(acos(u) - u*sqrt(1 - u*u));
end

   

Contact us