/********************************************************************
创建日期:2017年3月17日
文 件 名:PDSNurbs.h
原始作者:孙宇飞
描    述:非均匀有理B样条

修改记录:
  版本号       修改日期    作者            修改内容
*********************************************************************/
#ifndef _PDS_NURBS_H_
#define _PDS_NURBS_H_

#include <math.h>
#include <afxwin.h>
#include <afxtempl.h>

const double PDS_NURBS_ZERO = 0.000001;  // 小于该数字的值认为等于0

//全局函数

//判断两个双精度数是否相等
BOOL PDSNurbs_IsEqual(double d1, double d2);

//判断两个双精度数是否不等
BOOL PDSNurbs_NotEqual(double d1, double d2);

//判断给定的双精度数是否等于0
BOOL PDSNurbs_EqualZero(double d);

//判断给定的双精度数是否不为0
BOOL PDSNurbs_NotEqualZero(double d);

//计算非均匀有理B样条使用的点
class CPDSNurbsPoint : public CObject
{
public:
  double x;
  double y;
  double z;
  double weight;  /* 权重,默认值:1
                  // 举例:(1,1,1,1)等同于(4,4,4,4),当为某一个点设置权重的时候,一定要将xyz分量全都除以weight才能保证位置不变 */

  double knot;  // 该点对应的节点数值,默认值:0

  int type;  // 类型,=1 转折点,=0 曲线点,默认值:0,[20180322 syf]增加该变量
             // 当this表示控制点时,启用该变量,那么在反求样条曲线的控制点链时,从转折点出断开,该变量不参与等号operator==的判断

public:
  CPDSNurbsPoint(void);
  CPDSNurbsPoint(CPDSNurbsPoint& rhs);
  CPDSNurbsPoint(CPoint ptPoint);
  CPDSNurbsPoint(double x, double y, double z = 0.0, double weight = 1.0);
  virtual ~CPDSNurbsPoint(void);
  CPDSNurbsPoint& operator=(CPDSNurbsPoint& rhs);
  CPDSNurbsPoint operator+(CPDSNurbsPoint& rhs);
  CPDSNurbsPoint operator-(CPDSNurbsPoint& rhs);
  void operator+=(CPDSNurbsPoint& rhs);
  BOOL operator==(CPDSNurbsPoint& rhs);

  void Initial(void);

  //将this的坐标转换为CPoint并返回
  CPoint GetPoint(void);

  //设置this的坐标
  void SetPoint(CPoint ptPoint);

  //乘以给定的数值
  void Multiple(double dParam);

  //除以给定的数值
  void Divide(double dParam);

  //将每个分量除以权重,得到规范后的坐标
  void DivideWeight(void);
};

//带类型的点的链表
class CPDSNurbsPointList : public CList<CPDSNurbsPoint,CPDSNurbsPoint>
{
public:  //扩展的变量
public:  //扩展的方法

  CPDSNurbsPointList& operator=(CPDSNurbsPointList& rhs);  //赋值
  CPDSNurbsPoint& operator[](int iIndex);  //根据索引值检索,不可越界

  //根据listPoint设置this
  //特别说明:
  //  (1)this中每个点z=0,weight=1
  void SetPoint(CList<CPoint,CPoint>& listPoint);

  //将this的数据输出至listPoint
  void GetPoint(CList<CPoint,CPoint>& listPoint);

  POSITION AddWithoutDuplicate(CPDSNurbsPoint ptNurbsPoint, BOOL bAddTail = TRUE);

  void AddWithoutDuplicate(CPDSNurbsPointList& listNurbsPoint, BOOL bAddTail = TRUE);

  //根据给定的节点数值,查找索引值
  int FindKnotIndex(double dKnot);

  //获取长度
  double GetLength(void);

  //将this翻转
  void Reverse(void);
};

//节点链表
//对于一条NURBS曲线来说(3次),该链表的长度应该等于控制点链长度+4
class CPDSNurbsKnotList : public CList<double,double>
{
public:  //扩展的变量
public:  //扩展的方法

