#include "stdafx.h"

//***************************************************************************************************
//                                        坡屋顶求体积
//***************************************************************************************************
AcGePoint3d GetIntegralPointX(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
    AcGePoint3d ptRet;
	//modified by fxj 比较浮点数大小用两个数做差与1E-6比
    if(p1.x - p2.x > 1E-6 && p1.x - p3.x < -1E-6 || p1.x - p3.x > 1E-6 && p1.x - p2.x < -1E-6)
        ptRet = p1;
    else if(p2.x - p1.x > 1E-6 && p2.x - p3.x < -1E-6 || p2.x - p3.x > 1E-6 && p2.x - p1.x < -1E-6)
        ptRet = p2;
    else if(p3.x - p1.x > 1E-6 && p3.x - p2.x < -1E-6 || p3.x - p2.x > 1E-6 && p3.x - p1.x < -1E-6)
        ptRet = p3;
    else{
        if(fabs(p1.x - p2.x) < 1E-6){
            if(p1.y > p2.y)
                ptRet = p1;
            else
                ptRet = p2;
        }
        else if(fabs(p1.x - p3.x) < 1E-6){
            if(p1.y > p3.y)
                ptRet = p1;
            else
                ptRet = p3;
        }
        else if(fabs(p2.x - p3.x) < 1E-6){
            if(p2.y > p3.y)
                ptRet = p2;
            else
                ptRet = p3;
        }
    }
    return ptRet;
}


double GetIntegralX(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
    double dRet(0);
    AcGePoint3d pt0, pt1, pt2;
    pt0 = GetIntegralPointX(p1, p2, p3);
    if(pt0 == p1){
        if(p2.x < p3.x){
            pt1 = p2;
            pt2 = p3;
        }
        else{
            pt1 = p3;
            pt2 = p2;
        }
    }
    else if(pt0 == p2){
        if(p1.x < p3.x){
            pt1 = p1;
            pt2 = p3;
        }
        else{
            pt1 = p3;
            pt2 = p1;
        }
    }
    else if(pt0 == p3){
        if(p1.x < p2.x){
            pt1 = p1;
            pt2 = p2;
        }
        else{
            pt1 = p2;
            pt2 = p1;
        }
    }

    double x0 = pt0.x;
    double x1 = pt1.x;
    double x2 = pt2.x;
    double y0 = pt0.y;
    double y1 = pt1.y;
    double y2 = pt2.y;
  
    double d1 = x0*x0-x1*x1;
    double d2 = x2*x2-x0*x0;
    double d3 = x0*x0*x0-x1*x1*x1;
    double d4 = x2*x2*x2-x0*x0*x0;

    
    double d6 = (y1-y2)/(x1-x2);
    
    double nTmp3 = 1.0/3.0;
	

	//modified by fxj 用三个点的顺逆时针判断倒正三角
	XPoly2D poly(1);
	poly.AppendNode(pt0);
	poly.AppendNode(pt1);
	poly.AppendNode(pt2);
	double polyArea=poly.PathArea();
	if(polyArea<-1E-6){//正三角形
        double nTmp = 0;
        if(fabs(x0-x1) > 1E-6 && fabs(x1 - x2) > 1E-6)//有2个点x坐标相同，只积分一次
		{
			double d5 = (y0-y1)/(x0-x1);
		    nTmp = 0.5*d1*(y0 - y1 - d5*x0 + d6*x1) + nTmp3*d3*(d5 - d6);
		}
		
        if(fabs(x1 - x2) > 1E-6 && fabs(x2 - x0) > 1E-6)
		{
			double d7 = (y2-y0)/(x2-x0);
            nTmp += 0.5*d2*(y2 - y1 - d7*x2 + d6*x1) + nTmp3*d4*(d7 - d6);
		}
        dRet = nTmp;
    }
    else if(polyArea>1E-6){                            //倒三角形
        double nTmp = 0;
        if(fabs(x0-x1) > 1E-6 && fabs(x1 - x2) > 1E-6)
		{
			double d5 = (y0-y1)/(x0-x1);
            nTmp = 0.5*d1*(y1 - y0 - d6*x1 + d5*x0) + nTmp3*d3*(d6 - d5);
		}
        if(fabs(x1 - x2) > 1E-6 && fabs(x2 - x0) > 1E-6)
		{
			double d7 = (y2-y0)/(x2-x0);
            nTmp += 0.5*d2*(y1 - y2 - d6*x1 + d7*x2) + nTmp3*d4*(d6 - d7);
		}
        dRet = nTmp;
    }
    return dRet;
}



double GetIntegralY( AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
    AcGePoint3d ptTmp1, ptTmp2, ptTmp3;
 
    ptTmp1.x = p1.y;
    ptTmp1.y = p1.x;
    ptTmp1.z = 0;

    ptTmp2.x = p2.y;
    ptTmp2.y = p2.x;
    ptTmp2.z = 0;
    ptTmp3.x = p3.y;
    ptTmp3.y = p3.x;
    ptTmp3.z = 0;
    return GetIntegralX(ptTmp1, ptTmp2, ptTmp3);

}


//函数：通过积分，得到空间三角形投影到XY平面的体积
//                入口参数：

//p1, p2, p3 代表 投影三角形的3个顶点坐标
//Z0：积分水平基准面标高
//返回：体积积分

