【IDL】幾何圖形數學運算函式

地理遥感生态网平台發表於2024-06-28
函式

幾何形狀，分為點、線、多邊形（面）、體幾類，利用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
分享：

相關文章