% P=intp(X1,Y1,X2,Y2)
% 求两组离散点序列的交点
% X1,Y1为第一组的横纵坐标，X2,Y2为第二组横纵坐标，横坐标必须为有序数列。
% 返回值P为交点坐标，两列矩阵，第一列为横坐标，第二列为纵坐标。
function P=intp(X1,Y1,X2,Y2)
X1=X1(:);  % 变为列向量
X2=X2(:);
Y1=Y1(:);
Y2=Y2如何用matlab将已知点连线(:);
if max(X1)<min(X2) || max(X2)<min(X1)
    P=[];      % 两个区间没有重叠，不可能有交点
    return;
end
a=max(min(X1),min(X2));
b=min(max(X1),max(X2));
a1=find(X1>=a); a1=a1(1);
a2=find(X2>=a); a2=a2(1);
b1=find(X1<=b); b1=b1(end);
b2=find(X2<=b); b2=b2(end);
x=unique([X1(a1:b1); X2(a2:b2)]);
y1=interp1(X1,Y1,x,'linear');
y2=interp1(X2,Y2,x,'linear');  % 出公共部分
d=y1-y2;
ind0=find(d==0);  % 差为0的是交点。
P1=[x(ind0), y1(ind0)];
d1=sign(d);        % 求符号，负数为-1，正数为1, 0为0
d2=abs(diff(d1));
ind1=find(d2==2);  % 相邻符号相差2的，交点在此区间内
P2=zeros(length(ind1),2);
for k=1:length(ind1)
    i1=ind1(k);
    i2=ind1(k)+1;
    x1=x(i1);
    x2=x(i2);
    ya1=y1(i1);
    ya2=y1(i2);
    yb1=y2(i1);
    yb2=y2(i2);  % 两条线段四个端点坐标
    A=[ ya1-ya2, -(x1-x2)
        yb1-yb2, -(x1-x2)]; % 二元一次方程组系数矩阵
    B=[ (ya1-ya2)*x1-(x1-x2)*ya1
        (yb1-yb2)*x1-(x1-x2)*yb1];  % 常数项矩阵
    P2(k,:)=(A\B)';    % 解方程组，得到交点坐标
end
P=[P1;P2];            % 两种情况的交点合并
P=sortrows(P,1);      % 按横坐标排序
% 函数到此结束，保存到intp.m文件中
% 下面是接你的主程序来的,就是添在你画图那段程序后面
XA=AAAA;
YA=normpdf(AAAA,AA,AAA);
XB=BBBB;
YB=normpdf(BBBB,BB,BBB);
XC=CCCC;
YC=normpdf(CCCC,CC,CCC);
P1=intp(XA,YA,XB,YB);    % A,B交点，调用上面保存的那个intp函数
P2=intp(XB,YB,XC,YC);    % B,C交点
P3=intp(XC,YC,XA,YA);    % C,A交点
P=[P1;P2;P3];
% plot(XA,YA,XB,YB,XC,YC)
hold on
plot(P(:,1),P(:,2),'ro')
hold off
str=cell(size(P,1),1);
for k=1:size(P,1)
    str{k}=sprintf(' (%2.2g,%2.2g)',P(k,1),P(k,2));
end
text(P(:,1),P(:,2),str)
% 主程序到此结束
引言
    曾经思考过曲面求交，结果发现是学术界的一个难题，并且也想出了一个当前广泛使用方法原理一样的近似解法（追踪法）。当然网上也有很多方法，只不过那些方法非常粗糙，
无非就是meshgrid出离散网格，比较两曲面在某位置的坐标是否在某一精度范围内，然后标记显示之。这个方法仅仅当离散网格非常细的时候才比较精确。除此之外，还有个非常严重的问题：上面的“精度范围”不是你随心所欲给的，而且也没规律寻，当给得不恰当的时候，在格点处两曲面点作比较，会出很多个符合要求的点，或者一个也没有。这样就会使得交线非常曲折，甚至断裂等，严重影响精确度。
———————————————————分割线————————————————————————
    当然，既然有曲面求交，那么也有曲线求交，其基本结构就是两曲线求交。只是曲线求交问题，事先得澄清一些注意点：
    1. 数学分析层面求两曲线交点，其实就是方程组求解；
    2. “曲线”概念包括“直线”（处处曲率半径为无穷大）；
    3. Matlab的重点是离散点+矩阵运算，因此所有运算都是基于离散的，因而这里的曲线并不是绝对光滑的。
    4. 近似试探与未知函数表达式。