//modified by fxj
double TUTIL_IMPEXP GetTriangleVolumeByIntegral(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3,double Z0=0)
{
	p1.z-=Z0;
	p2.z-=Z0;
	p3.z-=Z0;

    double dRet(0);

	AcGeVector3d vt1(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
	AcGeVector3d vt2(p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
	AcGeVector3d vtNormal = vt1.crossProduct(vt2);
	vtNormal.normalize();
	AcGePlane    plane(p1, vtNormal); //投影三角形的平面
	double dA(0), dB(0), dC(0), dD(0); //dA, dB, dC,dD代表 aX + bY + cZ = D 空间三角形的平面方程
	plane.getCoefficients(dA, dB, dC, dD);
	dD = dD*(-1);
	AcGeVector3d vtNormal2 = plane.normal();
	
	BOOL bDzero=FALSE;
	if(fabs(dD)>1E-6)//平面方程转化为 aX + bY + cZ = 1 形式
	{
		dA = dA / dD ;
		dB = dB / dD ;
		dC = dC / dD ;
		
		bDzero=FALSE;
	}
	else
		bDzero=TRUE;
				
	if(fabs(dC)<1E-6)
		return dRet;
    
    XPoly2D poly;
    poly.AppendNode(p1);
    poly.AppendNode(p2);
    poly.AppendNode(p3);
    poly.AppendNode(p1);
	poly.SetElevation(0);

    double nProjectArea=poly.PathArea();
    double nCalcTmp = 1/dC*nProjectArea;
    double nX = dA / dC * GetIntegralX(p1, p2, p3);
    double nY = dB / dC *GetIntegralY(p1, p2, p3);

	if(nProjectArea<-1E-6)
	{
		nX=-nX;
		nY=-nY;
	}//积分是有方向的,与数学上积分的概念有所不同

	
    dRet= - nX - nY;
    if(bDzero!=TRUE)
		dRet += nCalcTmp;
    return dRet;
}
//对实心实体的体积计算 add by fxj 
double VolumeEntity(AcDbObjectId id)
{
	TMasterProperties prop;
    AcDbVoidPtrArray entSet;
	double dRet=0;
	AcDbEntity *pEnt=NULL;
	acdbOpenObject(pEnt,id,AcDb::kForRead);
	if(pEnt!=NULL) 
	{
		int nOldMode=DocGetDrawMode();
		DocSetDrawMode(DRAW_3D);
	   	
        AcDbVoidPtrArray TriangleSet;
		InitModel(prop, pEnt, NULL, entSet);
        DocSetDrawMode(nOldMode);
        AcGeVector3d vtNormalZ(0,0,1);

        int nLen=entSet.length();
        for(int n=0; n<nLen; n++) {
            AcDbEntity *pEnt=(AcDbEntity *)entSet[n];
            if(pEnt->isKindOf(AcDbFace::desc())) {
                AcDbFace *pFace = AcDbFace::cast(pEnt);
                AcGePoint3d p1,p2,p3,p4;
                pFace->getVertexAt(0, p1);
                pFace->getVertexAt(1, p2);
                pFace->getVertexAt(2, p3);
                pFace->getVertexAt(3, p4);
                AcGeVector3d vt1(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
                AcGeVector3d vt2(p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
                AcGeVector3d vtNormal = vt1.crossProduct(vt2);
                vtNormal.normalize();
                if(p1 == p4)
				{
                    TriangleSet.append(pEnt->clone());
                }
                else
				{
                    //注意保证逆时针方向
                    AcDbFace *pFaceNew1 = new AcDbFace();
                    pFaceNew1->setVertexAt(0, p1);
                    pFaceNew1->setVertexAt(1, p2);
                    pFaceNew1->setVertexAt(2, p3);
                    pFaceNew1->setVertexAt(3, p1);
                    TriangleSet.append(pFaceNew1);
                    AcDbFace *pFaceNew2 = new AcDbFace();
                    pFaceNew2->setVertexAt(0, p1);
                    pFaceNew2->setVertexAt(1, p3);
                    pFaceNew2->setVertexAt(2, p4);
                    pFaceNew2->setVertexAt(3, p1);
                    TriangleSet.append(pFaceNew2);
                }
            }
        }
        ReleaseEntSet(entSet);
        nLen=TriangleSet.length();
        
        for(n=0; n<nLen; n++) {
            AcDbEntity *pEnt=(AcDbEntity *)TriangleSet[n];
            if(pEnt->isKindOf(AcDbFace::desc())) {
                AcDbFace *pFace = AcDbFace::cast(pEnt);
                AcGePoint3d p1,p2,p3;
                pFace->getVertexAt(0, p1);
                pFace->getVertexAt(1, p2);
                pFace->getVertexAt(2, p3);
                double nV=0;
				
				nV=GetTriangleVolumeByIntegral(p1, p2, p3);//三角形积分计算
			
                dRet +=nV;
			}
        }
        ReleaseEntSet(TriangleSet);
     
		pEnt->close();
    }
	return -dRet;//与法向量有关
}
void VolumeRoofCmd()
{
    Entity ent;
    RBList rbMask(RTDXF0, _T("3DSOLID"),  NULL);
	if(SelectEntity(_T("\n请选择三维实体<退出>:"), ent, rbMask) != RTNORM)
        return;
    AcDbObjectId id = NULL;
    acdbGetObjectId(id, ent);
    double dV = VolumeEntity(id);
    adsout << _T("\n三维实体的体积为：") << dV*1e-9 << _T("立方米\n");
}