【IDL】几何图形数学运算函数

几何形状，分为点、线、多边形（面）、体几类，利用IDL对这些形状的几何运算，大致分下面几个部分。

1、 点集运算

点与点之间求距离：DISTANCE_MEASURE（IDL自带）

2、 线相关

计算点到直线的距离PNT_LINE（IDL自带）或CalDistancePtoLine.pro

计算两直线的交点CAL2LINESINTERSECTPOINT.pro

线段与坐标轴的夹角（数学坐标系）cal2pointsangle.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自带）

求两个平面的夹角cal2planeangle.pro

四面体与平面相交：TETRA_CLIP（IDL自带）

4、 体相关

体数据任意方向切面：EXTRACT_SLICE（IDL自带）

附源码：

;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;计算两个点的距离
;Write By DYQ
;可用distance_measure函数替代
;
function CalDistance, point1, point2
    compile_opt idl2
    ;
    Return, SQRT((point1[0]-point2[0])^2+(point1[1]-point2[1])^2+(point1[2]-point2[2])^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,linePos2

    a = CalDistance(point0,linePos1)
    b = CalDistance(point0,linePos2);SQRT((point[0]-linePos2[0])^2+(point[1]-linePos2[1])^2)

    c = CalDistance(linePos1,linePos2);SQRT((linePos2[0]-linePos1[0])^2+(linePos2[1]-linePos1[1])^2)

    p = (a+b+c)*0.5
    s = SQRT(p*(p-a)*(p-b)*(p-c))

    Return, s*2/c
end

;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;求集合的并集
; copy From DAVID's Code  ^_^
;
FUNCTION SETUNION, a, b
  ;
  COMPILE_OPT StrictArr
  IF N_ELEMENTS(a) EQ 0 THEN RETURN, b    ;A union NULL = a
  IF N_ELEMENTS(b) EQ 0 THEN RETURN, a    ;B union NULL = b
  RETURN, WHERE(HISTOGRAM([a,b], OMin = omin)) + omin ; Return combined set
END

;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;
;去除2维数组中重复的元素
;
; 2010年3月25日
; Write By DYQ
;
FUNCTION UNIQARRAY, inArray
  ;
  num = N_ELEMENTS(inArray[0,*])
  xArray = REFORM(inArray[0,*])
  yArray = REFORM(inArray[1,*])
  ;
  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((point1[0]-point2[0])^2+(point1[1]-point2[1])^2)
END

;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;计算点到直线的距离
;
;Write By DYQ
;其实可以用 PNT_LINE 函数，没办法，写好了才发现IDL自带o(∩_∩)o
;
FUNCTION CALDISTANCEPTOLINE, point0,linePos1,linePos2

  a = CALDISTANCE(point0,linePos1)
  b = CALDISTANCE(point0,linePos2);SQRT((point[0]-linePos2[0])^2+(point[1]-linePos2[1])^2)

  c = CALDISTANCE(linePos1,linePos2);SQRT((linePos2[0]-linePos1[0])^2+(linePos2[1]-linePos1[1])^2)

  p = (a+b+c)*0.5
  s = SQRT(p*(p-a)*(p-b)*(p-c))

  RETURN, s*2/c
END
;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;计算两直线的交点
;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 BEGIN
    IF ARRAY_EQUAL(P1S, P2S ) THEN BEGIN
     IF P1 EQ P2 THEN RETURN,P1
     RETURN, -1
    ENDIF ELSE BEGIN
     IF CALDISTANCEPTOLINE(P1,P1S,P2S) EQ 0 THEN RETURN, P1
     RETURN,-1
   ENDELSE
  ;线1的点不重合
  ENDIF ELSE BEGIN
    IF ARRAY_EQUAL(P1S ,P2S) THEN BEGIN
     RETURN, P1S
    ENDIF ELSE BEGIN
     ;如果第一条直线垂直x轴
     IF (P1[0]-P2[0]) EQ 0 THEN BEGIN
       ipX = p1[0];
       ;第二条也垂直
       IF (P1S[0]-P2S[0]) EQ 0 THEN BEGIN
         ;不相交
         RETURN, -1
       ENDIF ELSE BEGIN
         ;
         k2 = FLOAT(P1S[1]-P2S[1])/(P1S[0]-P2S[0])
         b2 = FLOAT(P1S[0]*P2S[1]-P2S[0]*P1S[1])/(P1S[0]-P2S[0])
         ;
         ipY = k2*ipX+b2
         RETURN,[ipX,ipY]
       ENDELSE
     ;第二条直线垂直X轴
     ENDIF ELSE IF (P1S[0]-P2S[0]) EQ 0 THEN BEGIN
       ipX = p2s[0];
       k1 = FLOAT(P1[1]-P2[1])/(P1[0]-P2[0])
       b1 = FLOAT(P1[0]*P2[1]-P2[0]*P1[1])/(P1[0]-P2[0])
       ipY = k1*ipX+b1
       RETURN,[ipX,ipY]

     ;都不垂直
     ENDIF ELSE BEGIN
       k1 = FLOAT(P1[1]-P2[1])/(P1[0]-P2[0])
       b1 = FLOAT(P1[0]*P2[1]-P2[0]*P1[1])/(P1[0]-P2[0])
       ;
       k2 = FLOAT(P1S[1]-P2S[1])/(P1S[0]-P2S[0])
       b2 = FLOAT(P1S[0]*P2S[1]-P2S[0]*P1S[1])/(P1S[0]-P2S[0])
       ;如果都垂直Y轴
       IF (K2 EQ 0) AND(K1 EQ 0) THEN RETURN,-1
       ipX = (b2-b1)/(k1-k2)
       ipY = (k1*b2-k2*b1)/(k1-k2)
       RETURN,[ipX,ipY]
     ENDELSE
   ENDELSE
  ENDELSE
;
END

;≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌≌
;
;
;计算平面上直线是否与多边形相交，有交点可返回
;2010年3月25日
; Write By DYQ
;
;调用格式result = CalLineInPolygon(linePoints, polyPoints) ;
;LinePoints: 两个点的二维坐标[2,2]，否则返回-1;
;PolyPoints: 多边形点的二维坐标[2,n]，否则返回-1
;
FUNCTION CALLINEINPOLYGON, linePoints, polyPoints
  COMPILE_OPT idl2
  ;
  IF ARRAY_EQUAL(SIZE(inPoints,/dimension), [2,2]) THEN RETURN,-1
  ;点个数
  pointNum = (SIZE(polyPoints,/dimension))[1]
  ;相交标识
  sign = -1
  ipoints = [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 BEGIN
     ipoints = [[iPoints],[REFORM(point)]]
     sign = 1
    ENDIF
  ENDFOR
  ;
  IF sign EQ -1 THEN RETURN, sign $
  ELSE result = UNIQARRAY(iPoints[*,1:N_ELEMENTS(iPoints[0,*])-1])
  ;判断交点是否在多边形内
  ;
  oROI = OBJ_NEW('IDLanROI', polyPoints)

  conArr = [0,0]
  FOR i=0, N_ELEMENTS(result)/2-1 DO BEGIN
    ;
    psign= oROI->CONTAINSPOINTS(result[*,i])
   ;有在多边形内的点
    IF psign NE 0 THEN BEGIN
     conArr = [[conArr], [result[*,i]]]
    ENDIF
  ENDFOR
  OBJ_DESTROY,oROI
  ;
  IF N_ELEMENTS(conArr)/2 EQ 1 THEN RETURN, -1 ELSE RETURN,conArr[*,1:N_ELEMENTS(conArr)/2-1]
END

;
;;
;根据经过两点线段与坐标轴的夹角(弧度)
;Point1到point2线段
; Write By DYQ
; 2007-10-27
;
FUNCTION Cal2pointsangle ,point1,point2
  ;
  point1 = FLOAT(point1)
  point2 = FLOAT(point2)
  ;
  dx = point2[0]- point1[0]
  dy = point2[1]- point1[1]

  IF dx EQ 0 THEN BEGIN ; x轴
    IF dy GT 0 THEN result = !PI/2 ELSE $
     result = 3*!PI/2

  ENDIF ELSE BEGIN
    IF dy EQ 0 THEN BEGIN ;x轴
     IF dx GT 0 THEN result = 0 ELSE $
       result = !PI
    ENDIF ELSE BEGIN
     result = ATAN(dy/dx)
     IF result GT 0 THEN IF dx LT 0 THEN result = result +!Pi
     IF result LT 0 THEN IF dx LT 0 THEN result = result +!Pi ELSE $
       result = result + 2*!PI
   ENDELSE
  END

  RETURN,result MOD (2*!PI)
END

;
;求两个平面的夹角
;平面方程：
; a1*x+b1*y+c1*z+d1 = 0
; a2*x+b2*y+c2*z+d2 = 0
; 返回角度（锐角）
function cal2PlaneAngle, plane1,plane2
  a1 = plane1[0]
  b1 = plane1[1]
  c1 = plane1[2]
  d1 = plane1[3]
  ;
  a2 = plane2[0]
  b2 = plane2[1]
  c2 = plane2[2]
  d2 = plane2[3]
 
  ;
  n1 = SQRT(a1^2+b1^2+c1^2)
  n2 = SQRT(a2^2+b2^2+c2^2)
  angleCos = ABS(a1*a2+b1*b2+c1*c2)/(n1*n2)
  angle = aCos(angleCos)/!dtor
  return ,angle
end
分享：