  CPDSNurbsKnotList& operator=(CPDSNurbsKnotList& rhs);  //赋值
  double& operator[](int iIndex);  //根据索引值检索,不可越界

  //根据listKnot设置this
  void SetKnot(CList<double,double>& listKnot);

  //将this的数据输出至listKnot
  void GetKnot(CList<double,double>& listKnot);

  //添加新的节点,并且相邻的节点不能相同
  POSITION AddWithoutDuplicate(double dKnot, BOOL bAddTail = TRUE);

  //查找给定的节点在this中的重复次数
  int FindRepeatTime(double dKnot);
};

//矩阵
class CPDSMatrix : public CObject
{
protected:
  double** m_dMatrix;  // 二维数组的指针,为this私有,不同的对象之间不能共享m_dMatrix
  int m_iRowCount;  // 行数
  int m_iColumnCount;  // 列数

public:
  CPDSMatrix(void);
  CPDSMatrix(CPDSMatrix& rhs);
  virtual ~CPDSMatrix(void);
  CPDSMatrix& operator=(CPDSMatrix& rhs);
  CPDSMatrix operator*(CPDSMatrix& rhs);  // this的列数必须等于rhs的行数,否则不予计算

  //初始化,用于释放内存
  void Initial(void);

  //设置矩阵的大小
  //特别说明:
  //  在使用矩阵的数据之前,一定要首先设置大小
  void SetSize(int iRowCount, int iColumnCount);

  //获取矩阵的数据
  void GetMatrix(double**& dMatrix, int& iRowCount, int& iColumnCount);

  //获取二维数组的指针
  double** GetArray(void);

  //获取逆矩阵
  //特别说明:
  //  (1)this必须是方阵
  //  (2)本函数根据初等变换计算
  CPDSMatrix GetInverseMatrix(void);

  //输出调试信息
  void OutputDebugInfo(void);

  //将两行的数据加起来
  //输入参数:
  //  dArray1  数组1,它的数据被增加
  //  dArray2  数组2,其中的数据被加到dArray1,自身的数据不改变
  //  iCount   两个数组的长度
  //  dProp    比例值,dArray2[]*dProp被增加到dArray1
  void AddRow(double* dArray1, double* dArray2, int iCount, double dProp = 1.0);

  //将数组中的每一项除以指定的参数
  //输入参数:
  //  dArray  数组
  //  iCount  数组的长度
  //  dParam  参数,数组的每一项都要除以该值
  void DivideRow(double* dArray, int iCount, double dParam);
};

//NURBS曲线类
//特别说明:
//  (1)特指3次非均匀有理B样条曲线
//  (2)什么是型值点?
//     指的是经过的点,在我们的系统中称之为"控制点",实际上规范化的定义叫做型值点
//     型值点也叫做拟合点
//     英文名:data point
//  (3)什么是控制点?
//     用于计算实际点链的点,部分版师形象的称之为"手臂",曲线只会经过第一点与最后一点,不会经过中间的其他点
//     控制点链,也叫做特征多边形,它是非均匀有理B样条的凸包,曲线不会超过它的范围
//     英文名:control point
//  (4)NURBS曲线的方程由递推公式给出,详见相关文档
//  (5)什么是节点?
//     一组单调不减的数字,标识每一个control point的作用范围,对于相同的特征多边形来说,选择不同的节点,最终的效果不同
//  (6)节点有多少个?
//     对于degree=3来说,节点个数=特征多边形个数+degree+1=特征多边形个数+4
//     并且前degree+1=4个节点相等(未必等于0),末尾degree+1=4个节点相等(未必等于1)
//     例如:{A,A,A,A,B,C,D,E,E,E,E}就是7个control point所形成的节点,每一个{A,B,C,D,E}都对应于一个data point,而这个data point在屏幕上可以看到的点
//     如果{A,B,C,D,E}构成等差数列,曲线叫做均匀有理B样条,否则叫做非均匀有理B样条
//  (7)如果只有3个控制点,那么只能计算2次曲线
//  (8)如果只有2个控制点,那么只能得到一条线段
class CPDSNurbsCurve : public CObject
{
public:
  //将4条链表封装在一起,方便写代码,NURBS曲线的计算函数并不使用成员变量