对于1，我想说的是，如果你想要求得两曲线的精确交点，并且一个不漏，那就直接求解方
程组，不用看本帖下文；
对于2，直线在Matlab里面是两个点确定，因此交点如果是一段线（无穷个点）的情况，可能只是显示两端点为交点；
对于3，很简单的例子，参数方程 x=cos(t),y=sin(t) (0<=t<=2*pi) 在数学分析（即连续空间）层面上是个圆，但是如果你在离散t的时候，间距比较大，那么最后Matlab绘制的图像不是圆，而是正多边形了。因此，此时我们讨论曲线交点是这个离散点连线的图形与其他图形的交点，而非圆与其他交点。这也是我在标题中加了“离散点连成”的修饰词，防止被误会。
对于4，既然是求曲线交点，那么本方法可以作为求方程组的近似解。当然，如果离散点够多，解的精确度可以保证，不过不能保证一个不漏。另外就是，对于一组离散点构成的曲线，很难知道它们的解析表达式，因此想通过非线性方程组求解的方法来求交点，就不大可能了（不过你可以用曲线拟合出函数解析式），因此，本帖的方法将会是一个较为有效求交点的方法。
    废话了那么多，下面就说说曲线求交点的方法吧。除了求解方程组，很多人想到的方法就是“离散点+判断距离是否足够接近”，这个方法原理跟引言中曲面求交的方法是一样的。
因此缺点也是一样的——太粗糙了。网上这种方法的代码也很多，这里就不上了。
下面将阐述我的方法以及给出例子代码。
    我有两种思路，一种是高级绘图层面的（不涉及到底层操作），一种是底层的。我只给出了第一种的代码，因为我不会底层操作。
    思路一：既然matlab曲线绘图是通过有序离散点依次连线形成，也就是说，通过“以直代曲”的过程，那么曲线交点无非就是离散点（结点）或者两线段交点。这比上面直接用交点附近的结点替代交点的方法要精确得多了。而两直线交点很容易求，只要知道四个点坐标，那么交点精确坐标自然可以表示出来。这就是求交点的原理。只是还有一些细节处理和要注意的地方，我会留到后面再详细说。
    思路二：仔细观察两曲线交点的特性，很容易发现，其实交点就是操作系统底层绘图重叠的那些像素点。因此，只要给要绘制的像素点做个标记，将那些重合的点突出显示（比如换个颜），那么就相当于显示出交点了。这种方法由于是本质性的，因此不会遗漏任何交点，而且精确度极高，适用范围广。Matlab提供的plot plot3 surf等绘图函数都属于高级绘图，底层绘图（或称低级绘图）只有line surface以及patch等少数函数。但是，这里的“底层”并非真正的底层，因为它还是经过封装了的，而C++的MFC里面直接用刷子绘图，
那才是依靠操作系统完成的真正的“底层”绘图操作（包括所有窗口都是操作系统绘制的）。这里扯远了，想要说明的就是底层绘图的概念而已。只是我不会用matlab实现这些底层绘图。
    上面说了思路，下面就详细说说一些注意点和需要处理的细节。
    为了算法的健壮性，就必须考虑各种奇异的情况，防止bug。我们要考虑曲线有分支（很多代数曲线是这样的，代数几何里面研究的东西）、间断跳跃（有绝对值函数或者存在渐近线情况）、首尾是交点、在切点相交，等等这些情况。而且对于定位交点处附近的四个最近端点也是个问题（因为这里存在一个情况，如果曲线1上的一条线段与曲线2上的两条或者以上的线段相交，我的程序因为这个问题没能有效解决，出现在一些非常特殊的情况下会遗漏部分交点）。上面的情况如果不考虑，那么你的程序就会出现各种各样的问题。
    对于通常情况，我考虑使用变号法则来判断交点（也就是高数里面“连续函数变号端点内存在零点”），对于上面说的特殊情况，那么预先处理，比如先看是否存在eps内的，或者为零的结点，有则直接记录，没有的话，通过两线段求交来确定交点。至于遍历顺序的问题，为了简便，我指考虑两曲线离散点个数相同的情况（因为不同的话，会出现一些无法
