平面射影変換

 投稿者:しばっち  投稿日:2015年 5月30日(土)22時16分41秒
  !' 平面射影変換
FILE GETNAME F$,"BMP,JPG,GIF,PNGファイル|*.BMP;*.JPG;*.GIF;*.PNG"
IF F$="" THEN STOP
CALL PICTURELOAD(F$,XSIZE,YSIZE)
DIM M(0 TO XSIZE-1,0 TO YSIZE-1),A(8,8),D(8)
ASK PIXEL ARRAY(0,0) M
!'入力画像
LET XA=XSIZE*.3 !'任意の四角形 左上(XA,YA)-右上(XB,YB)-右下(XC,YC)-左下(XD,YD)
LET YA=YSIZE*.1
LET XB=XSIZE*.8
LET YB=YSIZE*.05
LET XC=XSIZE*.9
LET YC=YSIZE*.8
LET XD=XSIZE*.1
LET YD=YSIZE*.9
FOR Y=0 TO YSIZE-1
   FOR X=0 TO XSIZE-1
      IF AREA4(XA,YA,XB,YB,XC,YC,XD,YD,X,Y)=0 THEN CALL PSET(X,Y,0,0,0)
   NEXT X
NEXT Y
!'出力画像
LET X1=XSIZE*.1 !'任意の四角形 左上(X1,Y1)- 右上(X2,Y2)-右下(X3,Y3)-左下(X4,Y4)
LET Y1=YSIZE*.3
LET X2=XSIZE*.9
LET Y2=YSIZE*.1
LET X3=XSIZE*.7
LET Y3=YSIZE*.8
LET X4=XSIZE*.2
LET Y4=YSIZE*.9
!'INPUT  PROMPT "ROTATE (0 - 3)=":N
!'CALL ROTATE(N,X1,X2,X3,X4) !'座標回転 (X1,Y1)-(X2,Y2)-(X3,Y3)-(X4,Y4) → (X2,Y2)-(X3,Y3)-(X4,Y4)-(X1,Y1) →
!'CALL ROTATE(N,Y1,Y2,Y3,Y4) !'         (X3,Y3)-(X4,Y4)-(X1,Y1)-(X2,Y2) → (X4,Y4)-(X1,Y1)-(X2,Y2)-(X3,Y3)
WAIT DELAY 2
CLEAR
SET COLOR COLORINDEX (1,0,0)
PLOT AREA : X1,Y1;X2,Y2;X3,Y3;X4,Y4
!'  平面射影変換
!'  X'=(A*X+B*Y+C)/(G*X+H*Y+1)
!'  Y'=(D*X+E*F+C)/(G*X+H*Y+1)

