已知三角形的端点和圆的圆心半径，求相交的面积

题目大概是已知三角形的三点坐标（ax,ay),(bx,by),(cx,cy)和圆的圆心（ox，oy）半径r，求三角形和圆的相交面积

这题很明显可以求得确定的面积，但是要考虑特别多种情况，没有一开始想的那么简单，代码是很早以前写的了，挺费脑子的，贴出来留个纪念。如果有大佬有更简单的算法的请告知
get_S

function res= get_S(ax,ay,bx,by,cx,cy,ox,oy,r)
    
         ans=0.0;
        %将圆心移动到原点 
        ax=ax-ox;
        ay=ay-oy;
        bx=bx-ox;
        by=by-oy;
        cx=cx-ox;
        cy=cy-oy;
        ox=0;
        oy=0;
        ans=get_area(ax,ay,bx,by,ox,oy,r)+get_area(bx,by,cx,cy,ox,oy,r)+get_area(cx,cy,ax,ay,ox,oy,r);
       
        res = abs(ans);
end

get_area

function res = get_area(ax,ay,bx,by,ox,oy,r)
sign=get_sign(ax,ay,bx,by);%判断有向面积的符号
ans=0;
oa=dis(ox,oy,ax,ay),ob=dis(ox,oy,bx,by),ab=dis(ax,ay,bx,by);
l=dis_line(ax,ay,bx,by);
if oa==0 || ob==0 %有一边边长为0则面积为0
    res= 0;
    return;
end
%第一种情况,两点都在圆内
if oa<=r&&ob<=r
    
    ans=abs(multi(ax,ay,bx,by))/2.0;%相交面积即为三角形oab的面积
    res =  sign*ans;
    return
    
    %第二种情况,两点都在圆外且圆心到直线ab距离大于半径
elseif (oa>=r&&ob>=r&&l>=r)
    
    [t1x,t1y]=get_point(ax,ay,r);%oa与圆的交点
    [t2x,t2y]=get_point(bx,by,r);%ob与圆的交点
    d=dis(t1x,t1y,t2x,t2y);
    ang=acos(get_angle(r,r,d));%角t1ot2
    ans=abs(ang*r*r/2.0);%相交面积即为扇形ot1t2的面积
    res= sign*ans;
    return
    
    %第三种情况,两点都在圆外且圆心到直线距离小于半径但三角形oab有一底角是钝角即直线ab与圆没有交点
    elsif (oa>=r&&ob>=r&&l<=r&&(get_angle(ab,oa,ob)<=0||get_angle(ab,ob,oa)<=0))
    
    [t1x,t1y]=get_point(ax,ay,r);%oa与圆的交点
    [t2x,t2y]=get_point(bx,by,r);%ob与圆的交点
    dist=dis(t1x,t1y,t2x,t2y);%线段t1t2长度
    ang=acos(get_angle(r,r,dist));%角t1ot2
    ans=abs(ang*r*r/2.0);%相交面积即为扇形ot1t2的面积
    res= sign*ans;
    return;
    
    %第四种情况,两点都在圆外且原点到直线ab距离小于半径且三角形oab底角均为锐角即直线ab与圆有两个交点
elseif(oa>=r&&ob>=r&&l<=r&&get_angle(ab,oa,ob)>0&&get_angle(ab,ob,oa)>0)
    
    %c,d为直线ab与圆的两个交点
    if(ax~=bx)%直线ab斜率存在
        
        k=(ay-by)/(ax-bx);%直线ab斜率
        h=ay-k*ax;%直线ab截距
        %解一元二次方程求圆与直线ab的交点
        a0=1.0+k*k;
        b0=2.0*k*h;
        c0=h*h-r*r;
        cx=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        cy=k*c.x+h;
        dx=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        dy=k*d.x+h;
        
    else%直线ab斜率不存在
        
        cx=ax;
        dx=ax;
        cy=sqrt(r*r-ax*ax);
        dy=-sqrt(r*r-ax*ax);
    end
    
    [t1x,t1y]=get_point(ax,ay,r);%oa与圆的交点
    [t2x,t2y]=get_point(bx,by,r);%ob与圆的交点
    d1=dis(cx,xy,dx,dy);%线段cd长度
    d2=dis(t1x,t1y,t2x,t2y);%线段t1t2长度
    ang1=acos(get_angle(r,r,d1));%角cod
    ang2=acos(get_angle(r,r,d2));%角t1ot2
    s1=abs(ang1*r*r/2.0);%小扇形ocd面积
    s2=abs(ang2*r*r/2.0);%小大扇形ot1t2面积
    s3=abs(multi(cx,xy,dx,dy))/2.0;%三角形ocd面积
    ans=s2+s3-s1;%相交面积即为三角形ocd面积加上两个扇形面积之差
    res= sign*ans;
    return
    
    %第五种情况,a点在圆外,b点在圆内