  CPDSNurbsPointList m_listCtrlPoint;  // 控制点链,也叫特征多边形,最终的实际点链只经过头尾点,而不经过中间的其他点
  CPDSNurbsKnotList m_listKnot;  // 节点链表,长度等于m_listCtrlPoint.GetCount()+4
  CPDSNurbsPointList m_listDataPoint;  // 型值点链,经过的那些点
  CPDSNurbsPointList m_listRealPoint;  // 实际点链,插值后的点

public:
  CPDSNurbsCurve(void);
  CPDSNurbsCurve(CPDSNurbsCurve& rhs);
  virtual ~CPDSNurbsCurve(void);
  CPDSNurbsCurve& operator=(CPDSNurbsCurve& rhs);

  void Initial(void);

  //根据型值点计算开口线
  //输入参数:
  //  listDataPoint    型值点
  //  dMaxError        最大误差值
  //  bTangentH        =true 头点有切线,=false 头点没有切线
  //  ptTangentPointH  头点的切线方向,指的是实际点链GetPrev的方向
  //  bTangentT        =true 尾点有切线,=false 尾点没有切线
  //  ptTangentPointT  尾点的切线方向,指的是实际点链GetNext的方向
  //输出参数:
  //  listOutputPoint  实际点链
  //返回值:
  //  =true 计算成功,=false 计算失败
  BOOL Calculate_DataPoint_Open(CList<CPoint,CPoint>& listDataPoint, double dMaxError, CList<CPoint,CPoint>& listOutputPoint, BOOL bTangentH = FALSE, CPoint ptTangentPointH = CPoint(0,0), BOOL bTangentT = FALSE, CPoint ptTangentPointT = CPoint(0,0));

  //根据型值点和切线点计算闭合线
  //输入参数:
  //  listDataPoint    型值点,尾点可以等于头点,也可以不等,本函数自动处理
  //  dMaxError        最大误差值
  //输出参数:
  //  listOutputPoint  实际点链,输出时尾点与头点重合
  //返回值:
  //  =true 计算成功,=false 计算失败
  BOOL Calculate_DataPoint_Close(CList<CPoint,CPoint>& listDataPoint, double dMaxError, CList<CPoint,CPoint>& listOutputPoint);

  //根据实际点链反求型值点
  //输入参数:
  //  listOldPoint     旧的实际点链
  //  listDataPoint    型值点,如果输入时不为空,表示必须经过这些点
  //  dSplineError     计算样条曲线的最大误差值
  //  dMaxError        反求之后的曲线,与listOldPoint的最大误差值
  //  bTangentH        =true 头点有切线,=false 头点没有切线
  //  ptTangentPointH  头点的切线方向,指的是实际点链GetPrev的方向
  //  bTangentT        =true 尾点有切线,=false 尾点没有切线
  //  ptTangentPointT  尾点的切线方向,指的是实际点链GetNext的方向
  //输出参数:
  //  listDataPoint  型值点,如果输入时不为空,表示必须经过这些点
  //返回值:
  //  =true 计算成功,=false 计算失败
  BOOL ReCalculate_DataPoint(
    CList<CPoint,CPoint>& listOldPoint, CList<CPoint,CPoint>& listDataPoint, double dSplineError, double dMaxError,
    BOOL bTangentH = FALSE, CPoint ptTangentPointH = CPoint(0,0), BOOL bTangentT = FALSE, CPoint ptTangentPointT = CPoint(0,0));