处理的情况），而且优先考虑离散点的坐标值中x或者y都相同的情况（比如x=0:0.1:pi; y1=sin(x), y2=x.^2这两条曲线的x值相同分布）。
下面是曲线y=cos(2*x).*exp(sin(x))与y2=sin(x).^2+cos(x)在[0:pi/18:2*pi]区间内的交点的代码：
注意：我没有写成接口的形式，虽然对于比那些较懒的人来说不太方便，但是这样做是为了让你能更好弄懂原理，并能自己改造代码。因此，下面的代码可以稍作修改，就能解决别的曲线求交点。这样，不愿思考的懒人就没法达到自己的目的了~
% 绘制两离散曲线的交点
% 注意：
%  1. 这里的“交点”指的是离散点连线绘出的图形的交点，而非函数或者方程理论分析上的交点，
%      因此，这个程序不能作为求根来用。
%  2. 要求两曲线的离散点的个数一样。
%  3. 两个曲线出现参数方程的话，大多数情况正常。但是经测试发现，对于某些非常特殊的情况会出现bug,
%      除非调用ezplot的数据（xdata,ydata）。
%
%  by kastin @Mar 21, 2012
clear;
debug=false; %关闭显示求交点过程
% 曲线1
x=0:pi/18:2*pi;
y=cos(2*x).*exp(sin(x));
% 曲线2
[x1 N]=sort(x);  %此处对于C1参数方程，C2为显式函数；或者均为参数方程时候有用
% 下面几句代码在本个案下没有什么特殊作用，但是当出现参数方程的时候，下面的方法改动一下就会有用。
y1=sin(x1).^2+cos(x1); %用于作图
x2=x;
y2=sin(x).^2+cos(x); %用于寻点
h=plot(x1,y1,'b',x,y,'c');
y(abs(y)<=eps)=0; y2(abs(y2)<=eps)=0;%对于三角函数关于零点的部分处理,但是发现sin(k*pi)不一定全在eps范围内
cy=y2-y; %作差
%符号记录
pos=cy>0;
neg=cy<=0;
%确定变号位置
fro=diff(pos)~=0; %变号的前导位置
rel=diff(neg)~=0; %变号的尾巴位置
zpf=find(fro==1); %记录索引
zpr=find(rel==1)+1; %记录索引
zpfr=[zpf; zpr];
hold on
% 观看求交点过程
if debug, hp=plot(x(zpfr),y(zpfr),'r.-',x2(zpfr),y2(zpfr),'g.-'); end
%线性求交
x0=(x(zpr).*(y2(zpf)-y(zpf))-x(zpf).*(y2(zpr)-y(zpr)))./(y(zpr)+y2(zpf)-y(zpf)-y2(zpr));
y0=y(zpf)+(x0-x(zpf)).*(y(zpr)-y(zpf))./(x(zpr)-x(zpf));
if any(isnan(y0)), y0=y2(zpf); end
%加入已经判断为零的位置
x0=[x(abs(cy)<=eps) x0].';
y0=[y(abs(cy)<=eps) y0].';
hc=plot(x0,y0,'k.'); %绘制交点
if debug, legend([h;hc;hp(1);hp(end)],'C1','C2','交点','微线段1','微线段2',0); end
legend([h;hc],'C1','C2','交点',0)
xlabel('x'), ylabel('y'), zlabel('z');
title('平面曲线交点')
axis equal
hold off
disp('交点坐标[x,y]为：')
disp(unique([x0,y0],'rows')) %排除重复的点
复制代码
经测试十几种奇怪的曲线相交（包括参数方程形式的曲线），目前发现上述代码的方法有四种情况会出现遗漏一两个交点。（其实上面代码本意是求显式函数的曲线交点，或者未知表达式的离散点曲线的交点，并未针对参数方程，隐函数方程做优化，但是可以凑合着用用。）

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。