|
|
!'離散WAVELET変換(多重解像度解析)
OPTION BASE 0
PUBLIC NUMERIC XSIZE,YSIZE,DAUBECHIES,P(7),Q(7)
FILE GETNAME F$, "BMP,JPG,GIFファイル|*.BMP;*.JPG;*.GIF"
IF F$="" THEN STOP
CALL PICTURELOAD(F$,XSIZE,YSIZE)
LET LEV=2 !'等分割 (4^LEV 分割)
!'16分割時の帯域分割(L:低周波 H:高周波)
!'LLLL LLHL HLLL HLHL
!'LLLH LLHH HLLH HLHH
!'LHLL LHHL HHLL HHHL
!'LHLH LHHH HHLH HHHH
LET DAUBECHIES=2
LET SIZEX=INT((XSIZE+2^LEV-1)/2^LEV)*2^LEV !'補正
LET SIZEY=INT((YSIZE+2^LEV-1)/2^LEV)*2^LEV
DIM R1(SIZEX,SIZEY),G1(SIZEX,SIZEY),B1(SIZEX,SIZEY) !'大量のメモリーが必要
DIM R2(SIZEX,SIZEY),G2(SIZEX,SIZEY),B2(SIZEX,SIZEY) !'再帰呼び出しも行っている
SELECT CASE DAUBECHIES
!'ドペシィーの数列
CASE 1
LET P(0) = SQR(2)/2
LET P(1) = SQR(2)/2
CASE 2
LET P(0) = 0.482962913144534143374871599864
LET P(1) = 0.836516303737807905575293780917
LET P(2) = 0.22414386804201338102597276224
LET P(3) = -0.129409522551260381174449418812
CASE 3
LET P(0) = 0.332670552950082615997608286301
LET P(1) = 0.806891509311092576493086842682
LET P(2) = 0.459877502118491570095951630688
LET P(3) = -0.135011020010254588694583301022
LET P(4) = -0.854412738820266616927155548839E-1
LET P(5) = 0.352262918857095366023408204452E-1
END SELECT
FOR I=0 TO 2*DAUBECHIES-1
LET Q(I)=(-1)^I*P(2*DAUBECHIES-I-1)
NEXT I
PRINT "画像データ読み込み"
FOR X=0 TO XSIZE-1
FOR Y=0 TO YSIZE-1
CALL GETPOINT(X,Y,R1(X,Y),G1(X,Y),B1(X,Y))
NEXT Y
NEXT X
PRINT "ウェーブレット変換"
CALL WAVELET(LEV,0,0,SIZEX,SIZEY,R1,G1,B1,R2,G2,B2)
PRINT "マウスでクリアしたい帯域を左クリック"
PRINT "右クリックで再構成します"
DO
DO
MOUSE POLL XX,YY,L,R
LOOP UNTIL L<>0 OR R<>0
IF L<>0 THEN CALL CLEARAREA(LEV,XX,YY,SIZEX,SIZEY,R1,G1,B1) !'クリックされた帯域内をクリア
LOOP UNTIL R<>0
!'再構成(逆ウェーブレット変換)
LET SIZEX=INT((XSIZE+2^LEV-1)/2^LEV)*2^LEV !'補正
LET SIZEY=INT((YSIZE+2^LEV-1)/2^LEV)*2^LEV
MAT R2=ZER !'データクリア
MAT G2=ZER
MAT B2=ZER
PRINT "再構成"
CALL INVWAVELET(LEV,0,0,SIZEX,SIZEY,R1,G1,B1,R2,G2,B2) !'逆ウェーブレット変換
END
EXTERNAL SUB WAVELET(N,XS,YS,SIZEX,SIZEY,R1(,),G1(,),B1(,),R2(,),G2(,),B2(,))
IF N>0 THEN
FOR Y=0 TO SIZEY-1 !'X方向に帯域分割
FOR X=0 TO SIZEX-1 STEP 2
LET R2(XS+X/2,Y+YS)=0
LET G2(XS+X/2,Y+YS)=0
LET B2(XS+X/2,Y+YS)=0
LET R2(XS+X/2+SIZEX/2,Y+YS)=0
LET G2(XS+X/2+SIZEX/2,Y+YS)=0
LET B2(XS+X/2+SIZEX/2,Y+YS)=0
FOR NN=0 TO DAUBECHIES-1
LET R2(XS+X/2,Y+YS)=R2(XS+X/2,Y+YS)+P(2*NN)*R1(X+XS,Y+YS)+P(2*NN+1)*R1(X+1+XS,Y+YS) !'低周波
LET G2(XS+X/2,Y+YS)=G2(XS+X/2,Y+YS)+P(2*NN)*G1(X+XS,Y+YS)+P(2*NN+1)*G1(X+1+XS,Y+YS)
LET B2(XS+X/2,Y+YS)=B2(XS+X/2,Y+YS)+P(2*NN)*B1(X+XS,Y+YS)+P(2*NN+1)*B1(X+1+XS,Y+YS)
LET R2(XS+X/2+SIZEX/2,Y+YS)=R2(XS+X/2+SIZEX/2,Y+YS)+Q(2*NN)*R1(X+XS,Y+YS)+Q(2*NN+1)*R1(X+1+XS,Y+YS) !'高周波
LET G2(XS+X/2+SIZEX/2,Y+YS)=G2(XS+X/2+SIZEX/2,Y+YS)+Q(2*NN)*G1(X+XS,Y+YS)+Q(2*NN+1)*G1(X+1+XS,Y+YS)
LET B2(XS+X/2+SIZEX/2,Y+YS)=B2(XS+X/2+SIZEX/2,Y+YS)+Q(2*NN)*B1(X+XS,Y+YS)+Q(2*NN+1)*B1(X+1+XS,Y+YS)
NEXT NN
CALL PSET(XS+X/2,Y+YS,R2(X/2+XS,Y+YS),G2(X/2+XS,Y+YS),B2(X/2+XS,Y+YS))
CALL PSET(XS+X/2+SIZEX/2,Y+YS,R2(X/2+SIZEX/2+XS,Y+YS),G2(X/2+SIZEX/2+XS,Y+YS),B2(X/2+SIZEX/2+XS,Y))
NEXT X
NEXT Y
FOR X=0 TO SIZEX-1
FOR Y=0 TO SIZEY-1 STEP 2 !'Y方向に帯域分割
LET R1(XS+X,YS+Y/2)=0
LET G1(XS+X,YS+Y/2)=0
LET B1(XS+X,YS+Y/2)=0
LET R1(XS+X,YS+Y/2+SIZEY/2)=0
LET G1(XS+X,YS+Y/2+SIZEY/2)=0
LET B1(XS+X,YS+Y/2+SIZEY/2)=0
FOR NN=0 TO DAUBECHIES-1
LET R1(XS+X,YS+Y/2)=R1(XS+X,YS+Y/2)+P(2*NN)*R2(XS+X,Y+YS)+P(2*NN+1)*R2(XS+X,Y+YS+1) !'低周波
LET G1(XS+X,YS+Y/2)=G1(XS+X,YS+Y/2)+P(2*NN)*G2(XS+X,Y+YS)+P(2*NN+1)*G2(XS+X,Y+YS+1)
LET B1(XS+X,YS+Y/2)=B1(XS+X,YS+Y/2)+P(2*NN)*B2(XS+X,Y+YS)+P(2*NN+1)*B2(XS+X,Y+YS+1)
LET R1(XS+X,YS+Y/2+SIZEY/2)=R1(XS+X,YS+Y/2+SIZEY/2)+Q(2*NN)*R2(XS+X,Y+YS)+Q(2*NN+1)*R2(XS+X,Y+YS+1) !'高周波
LET G1(XS+X,YS+Y/2+SIZEY/2)=G1(XS+X,YS+Y/2+SIZEY/2)+Q(2*NN)*G2(XS+X,Y+YS)+Q(2*NN+1)*G2(XS+X,Y+YS+1)
LET B1(XS+X,YS+Y/2+SIZEY/2)=B1(XS+X,YS+Y/2+SIZEY/2)+Q(2*NN)*B2(XS+X,Y+YS)+Q(2*NN+1)*B2(XS+X,Y+YS+1)
NEXT NN
CALL PSET(XS+X,YS+Y/2,R1(XS+X,Y/2+YS),G1(XS+X,Y/2+YS),B1(XS+X,Y/2+YS))
CALL PSET(XS+X,YS+Y/2+SIZEY/2,R1(XS+X,Y/2+SIZEY/2+YS),G1(XS+X,Y/2+SIZEY/2+YS),B1(XS+X,Y/2+SIZEY/2+YS))
NEXT Y
NEXT X
CALL WAVELET(N-1,0,0,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL WAVELET(N-1,SIZEX/2,0,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL WAVELET(N-1,0,SIZEY/2,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL WAVELET(N-1,SIZEX/2,SIZEY/2,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
END IF
END SUB
EXTERNAL SUB INVWAVELET(N,XS,YS,SIZEX,SIZEY,R1(,),G1(,),B1(,),R2(,),G2(,),B2(,))
IF N>0 THEN
CALL INVWAVELET(N-1,0,0,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL INVWAVELET(N-1,SIZEX/2,0,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL INVWAVELET(N-1,0,SIZEY/2,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
CALL INVWAVELET(N-1,SIZEX/2,SIZEY/2,SIZEX/2,SIZEY/2,R1,G1,B1,R2,G2,B2)
FOR Y=0 TO SIZEY/2-1
FOR X=0 TO SIZEX-1
LET R2(X+XS,Y*2+YS)=0
LET G2(X+XS,Y*2+YS)=0
LET B2(X+XS,Y*2+YS)=0
LET R2(X+XS,Y*2+1+YS)=0
LET G2(X+XS,Y*2+1+YS)=0
LET B2(X+XS,Y*2+1+YS)=0
FOR NN=0 TO DAUBECHIES-1
LET R2(X+XS,Y*2+YS)=R2(X+XS,Y*2+YS)+P(2*NN)*R1(X+XS,Y+YS)+Q(2*NN)*R1(X+XS,Y+SIZEY/2+YS)
LET G2(X+XS,Y*2+YS)=G2(X+XS,Y*2+YS)+P(2*NN)*G1(X+XS,Y+YS)+Q(2*NN)*G1(X+XS,Y+SIZEY/2+YS)
LET B2(X+XS,Y*2+YS)=B2(X+XS,Y*2+YS)+P(2*NN)*B1(X+XS,Y+YS)+Q(2*NN)*B1(X+XS,Y+SIZEY/2+YS)
LET R2(X+XS,Y*2+1+YS)=R2(X+XS,Y*2+1+YS)+P(2*NN+1)*R1(X+XS,Y+YS)+Q(2*NN+1)*R1(X+XS,Y+SIZEY/2+YS)
LET G2(X+XS,Y*2+1+YS)=G2(X+XS,Y*2+1+YS)+P(2*NN+1)*G1(X+XS,Y+YS)+Q(2*NN+1)*G1(X+XS,Y+SIZEY/2+YS)
LET B2(X+XS,Y*2+1+YS)=B2(X+XS,Y*2+1+YS)+P(2*NN+1)*B1(X+XS,Y+YS)+Q(2*NN+1)*B1(X+XS,Y+SIZEY/2+YS)
NEXT NN
CALL PSET(X+XS,Y*2+YS,R2(X+XS,Y*2+YS),G2(X+XS,Y*2+YS),B2(X+XS,Y*2+YS))
CALL PSET(X+XS,Y*2+1+YS,R2(X+XS,Y*2+1+YS),G2(X+XS,Y*2+1+YS),B2(X+XS,Y*2+1+YS))
NEXT X
NEXT Y
FOR X=0 TO SIZEX/2-1
FOR Y=0 TO SIZEY-1
LET R1(X*2+XS,Y+YS)=0
LET G1(X*2+XS,Y+YS)=0
LET B1(X*2+XS,Y+YS)=0
LET R1(X*2+1+XS,Y+YS)=0
LET G1(X*2+1+XS,Y+YS)=0
LET B1(X*2+1+XS,Y+YS)=0
FOR NN=0 TO DAUBECHIES-1
LET R1(X*2+XS,Y+YS)=R1(X*2+XS,Y+YS)+P(2*NN)*R2(X+XS,Y+YS)+Q(2*NN)*R2(X+SIZEX/2+XS,Y+YS)
LET G1(X*2+XS,Y+YS)=G1(X*2+XS,Y+YS)+P(2*NN)*G2(X+XS,Y+YS)+Q(2*NN)*G2(X+SIZEX/2+XS,Y+YS)
LET B1(X*2+XS,Y+YS)=B1(X*2+XS,Y+YS)+P(2*NN)*B2(X+XS,Y+YS)+Q(2*NN)*B2(X+SIZEX/2+XS,Y+YS)
LET R1(X*2+1+XS,Y+YS)=R1(X*2+1+XS,Y+YS)+P(2*NN+1)*R2(X+XS,Y+YS)+Q(2*NN+1)*R2(X+SIZEX/2+XS,Y+YS)
LET G1(X*2+1+XS,Y+YS)=G1(X*2+1+XS,Y+YS)+P(2*NN+1)*G2(X+XS,Y+YS)+Q(2*NN+1)*G2(X+SIZEX/2+XS,Y+YS)
LET B1(X*2+1+XS,Y+YS)=B1(X*2+1+XS,Y+YS)+P(2*NN+1)*B2(X+XS,Y+YS)+Q(2*NN+1)*B2(X+SIZEX/2+XS,Y+YS)
NEXT NN
CALL PSET(X*2+XS,Y+YS,R1(X*2+XS,Y+YS),G1(X*2+XS,Y+YS),B1(X*2+XS,Y+YS))
CALL PSET(X*2+1+XS,Y+YS,R1(X*2+1+XS,Y+YS),G1(X*2+1+XS,Y),B1(X*2+1+XS,Y+YS))
NEXT Y
NEXT X
END IF
END SUB
EXTERNAL SUB CLEARAREA(N,X,Y,SIZEX,SIZEY,R1(,),G1(,),B1(,))
LET XX=INT(SIZEX/N/2)
LET YY=INT(SIZEY/N/2)
LET XS=INT(X/XX)*XX
LET YS=INT(Y/YY)*YY
CALL CLEAR(XS,YS,XS+XX-1,YS+YY-1,R1,G1,B1)
END SUB
EXTERNAL SUB CLEAR(XS,YS,XE,YE,R(,),G(,),B(,))
FOR XX=XS TO XE
FOR YY=YS TO YE
LET R(XX,YY)=0 !'データクリア
LET G(XX,YY)=0
LET B(XX,YY)=0
IF XX<=XSIZE-1 AND YY<=YSIZE-1 THEN CALL PSET(XX,YY,0,0,0) !'画像もクリア
NEXT YY
NEXT XX
END SUB
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 GETPOINT(X,Y,R,G,B)
ASK PIXEL VALUE (X,Y) C
CALL RGB(C,R,G,B)
END SUB
EXTERNAL SUB GINIT(XSIZE,YSIZE)
SET BITMAP SIZE XSIZE,YSIZE
SET COLOR MODE "NATIVE"
CLEAR
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(MOD(INT(R)+256,256)/255,MOD(INT(G)+256,256)/255,MOD(INT(B)+256,256)/255)
PLOT POINTS: X , Y
END SUB
EXTERNAL SUB PICTURELOAD(N$,XSIZE,YSIZE)
SET COLOR MODE "NATIVE"
CLEAR
SET POINT STYLE 1
GLOAD N$
LET XSIZE=PIXELX(1)+1
LET YSIZE=PIXELY(1)+1
SET BITMAP SIZE XSIZE,YSIZE
SET WINDOW 0,XSIZE-1,YSIZE-1,0
END SUB
|
|