  //根据给定的控制点链与节点链计算实际点链
  //特别说明:
  //  (1)本函数将会自动设置m_listCtrlPoint与m_listKnot,无需事先设置
  //  (2)输入时要保证listKnot.GetCount()=listCtrlPoint.GetCount()+4,这是NURBS曲线的要求,并且要保证节点数的正确性:单调不减,否则计算的结果就是错误的
  //  (3)从DXF等文件中读入的NURBS曲线,可以直接调用本函数来计算实际点链
  //输入参数:
  //  listCtrlPoint  控制点链
  //  listKnot       节点链表
  //  dMaxError      最大误差值
  //输出参数:
  //  listDataPoint  型值点链,即:曲线经过的那些点
  //  listRealPoint  实际点链
  //返回值:
  //  =true 计算成功,=false 计算失败
  BOOL Calculate_CtrlPoint(CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, double dMaxError, CPDSNurbsPointList& listDataPoint, CPDSNurbsPointList& listRealPoint);

  //根据给定的控制点链与节点链计算型值点
  //输入参数:
  //  listCtrlPoint  控制点链
  //  listKnot       节点链表
  //输出参数:
  //  listDataPoint  型值点链,即:曲线经过的那些点
  //返回值:
  //  =true 计算成功,=false 计算失败
  BOOL Calculate_CtrlPoint(CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, CPDSNurbsPointList& listDataPoint);

protected:

  //根据给定的型值点计算控制点(计算开口线)
  //输入参数:
  //  listDataPoint    型值点链,指的是曲线经过的那些点
  //  bTangentH        =true 头点有切线,=false 头点没有切线
  //  ptTangentPointH  头点的切线方向,指的是实际点链GetPrev的方向
  //  bTangentT        =true 尾点有切线,=false 尾点没有切线
  //  ptTangentPointT  尾点的切线方向,指的是实际点链GetNext的方向
  //输出参数:
  //  listCtrlPoint  控制点链,它是NURBS曲线的控制多边形或者说特征多边形,其长度=listDataPoint.GetCount()+2
  //  listKnot       节点链表
  void CalculateCtrlPoint_Open(CPDSNurbsPointList& listDataPoint, BOOL bTangentH, CPDSNurbsPoint ptTangentPointH, BOOL bTangentT, CPDSNurbsPoint ptTangentPointT, CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot);

  //根据给定的型值点计算控制点(闭合线)
  //特别说明:
  //  (1)listInputPoint的尾点可以等于头点,也可以不等于
  //  (2)最终输出的listCtrlPoint头尾点相等,表示再次回到头点上
  //输入参数:
  //  listInputPoint  型值点链,指的是曲线经过的那些点
  //输出参数:
  //  listCtrlPoint  控制点链,它是NURBS曲线的控制多边形或者说特征多边形
  //  listKnot       节点链表
  void CalculateCtrlPoint_Close(CPDSNurbsPointList& listInputPoint, CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot);

  //获取给定线段的长度
  //输入参数:
  //  listCtrlPoint  控制点链
  //  iIndex         段的索引值,指的是(listCtrlPoint[iIndex], [iIndex+1])所形成的线段,例如:第0段指的是([0],[1])
  //返回值:
  //  长度,如果传入的索引值越界则返回0
  double GetSegmentLength(CPDSNurbsPointList& listCtrlPoint, int iIndex);

  //查找给定节点的索引值
  //特别说明:
  //  (1)如果给定的节点位于[listKnot[i],listKnot[i+1]),那么返回i
  //  (2)注意可以等于listKnot[i],一定要严格小于[i+1]
  //  (3)从NURBS曲线的定义可以知道,由于{A,A,A,A...}前四个数据重复,所以该函数的返回值从3开始
  //输入参数:
  //  listKnot  节点链表,单调不减
  //  dKnot     节点
  //返回值:
  //  索引值,从0开始,若为-1则表示查找失败
  int FindKnotIndex(CPDSNurbsKnotList& listKnot, double dKnot);

  //根据给定的节点数值,计算B样条基函数的值
  //特别说明:
  //  (1)这是一个递归函数,需要调用iDegree-1的数据进行计算
  //  (2)在计算的过程中可能出现分母为0/0的情况,此处我们定义0/0=0
  //输入参数:
  //  listKnot  节点链表
  //  dKnot     给定的节点数值
  //  iIndex    dKnot所在的索引值
  //  iDegree   NURBS曲线的次数,一般等于3,如果只有三个控制点,那么只能生成2次曲线
  //返回值:
  //  B样条基函数的值
  double CalculateKnot(CPDSNurbsKnotList& listKnot, double dKnot, int iIndex, int iDegree);