!'  A*X+B*Y+C-X'*(G*X+H*Y)=X'
!'  D*X+E*Y+F-Y'*(G*X+H*Y)=Y'
!' 八元連立方程式
!'(A)(XA, YA, 1, 0, 0, 0,-X1*XA,-X1*YA) (X1)
!'(B)(XB, YB, 1, 0, 0, 0,-X2*XB,-X2*YB) (X2)
!'(C)(XC, YC, 1, 0, 0, 0,-X3*XC,-X3*YC) (X3)
!'(D)(XD, YD, 1, 0, 0, 0,-X4*XD,-X4*YD) (X4)
!'(E)( 0,  0, 0,XA,YA, 1,-Y1*XA,-Y1*YA)=(Y1)
!'(F)( 0,  0, 0,XB,YB, 1,-Y2*XB,-Y2*YB) (Y2)
!'(G)( 0,  0, 0,XC,YC, 1,-Y3*XC,-Y3*YC) (Y3)
!'(H)( 0,  0, 0,XD,YD, 1,-Y4*XD,-Y4*YD) (Y4)
LET A(1,1)=XA
LET A(1,2)=YA
LET A(1,3)=1
LET A(1,7)=-X1*XA
LET A(1,8)=-X1*YA
LET D(1)=X1
LET A(2,1)=XB
LET A(2,2)=YB
LET A(2,3)=1
LET A(2,7)=-X2*XB
LET A(2,8)=-X2*YB
LET D(2)=X2
LET A(3,1)=XC
LET A(3,2)=YC
LET A(3,3)=1
LET A(3,7)=-X3*XC
LET A(3,8)=-X3*YC
LET D(3)=X3
LET A(4,1)=XD
LET A(4,2)=YD
LET A(4,3)=1
LET A(4,7)=-X4*XD
LET A(4,8)=-X4*YD
LET D(4)=X4
LET A(5,4)=XA
LET A(5,5)=YA
LET A(5,6)=1
LET A(5,7)=-Y1*XA
LET A(5,8)=-Y1*YA
LET D(5)=Y1
LET A(6,4)=XB
LET A(6,5)=YB
LET A(6,6)=1
LET A(6,7)=-Y2*XB
LET A(6,8)=-Y2*YB
LET D(6)=Y2
LET A(7,4)=XC
LET A(7,5)=YC
LET A(7,6)=1
LET A(7,7)=-Y3*XC
LET A(7,8)=-Y3*YC
LET D(7)=Y3
LET A(8,4)=XD
LET A(8,5)=YD
LET A(8,6)=1
LET A(8,7)=-Y4*XD
LET A(8,8)=-Y4*YD
LET D(8)=Y4
MAT A=INV(A)
MAT D=A*D
LET A1=D(1)
LET B1=D(2)
LET C1=D(3)
LET D1=D(4)
LET E1=D(5)
LET F1=D(6)
LET G1=D(7)
LET H1=D(8)
FOR Y=0 TO YSIZE-1
   FOR X=0 TO XSIZE-1
   !' X,Yについて解く(逆変換)
   !'  X'=(A*X+B*Y+C)/(G*X+H*Y+1)
   !'  Y'=(D*X+E*Y+F)/(G*X+H*Y+1)

   !'  A*X+B*Y+C-X'(G*X+H*Y+1)=0
   !'  D*X+E*Y+F-Y'(G*X+H*Y+1)=0
   !'  X,Yについてまとめる
   !'  X(A-G*X')+Y(B-X'*H)=X'-C
   !'  X(D-G*Y')+Y(E-Y'*H)=Y'-F
      IF AREA4(X1,Y1,X2,Y2,X3,Y3,X4,Y4,X,Y)<>0 THEN
         CALL SOLVE(A1-G1*X,B1-X*H1,X-C1,D1-G1*Y,E1-Y*H1,Y-F1,XA,YA)
         IF XA>=0 AND YA>=0 AND XA<=XSIZE-1 AND YA<=YSIZE-1 THEN
            CALL RGB(M(XA,YA),RR,GG,BB)
            CALL PSET(X,Y,RR,GG,BB)
         END IF
      END IF
   NEXT X
NEXT Y
END

EXTERNAL FUNCTION AREA3(X1,Y1,X2,Y2,X3,Y3,PX,PY)
LET T=TRIANGLE(X1,Y1,X2,Y2,X3,Y3)
LET A=TRIANGLE(X1,Y1,X2,Y2,PX,PY)
LET B=TRIANGLE(X2,Y2,X3,Y3,PX,PY)
LET C=TRIANGLE(X1,Y1,X3,Y3,PX,PY)
IF ABS(A+B+C-T)<3 THEN LET AREA3=-1 ELSE LET AREA3=0
END FUNCTION

EXTERNAL FUNCTION AREA4(X1,Y1,X2,Y2,X3,Y3,X4,Y4,PX,PY)
LET A=AREA3(X1,Y1,X2,Y2,X3,Y3,PX,PY)
LET B=AREA3(X1,Y1,X4,Y4,X3,Y3,PX,PY)
IF A<>0 OR B<>0 THEN LET AREA4=-1 ELSE LET AREA4=0
END FUNCTION

EXTERNAL SUB SOLVE(A1,B1,C1,A2,B2,C2,X,Y)
!'A1*X+B1*Y=C1
!'A2*X+B2*Y=C2
LET DD=A1*B2-A2*B1
IF DD=0 THEN
   LET X=-1
   LET Y=-1
ELSE
   LET X=(C1*B2-C2*B1)/DD
   LET Y=(A1*C2-A2*C1)/DD
END IF
END SUB

EXTERNAL FUNCTION TRIANGLE(X1,Y1,X2,Y2,X3,Y3)
LET TRIANGLE=ABS(X1*Y2+X2*Y3+X3*Y1-X2*Y1-X3*Y2-Y3*X1)/2
END FUNCTION

EXTERNAL SUB RGB(X,R,G,B)
LET B=MOD(INT(X/65536),256)
LET G=MOD(INT(X/256),256)
LET R=MOD(X,256)
END SUB

EXTERNAL SUB GINIT(XSIZE,YSIZE)
SET BITMAP SIZE XSIZE,YSIZE
CLEAR
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
SET WINDOW 0,XSIZE-1,YSIZE-1,0
END SUB

EXTERNAL SUB PSET(X,Y,R,G,B)
SET COLOR COLORINDEX(R/255,G/255,B/255)
PLOT POINTS:X,Y
END SUB

EXTERNAL SUB PICTURELOAD(N$,XSIZE,YSIZE)
CLEAR
SET COLOR MODE "NATIVE"
SET POINT STYLE 1
GLOAD N$
LET XSIZE=PIXELX(1)
LET YSIZE=PIXELY(1)
SET BITMAP SIZE XSIZE,YSIZE
SET WINDOW 0,XSIZE-1,YSIZE-1,0
END SUB

EXTERNAL  SUB ROTATE(N,A,B,C,D)
FOR I=1 TO N
   LET T=B
   LET B=A
   LET A=D
   LET D=C
   LET C=T
NEXT I
END SUB
 

戻る