ピックの定理

 投稿者：山中和義  投稿日：2010年12月 7日(火)15時59分40秒
  !ピックの定理（Pick's theorem）で、格子多角形の面積を求める

LET N=8 !頂点の数

DATA 4,0 !頂点の座標　※反時計まわり
DATA 2,1
DATA 1,3
DATA 0,1
DATA -3,2
DATA -2,-2
DATA 0,-1
DATA 2,-3


DIM X(N+1),Y(N+1) !頂点の座標
FOR i=1 TO N
   READ X(i),Y(i)
NEXT i

LET X(N+1)=X(1) !X[n+1]=X[1]で閉じる
LET Y(N+1)=Y(1)


LET Xmin=X(1) !4辺が軸と平行である制限を持った矩形（Axis Aligned Rectangle）を得る
LET Xmax=X(1)
LET Ymin=Y(1)
LET Ymax=Y(1)
FOR i=2 TO N
   IF X(i)<Xmin THEN LET Xmin=X(i)
   IF X(i)>Xmax THEN LET Xmax=X(i)
   IF Y(i)<Ymin THEN LET Ymin=Y(i)
   IF Y(i)>Ymax THEN LET Ymax=Y(i)
NEXT i
PRINT Xmin;Xmax; Ymin;Ymax !debug

IF Xmax-Xmin>Ymax-Ymin THEN LET SZ=Xmax-Xmin ELSE LET SZ=Ymax-Ymin

SET WINDOW Xmin-1,Xmin+SZ+1, Ymin-1,Ymin+SZ+1 !表示領域
DRAW grid


!AABB（Axis Aligned Bounding Box）を描く
SET LINE COLOR 2
PLOT LINES: Xmin,Ymin; !左下
PLOT LINES: Xmax,Ymin; !右下
PLOT LINES: Xmax,Ymax; !右上
PLOT LINES: Xmin,Ymax; !左上
PLOT LINES: Xmin,Ymin !閉じる
SET LINE COLOR 1


!多角形を描く
FOR i=1 TO N+1 !辺
   PLOT LINES: X(i),Y(i);
NEXT i
PLOT LINES
FOR i=1 TO N+1 !頂点
   DRAW disk WITH SCALE(0.05)*SHIFT(X(i),Y(i))
NEXT i


SET AREA COLOR 2

LET B=0 !辺上にある格子点の数、ただし、頂点は除く
FOR i=1 TO N
   LET s=ABS(X(i)-X(i+1)) !線分P1,P2上の格子点（P1,P2を除く）は、|x1-x2|と|y1-y2|の最大公約数-1
   LET t=ABS(Y(i)-Y(i+1))
   IF s=0 AND t=0 THEN !０個
   ELSE
      LET B=B+GCD(s,t)-1 !GCD(s,t)-1個

      LET w=GCD(s,t) !点を描く
      IF w>1 THEN
         LET dx=(X(i)-X(i+1))/w !増分
         LET dy=(Y(i)-Y(i+1))/w
         FOR j=1 TO w-1
            DRAW disk WITH SCALE(0.1)*SHIFT(X(i)-dx*j,Y(i)-dy*j)
         NEXT j
      END IF
   END IF
NEXT i
PRINT B !debug


SET AREA COLOR 4

LET A=0 !内部にある格子点の数
FOR PY=Ymin+1 TO Ymax-1 !対象の格子点(PX,PY)は、AABB内の点である
   FOR PX=Xmin+1 TO Xmax-1
      FOR i=1 TO N
         IF PX=X(i) AND PY=Y(i) THEN EXIT FOR
      NEXT i
      IF i>N THEN !頂点以外なら


         LET w=0 !多角形の内外判定　※辺上（頂点を含む）の点は内部にならない
         FOR i=1 TO N
            LET rx=PX-X(i)
            LET nx=PX-X(i+1)
            IF (rx>=0 AND nx<0) OR (rx<0 AND nx>=0) THEN
               LET ry=PY-Y(i)
               LET dx=X(i+1)-X(i)
               LET dy=Y(i+1)-Y(i)
               IF rx*dy<ry*dx THEN LET w=w+1 ELSE LET w=w-1 !※rx*dy=ry*dxなら、辺上となる
            END IF
         NEXT i
         !!PRINT PX;PY; w !debug
         IF w<>0 THEN !０なら外

            LET A=A+1
            DRAW disk WITH SCALE(0.1)*SHIFT(PX,PY) !点を描く

         END IF


      END IF
   NEXT PX
NEXT PY
PRINT A !debug


PRINT "面積="; A+(B+N)/2-1

END


EXTERNAL FUNCTION GCD(a,b) !最大公約数
DO WHILE b<>0
   LET t=b
   LET b=MOD(a,b)
   LET a=t
LOOP
LET GCD=a
END FUNCTION
 		 

戻る