  //根据给定的节点数值,计算NURBS曲线上的点坐标
  //输入参数:
  //  listCtrlPoint  控制点链
  //  listKnot       节点链表
  //  dKnot          给定的节点数值
  //  iDegree        样条曲线的次数
  //返回值:
  //  计算出的点位置
  CPDSNurbsPoint CalculateNurbsPoint(CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, double dKnot, int iDegree);

  //根据给定的两个节点,计算弦高
  //输入参数:
  //  listCtrlPoint     控制点链
  //  listKnot          节点链表
  //  ptPrevNurbsPoint  前一个点的坐标
  //  dPrevKnot         前一个点的节点
  //  ptNextNurbsPoint  后一个点的坐标
  //  dNextKnot         后一个点的节点
  //  iDegree           样条曲线的次数
  //返回值:
  //  弦高
  double CalculateChordHeight(
    CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot,
    CPDSNurbsPoint ptPrevNurbsPoint, double dPrevKnot, 
    CPDSNurbsPoint ptNextNurbsPoint, double dNextKnot, int iDegree);

  //根据给定的数量计算等分后的节点
  //特别说明:
  //  (1)若dPrevKnot=dNextKnot,那么输出链表为空
  //输入参数:
  //  dPrevKnot  前一个节点
  //  dNextKnot  后一个节点
  //  iCount     给定节点之间的点数量,若=2,表示有2个点,输出链表为{dPrevKnot,节点1,节点2,dNextKnot}
  //输入参数:
  //  listDivideKnot  分割后的节点链表,包括dPrevKnot与dNextKnot
  void CalculateDivideKnot(double dPrevKnot, double dNextKnot, int iCount, CPDSNurbsKnotList& listDivideKnot);

  //计算两个节点之间的实际点链
  //特别说明:
  //  (1)若dPrevKnot=dNextKnot,那么输出链表只有一个点
  //输入参数:
  //  listCtrlPoint  控制点链
  //  listKnot       节点链表
  //  dPrevKnot      前一个节点
  //  dNextKnot      后一个节点
  //  iDegree        样条曲线的次数
  //  dMaxError      最大误差值,即:每条线段的弦高不能超过该值
  //输出参数:
  //  listOutputPoint  实际点链,包含dPrevKnot与dNextKnot对应的点
  void CalculateNurbsPoint_BetweenTwoKnot(
    CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, double dPrevKnot, double dNextKnot,
    int iDegree, double dMaxError, CPDSNurbsPointList& listOutputPoint);

  //根据等分后的节点链表计算实际点链
  //输入参数:
  //  listCtrlPoint   控制点链
  //  listKnot        节点链表
  //  listDivideKnot  等分后的节点链表
  //  iDegree         样条曲线的次数
  //  dMaxError       最大误差值,即:每条线段的弦高不能超过该值
  //输出参数:
  //  listOutputPoint  实际点链,如果计算成功,则不为空,如果计算失败,该链表为空
  void CalculateNurbsPoint_UseDivideKnot(
    CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, CPDSNurbsKnotList& listDivideKnot,
    int iDegree, double dMaxError, CPDSNurbsPointList& listOutputPoint);

  //计算中点
  //输入参数:
  //  ptPrevNurbsPoint  前一个点,将其weight去掉之后再求解
  //  ptNextNurbsPoint  后一个点,将其weight去掉之后再求解
  //返回值:
  //  中点,它的weight=1
  CPDSNurbsPoint GetMiddleNurbsPoint(CPDSNurbsPoint ptPrevNurbsPoint, CPDSNurbsPoint ptNextNurbsPoint);

