IDL中几何图形的数学运算函数

IDL中几何图形的数学运算函数几何形状，分为点、线、多边形（面）、体几类，利用IDL对这些形状的几何运算，大致分下面几个部分。1、 点集运算点与点之间求距离：DISTANCE_MEASURE（IDL自带）2、 线相关计算点到直线的距离PNT_LINE（IDL自带）或CalDistancePtoL计算两直线的交点CAL2LINESINTERSECTPOINT.pro线段与坐标轴的夹角（数学坐标系）3、 多边形相关（面）求多边形面积：poly_area、(IDLanROI)- ComputeGeometry, area = result（IDL自带）求多边形周长：(IDLanROI)- ComputeGeometry, PERIMETER = result（IDL自带）点坐标是否在多边形范围内：(IDLanROI)-ContainsPoints（IDL自带）三点求通过该系列三个点的圆心坐标和圆半径；CIR_3PNT（IDL自带）多边形与曲面交集：MESH_CLIP（IDL自带）两个多边形合并：MESH_MERGE（IDL自带）多边形是否空间闭合：MESH_ISSOLID（IDL自带）多边形包含的三角形个数：MESH_NUMTRIANGLES（IDL自带）复杂多边形正确显示：IDLgrTessellator（IDL自带）求两个平面的夹角四面体与平面相交：TETRA_CLIP（IDL自带）4、 体相关体数据任意方向切面：EXTRACT_SLICE（IDL自带）重点推荐MESH_OBJ，可以创建三角形、矩形、圆锥、圆柱、任意多边形旋转等。;Mesh_Obj,2, Vertex_List, Polygon_List, $; 调用格式 ; Mesh_Obj, Type, Vertex_List, Polygon_List , $; Array1 , Array2 , /CLOSED , /DEGREES , $; P1, P2, P3, P4, P5 = value;Type = 6 创建旋转后模型的点坐标及连接关系;Vertex_List 创建后的点坐标集;Polygon_List 旋转模型的点连接关系;Array1 创建旋转模型的起始数据集;Array2 Type 为6时不支持;/Closed 旋转后模型闭合;/DEGREES 旋转的输入参数为角度，默认是弧度;P1 旋转后的边数;P2 - 旋转的中心点;P3 - 旋转的方向0,0,1;P4 - 旋转起始角度;P5 - 旋转结束角度附提到的系列源码： CODE:;计算两个点的距离;Write By DYQ;可用distance_measure函数替代;function CalDistance, point1, point2compile_opt idl2;Return, SQRT(point10-point20)2+(point11-point21)2+(point12-point22)2)end; Write By ZML;计算点到直线的距离;s = SQRT(p*(p-a)*(p-b)*(p-c)p = (a+b+c)/2;Modified By DYQ;function CalDistancePtoLine, point0,linePos1,linePos2a = CalDistance(point0,linePos1)b = CalDistance(point0,linePos2);SQRT(point0-linePos20)2+(point1-linePos21)2)c = CalDistance(linePos1,linePos2);SQRT(linePos20-linePos10)2+(linePos21-linePos11)2)p = (a+b+c)*0.5s = SQRT(p*(p-a)*(p-b)*(p-c)Return, s*2/cend CODE:;求集合的并集; copy From DAVIDs Code _;FUNCTION SETUNION, a, b;COMPILE_OPT StrictArrIF N_ELEMENTS(a) EQ 0 THEN RETURN, b ;A union NULL = aIF N_ELEMENTS(b) EQ 0 THEN RETURN, a ;B union NULL = bRETURN, WHERE(HISTOGRAM(a,b, OMin = omin) + omin ; Return combined setEND;去除2维数组中重复的元素; 2010年3月25日; Write By DYQ;FUNCTION UNIQARRAY, inArray;num = N_ELEMENTS(inArray0,*)xArray = REFORM(inArray0,*)yArray = REFORM(inArray1,*);xUniq = UNIQ(xArray(SORT(xArray)yUniq = UNIQ(yArray(SORT(yArray);RETURN, inArray*,SETUNION(xUniq,yUniq)END;计算两个点的距离,XY水平面; 2008-4-28; Write By DYQ;其实可以用 DISTANCE_MEASURE 函数，没办法，写好了才发现IDL自带o(_)o;FUNCTION CALDISTANCE, point1, point2;point1 = point1*1.point2 = point2*1.RETURN, SQRT(point10-point20)2+(point11-point21)2)END;计算点到直线的距离;Write By DYQ;其实可以用 PNT_LINE 函数，没办法，写好了才发现IDL自带o(_)o;FUNCTION CALDISTANCEPTOLINE, point0,linePos1,linePos2a = CALDISTANCE(point0,linePos1)b = CALDISTANCE(point0,linePos2);SQRT(point0-linePos20)2+(point1-linePos21)2)c = CALDISTANCE(linePos1,linePos2);SQRT(linePos20-linePos10)2+(linePos21-linePos11)2)p = (a+b+c)*0.5s = SQRT(p*(p-a)*(p-b)*(p-c)RETURN, s*2/cEND;计算两直线的交点;2008-9-27;Write By DYQ; 输入; P1:直线1的起点坐标- P2:直线1的终点坐标; P1S:直线2的起点坐标- P2S:直线2的终y点坐标; 相交则返回交点，不相交则返回-1; Modified By DYQ; 2008-10-29 -增加了p1、p2及P1S、P2S重合的判断;Example:; 直线: -1,-1,1,1; 直线：1,0,0,1; point = Cal2linesintersectpoint(-1,-1,1,1,1,0,0,1)FUNCTION CAL2LINESINTERSECTPOINT, P1,P2,P1S,P2S;COMPILE_OPT IDL2;如果线1点重合了IF ARRAY_EQUAL(P1, P2) THEN BEGINIF ARRAY_EQUAL(P1S, P2S ) THEN BEGINIF P1 EQ P2 THEN RETURN,P1RETURN, -1ENDIF ELSE BEGINIF CALDISTANCEPTOLINE(P1,P1S,P2S) EQ 0 THEN RETURN, P1RETURN,-1ENDELSE;线1的点不重合ENDIF ELSE BEGINIF ARRAY_EQUAL(P1S ,P2S) THEN BEGINRETURN, P1SENDIF ELSE BEGIN;如果第一条直线垂直x轴IF (P10-P20) EQ 0 THEN BEGINipX = p10;第二条也垂直IF (P1S0-P2S0) EQ 0 THEN BEGIN;不相交RETURN, -1ENDIF ELSE BEGIN;k2 = FLOAT(P1S1-P2S1)/(P1S0-P2S0)b2 = FLOAT(P1S0*P2S1-P2S0*P1S1)/(P1S0-P2S0);ipY = k2*ipX+b2RETURN,ipX,ipYENDELSE;第二条直线垂直X轴ENDIF ELSE IF (P1S0-P2S0) EQ 0 THEN BEGINipX = p2s0;k1 = FLOAT(P11-P21)/(P10-P20)b1 = FLOAT(P10*P21-P20*P11)/(P10-P20)ipY = k1*ipX+b1RETURN,ipX,ipY;都不垂直ENDIF ELSE BEGINk1 = FLOAT(P11-P21)/(P10-P20)b1 = FLOAT(P10*P21-P20*P11)/(P10-P20);k2 = FLOAT(P1S1-P2S1)/(P1S0-P2S0)b2 = FLOAT(P1S0*P2S1-P2S0*P1S1)/(P1S0-P2S0);如果都垂直Y轴IF (K2 EQ 0) AND(K1 EQ 0) THEN RETURN,-1ipX = (b2-b1)/(k1-k2)ipY = (k1*b2-k2*b1)/(k1-k2)RETURN,ipX,ipYENDELSEENDELSEENDELSE;END;计算平面上直线是否与多边形相交，有交点可返回;2010年3月25日; Write By DYQ;调用格式result = CalLineInPolygon(linePoints, polyPoints) ;LinePoints: 两个点的二维坐标2,2，否则返回-1;PolyPoints: 多边形点的二维坐标2,n，否则返回-1;FUNCTION CALLINEINPOLYGON, linePoints, polyPointsCOMPILE_OPT idl2;IF ARRAY_EQUAL(SIZE(inPoints,/dimension), 2,2) THEN RETURN,-1;点个数pointNum = (SIZE(polyPoints,/dimension)1;相交标识sign = -1ipoints = 0,0;循环求解FOR i=0, pointNum -1 DO BEGIN;point = CAL2LINESINTERSECTPOINT(linePoints*,0, $linePoints*,1, polyPoints*,i,polyPoints*,(i+1) MOD pointNum)IF N_ELEMENTS(point) GT 1 THEN BEGINipoints = iPoints,REFORM(point)sign = 1ENDIFENDFOR;IF sign EQ -1 THEN RETURN, sign $ELSE result = UNIQARRAY(iPoints*,1:N_ELEMENTS(iPoints0,*)-1);判断交点是否在多边形内;oROI = OBJ_NEW(IDLanROI, polyPoints)conArr = 0,0FOR i=0, N_ELEMENTS(result)/2-1 DO BEGIN;psign= oROI-CONTAINSPOINTS(result*,i);有在多边形内的点IF psign NE 0 THEN BEGINconArr = conArr, result*,iENDIFENDFOROBJ_DESTROY,oROI;IF N_ELEMENTS(conArr)/2 EQ 1 THEN RETURN, -1 ELSE RETURN,conArr*,1:N_ELEMENTS(conArr)/2-1END;测试程序PRO TEST; restore,C:UsersAdministratorDesktopCalLineInPolygon.savlinePoints = 0,0,2,1polyPoints = 1,1,2,1,3,2,4,1points = CALLINEINPOLYGON(linePoints, polyPoints)PRINT,pointsEND; CODE:;根据经过两点线段与坐标轴的夹角(弧度);Point1到point2线段; Write By DYQ; 2007-10-27;FUNCTION Cal2pointsangle ,point1,point2;point1 = FLOAT(point1)point2 = FLOAT(point2);dx = point20- point10dy = point21- point11IF dx EQ 0 THEN BEGIN ; x轴IF dy GT 0 THEN result = !PI/2 ELSE $result = 3*!PI/2ENDIF ELSE BEGINIF dy EQ 0