elseif(oa>=r&&ob<=r)
    
    %c,d两点为直线ab与圆的两个交点,e点为介于ab两点之间的那个交点
    if(ax~=bx)%直线ab斜率存在
        
        k=(ay-by)/(ax-bx);%直线ab斜率
        h=ay-k*ax;%直线ab截距
        %解一元二次方程求圆与直线ab的交点
        a0=1.0+k*k;
        b0=2.0*k*h;
        c0=h*h-r*r;
        cx=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        cy=k*cx+h;
        dx=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        dy=k*dx+h;
        %交点应为c,d两点中横坐标介于(a.x,b.x)两点之间的那一个
        if(ax<=cx&&cx<=bx||ax>=cx&&cx>=bx)
            ex=cx;
            ey=cy;
        else
            ex=dx;
            ey=dy;
        end
        
        
    else%直线ab斜率不存在
        
        dx=ax;
        cx=dx;
        cy=-sqrt(r*r-a.x*a.x);
        dy=sqrt(r*r-a.x*a.x);
        %交点应为c,d两点中纵坐标介于(a.y,b.y)两点之间的那一个
        if(a.y<=c.y&&c.y<=b.y||a.y>=c.y&&c.y>=b.y)
            ex=cx;
            ey=cy;
        else
            ex=dx;
            ey=dy;
        end
    end
    [ t1x,t1y]=get_point(ax,ay,r);%oa与圆的交点
    dist=dis(t1x,t1y,ex,ey);%线段t1e长度
    ang=acos(get_angle(r,r,dist));%角eot1
    s1=abs(ang*r*r/2.0);%扇形oet1面积
    s2=abs(multi(bx,by,ex,ey))/2.0;%三角形obe面积
    ans=s1+s2;%相交面积即三角形obe面积加上扇形oet1面积
    res= sign*ans;
    return;
    
    %第五种情况,b点在圆外,a点在圆内
elseif(ob>=r&&oa<=r)
    
    %c,d两点为直线ab与圆的两个交点,e点为介于ab两点之间的那个交点
    if(ax~=bx)%直线ab斜率存在
        
        k=(ay-by)/(ax-bx);%直线ab斜率
        h=ay-k*ax;%直线ab截距
        %解一元二次方程求圆与直线ab的交点
        a0=1.0+k*k;
        b0=2.0*k*h;
        c0=h*h-r*r;
        cx=(-b0+sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        cy=k*cx+h;
        dx=(-b0-sqrt(b0*b0-4.0*a0*c0))/(2.0*a0);
        dy=k*dx+h;
        %交点应为c,d两点中横坐标介于(a.x,b.x)两点之间的那一个
        if(ax<=cx&&cx<=bx||ax>=cx&&cx>=bx)
            ex=cx;
            ey=cy;
        else
            ex=dx;
            ey=dy;
        end
        
    else%直线ab斜率不存在
        
        dx=ax;
        cx=dx;
        cy=-sqrt(r*r-a.x*c.x);
        dy=sqrt(r*r-a.x*a.x);
        %交点应为c,d两点中纵坐标介于(a.y,b.y)两点之间的那一个
        if(ay<=cy&&cy<=by||ay>=cy&&cy>=by)
            ex=cx;
            ey=cy;
        else
            ex=dx;
            ey=dy;
        end
    end
    [t1x,t1y]=get_point(bx,by,r);%oa与圆的交点
    dist=dis(t1x,t1y,ex,ey);%线段t1e长度
    ang=acos(get_angle(r,r,dist));%角eot1
    s1=abs(ang*r*r/2.0);%扇形oet1面积
    s2=abs(multi(ax,ay,ex,ey))/2.0;%三角形oae面积
    ans=s1+s2;%相交面积即三角形oae面积加上扇形oet1面积
    res= sign*ans;
    return
end

以下是辅助函数：

function [x,y] = get_point(ax,ay,r)
if ax~=0 %斜率存在
    
    k=ay/ax;
    x=abs(r)/sqrt(1.0+k*k);
    if(ax<0)
        x=-x;%判断ans.x的符号
    end
    y=k*x;
    
    
else%斜率不存在
    
    x=0;
    if(ay>0) 
        y=r;%判断ans.y的符号
    else
        y=-r;
    end
    
end
end

function res = get_sign(ax,ay,bx,by)
if multi(ax,ay,bx,by)>0
    res= 1;
end
res= -1;
end

function res = dis(ax,ay,bx,by)
res=sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by));
end

function res = get_angle(a,b,c)
res=(a*a+b*b-c*c)/(2.0*a*b);
end

function res=dis_line(ax,ay,bx,by)
res=abs(multi(ax,ay,bx,by))/dis(ax,ay,bx,by);

function [theta1,theta2] = get_theta(x0,y0,r,x1,y1)
l2=(y1-y0)^2+(x1-x0)^2;

syms x y
[x,y]=solve((x-x0)^2+(y-y0)^2-r^2,(x-x1)^2+(y-y1)^2-l2); %得到交点的解析解
v1=[x(1)-x1,y(1)-y1];
v2=[x(2)-x1,y(2)-y1];
theta1=eval(acos(dot(v1,[1,0])/norm(v1)));
if eval(v1(2))<0
    theta1=2*pi-theta1;
end
theta2=eval(acos(dot(v2,[1,0])/norm(v2)));
if eval(v2(2))<0
    theta2=2*pi-theta2;
end
if theta1>theta2
    t=theta1;
    theta1=theta2;
    theta2=t;
end