  //计算两个点之间的距离
  //输入参数:
  //  ptPrevNurbsPoint  前一个点,将其weight去掉之后再求解
  //  ptNextNurbsPoint  后一个点,将其weight去掉之后再求解
  //返回值:
  //  距离
  double TwoNurbsPointLength(CPDSNurbsPoint ptPrevNurbsPoint, CPDSNurbsPoint ptNextNurbsPoint);

  //根据给定点在线段上的参数,计算该点
  CPDSNurbsPoint ParamToNurbsPoint(CPDSNurbsPoint ptPrevNurbsPoint, CPDSNurbsPoint ptNextNurbsPoint, double dParam);

  //查找链表中的哪个点与给定点最近,返回其索引值
  //输入参数:
  //  listPoint   链表
  //  ptSelPoint  给定点
  //输出参数:
  //  dMinDistance  最近距离
  //返回值:
  //  索引值,从0开始
  int FindNearestIndex(CList<CPoint,CPoint>& listPoint, CPoint ptSelPoint, double& dMinDistance);

  //以累加弦长的方式生成节点
  //输入参数:
  //  listDataPoint  型值点,至少3个点
  //输出参数:
  //  listUniqueKnot  不重复的节点链表
  //  listKnot        重复的节点链表
  void GenerateKnot_ChordLength(CPDSNurbsPointList& listDataPoint, CPDSNurbsKnotList& listUniqueKnot, CPDSNurbsKnotList& listKnot);

  //计算垂足在矢量上的参数值
  //输入参数:
  //  ptNurbsPointS  矢量S
  //  ptNurbsPointE  矢量E
  //  ptNurbsPoint1  给定点
  //返回值:
  //  垂足点在矢量上的参数,如果SE点重合则直接返回0
  double GetNurbsPointVerticalParam(CPDSNurbsPoint ptNurbsPointS, CPDSNurbsPoint ptNurbsPointE, CPDSNurbsPoint ptNurbsPoint1);

  //获取给定点到链表的距离
  //输入参数:
  //  ptNurbsPoint1   给定点
  //  listNurbsPoint  链表
  //返回值:
  //  距离
  double DistanceOfNurbsPointPloyline(CPDSNurbsPoint ptNurbsPoint1, CPDSNurbsPointList& listNurbsPoint);

  //查找链表中的哪个点与给定点最近
  //输入参数:
  //  ptNurbsPoint1   给定点
  //  listNurbsPoint  链表
  //输出参数:
  //  dDistance  最近距离
  //返回值:
  //  最近点
  CPDSNurbsPoint FindNearestNurbsPoint(CPDSNurbsPoint ptNurbsPoint1, CPDSNurbsPointList& listNurbsPoint, double& dDistance);

  //判断新链表与旧链表的距离是否小于误差值
  //输入参数:
  //  listOldPoint  旧的链表
  //  listNewPoint  新的链表
  //  dMaxError     最大误差值
  //返回值:
  //  =true 新链表的每个点,到旧链表的距离都小于或者等于误差值,=false 某个点的距离大于误差值
  BOOL DistanceSmallerThanError(CList<CPoint,CPoint>& listOldPoint, CList<CPoint,CPoint>& listNewPoint, double dMaxError);

  //根据给定的控制点链与节点链计算实际点链,并且实际点链与旧点链的距离小于误差
  //输入参数:
  //  listCtrlPoint  控制点链
  //  listKnot       节点链表
  //  iSelIndex      给定的索引值,从0开始,从该索引值开始,向前和向后逐段计算,并比较偏差,该索引值附近的偏差是最大的
  //  dSplineError   计算NURBS曲线的最大误差值
  //  dMaxError      实际点链与旧点链的最大误差值
  //  listOldPoint   旧点链
  //返回值:
  //  =true 计算成功,每个点都小于等于误差值,=false 计算失败,某个点的距离大于误差值
  BOOL Calculate_CtrlPoint_SmallerThanError(CPDSNurbsPointList& listCtrlPoint, CPDSNurbsKnotList& listKnot, int iSelIndex, double dSplineError, double dMaxError, CList<CPoint,CPoint>& listOldPoint);
};

#endif  // _PDS_NURBS_H_