遊び tokidoki 仕事

数学と音楽と教育と遊び

| おしごと - きょういく - がくせい - ゼミ - イベント | すうがく - おんがく - 数理音楽 - DTM - かがく - scratch
| Art - photo - おきにー - Tips - ものもう - あれこれ | About - Top

超解像シミュレーション

痛すぎる卒論指導も終盤。
ひょんなことから一人のゼミ生には
超解像の原理についてまとめてもらっている。
原理は関数解析的な話になるがそれはさておき、
実際にそれができるものなのか試したくなった。
10進BASICで簡単なシミュレーションをやってみたので
そのプログラムを公開。
ただ、かなり汚い書き方なのでご容赦を。

【使い方】
とにかく実行すると超解像を順次繰り返します。
画像をファイルとして保存したいときは
!GSAVE
となっている行の ! を外してください。

こんなかんじ↓






10進BASICプログラム→[superResolution1D.bas]

REM [superResolution1D.bas]
REM 超解像のシミュレーション
REM 帯域制限された1D波形を回復するシミュレーション
REM

REM FFTおよびIFFTサブルーチンを使用するために必要な宣言
OPTION ARITHMETIC complex  ! 複素数モードに設定
DECLARE EXTERNAL SUB FFT   ! FFT関数
DECLARE EXTERNAL SUB IFFT   ! FFT関数

REM FFTのための配列および定数宣言
LET dn = 1024*16           ! データ総数
DIM d(dn)               ! 入出力バッファ用配列
DIM g(dn)               ! 帯域制限領域周波数用配列

REM 【元の関数の定義】
DEF f(x)=MIN(MAX(2-(x/.0032)^2,0),1.8)

LET  dt = .001   ! サンプリング間隔
LET  XB=.4   !空間幅(割合)
LET  FB=.4   !周波数制限(割合)
LET  ht=2.2 !元の表示windowの高さ
LET  wd=dt*dn/64 !元の表示wondowの幅
LET  fht=250 !周波数領域の表示windowの高さ
LET  win=.01
LET  boxx=2
LET  boxy=3
LET  bmp=300
SET bitmap SIZE bmp*(1+win*2)*boxx,bmp*(1+win*2)*boxy
LET  boxz=1/MAX(boxx,boxy)
SET AREA COLOR 0

REM 【元の関数の表示】====================================================
SET VIEWPORT win,boxz-win , 1-boxz+win,1-win
SET WINDOW -wd*1.1, wd*1.1, -ht, ht
DRAW axes0
LET t = 0
FOR i = 1 TO dn
   LET d(i) = f(i/dn-.5) 
   PLOT LINES: (i-dn/2)*dt, d(i);
   LET t = t + dt
NEXT i
PLOT TEXT ,AT -wd,ht*.9: "元の信号"

REM 【周波数リストの作成】(後半半分のデータは無視する)
dim freq(dn/2+1)
let maxFreq = 1/dt / 2  ! 最大周波数
for i=1 to dn/2+1
   let freq(i) = maxFreq * (i-1) / (dn/2)
NEXT i

REM 【帯域制限定数】
LET tauX=wd*XB/dt !空間制限
LET  tau=INT(maxFreq*FB)      !周波数制限


REM 【元信号の周波数スペクトラム】================================================
SET VIEWPORT boxz+win,boxx*boxz-win , 1-boxz+win,1-win
SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht
DRAW axes(maxFreq/20, fht/2)
LET t = 0
call FFT(d)
PLOT LINES
FOR i = 1 TO dn/2+1
   PLOT LINES: freq(i), re(d(i));
   IF i>tau THEN LET  d(i)=0
   LET  g(i)=d(i)
NEXT i
PLOT TEXT ,AT maxFreq*.02,fht*.9: "元の信号のフーリエ変換"

LET  c=0    !超解像の適用回数カウンタ
DO WHILE(c<1000)
 
   REM 【周波数帯域制限】====================================================
   SET VIEWPORT boxz+win,boxx*boxz-win , 1-boxz*2+win,1-boxz-win
   SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht
   PLOT AREA : -maxFreq,-fht;-maxFreq,fht; maxFreq,fht; maxFreq,-fht
   DRAW axes(maxFreq/20, fht/2)
   PLOT LINES
   FOR i = 1 TO dn/2+1
      IF i<=tau THEN LET d(i)=g(i)
      PLOT LINES: freq(i), re(d(i));
   NEXT i
   PLOT LINES
   SET LINE COLOR 2
   PLOT LINES: freq(tau),-fht*.9; freq(tau),fht*.9
   PLOT TEXT ,AT freq(tau), -fht*.4: "τ"
   PLOT TEXT ,AT maxFreq*.02,fht*.9: "帯域制限信号との合成"
    
   REM 【フーリエ逆変換して表示】============================================
   SET VIEWPORT win,boxz-win , 1-boxz*2+win,1-boxz-win
   SET WINDOW -wd*1.1, wd*1.1, -ht, ht
   PLOT AREA : -wd*1.1,-ht;-wd*1.1,ht; wd*1.1,ht; wd*1.1,-ht
   DRAW axes0
   LET t = 0
   PLOT lines
   CALL IFFT(d)
   SET LINE COLOR 1
   FOR i = 1 TO dn
      PLOT LINES: (i-dn/2)*dt, d(i);
      LET t = t + dt
   NEXT i
   PLOT TEXT ,AT -wd,ht*.9: "超解像"&STR$(c)&"回目"
    
   REM 【空間帯域制限】======================================================
   SET VIEWPORT win,boxz-win , win,1-boxz*2-win
   SET WINDOW -wd*1.1, wd*1.1, -ht, ht
   PLOT AREA : -wd*1.1,-ht;-wd*1.1,ht; wd*1.1,ht; wd*1.1,-ht
   DRAW axes0
   LET t = 0
   PLOT LINES
   FOR i = 1 TO dn
      IF ABS(i-dn/2)>tauX THEN LET  d(i)=0
      PLOT LINES: (i-dn/2)*dt, d(i);
      LET t = t + dt
   NEXT i
   PLOT LINES
   SET LINE COLOR 2
   PLOT LINES: tauX*dt,-ht*.9; tauX*dt,ht*.9
   PLOT LINES: -tauX*dt,-ht*.9; -tauX*dt,ht*.9
   PLOT TEXT ,AT -wd,ht*.9: "空間制限"
   PLOT TEXT ,AT -tauX*dt, -ht*.2: "-T"
   PLOT TEXT ,AT tauX*dt*1.1, -ht*.2: "T"
    
   REM 【超解像周波数スペクトラム】================================================
   SET VIEWPORT boxz+win,boxx*boxz-win , win,1-boxz*2-win
   SET WINDOW -maxFreq/100, maxFreq/10, -fht, fht
   PLOT AREA : -maxFreq,-fht;-maxFreq,fht; maxFreq,fht; maxFreq,-fht
   DRAW axes(maxFreq/20, fht/2)
   LET t = 0
   call FFT(d)
   PLOT LINES
   SET LINE COLOR 1
   FOR i = 1 TO dn/2+1
      PLOT LINES: freq(i), re(d(i));
   NEXT i
   PLOT TEXT ,AT maxFreq*.02,fht*.9: "空間制限信号のフーリエ変換"
    
   REM 【画像ファイルとして出力】==========================================
   LET  sw=0
   IF MOD(c,50)=0 THEN
      LET  sw=1
   ELSEIF (c<31 AND c>9) AND MOD(c,10)=0 THEN
      LET  sw=1
   ELSEIF c<10 THEN
      LET  sw=1
   END IF
   IF sw=1 THEN
      LET  f$="spr"&right$("000"&STR$(c),3)&".bmp"
      !gsave f$, "4bit"
   END IF
   LET  c=c+1
LOOP

END

REM ======================================================================
REM  [1D FFT/IFFT]
REM   元データ x(,) → FFT/IFFT → 変換後データ x(,)
REM   ユニタリ変換になるようFFT/IFFT両方を1/N
REM   Fk = 1/N * Σ[n=1...N] (fn * W^( (k-1)*(n-1)) )   :FFT
REM   fk = 1/N * Σ[n=1...N] (Fn * W^(-(k-1)*(n-1)) )  :IFFT
REM   W = exp(-��(-1) * 2π/N)
REM   配列添字は 1 から出発

REM 【1D FFT】
REM θ = -2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB FFT(x())
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x)
LET th = -2 * PI / sx
CALL FFTMAIN(x, th)
END SUB

REM 【1D IFFT】
REM θ = 2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB IFFT(x())
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x)
LET th = 2 * PI / sx
CALL FFTMAIN(x, th)
MAT x = (1/sx) * x
END SUB

REM [バタフライ演算]
REM 配列サイズを半分にしながら再帰呼び出しする
REM wk= exp( (k-1)θ )
REM x(k)    --- *1  -- x0(k)
REM         × 
REM x(k+N/2)--- *wk -- x1(k)

EXTERNAL SUB FFTMAIN(x(), th)
OPTION ARITHMETIC complex
LET sx = SIZE(x)
IF sx>1 THEN
   LET sx2 = sx/2
   DIM x0(sx2), x1(sx2)
   FOR k=1 TO sx2
      LET x0(k) = x(k) + x(k+sx2)
      LET wk = EXP( complex(0, th * (k-1)) )
      LET x1(k) = wk * ( x(k) - x(k+sx2) )
   NEXT k
   LET th0 = 2 * th
   CALL FFTMAIN(x0, th0)
   LET th1 = 2 * th
   CALL FFTMAIN(x1, th1)
   FOR k=1 TO sx2
      LET x(2*k-1) = x0(k)
      LET x(2*k) = x1(k)
   NEXT k
END IF
END SUB

2次元FFTで画像のフーリエ変換

引き続き、卒論用のプログラムシリーズ。
今回は、画像のフーリエ変換を表示する10進BASICプログラム。
たとえば円環領域のフーリエ変換は

といった感じに出力。


10進BASICプログラム→[FFT2D.bas]

REM [FFT2D.bas]
REM 2D FFT/IFFT
REM FFT/IFFTは外部関数として使う

OPTION ARITHMETIC complex
SET TEXT FONT "",6
SET AXIS COLOR 2

DECLARE EXTERNAL SUB FFT    ! FFT
DECLARE EXTERNAL SUB IFFT   ! IFFT

! 配列サイズ dn は必ず2^nで
LET dn = 256            ! データ総数
LET  dn2 = dn/2
LET  win=.02            ! windowの余白
DIM bmp(dn,dn)          ! 画像用配列 dn×dn
SET bitmap SIZE dn*2*(1+win*2),dn*(1+win*2)

REM 【元の関数の定義】
DEF f(x,y)=SGN(x^2+y^2-.45)-SGN(x^2+y^2-.5)-.1

LET  dt = 1   ! サンプリング間隔
LET  ct= 10^2  !元画像のコントラスト
LET  fct= 10^1  !フーリエスペクトルのコントラスト

SET POINT COLOR 10
SET POINT STYLE 1

REM 【元の画像の表示】

SET VIEWPORT win,.5-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
FOR y=-dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET  g=f(x/dn2,y/dn2)
      LET bmp(x+dn2+1,y+dn2+1) = g
      LET  g=ATN(g*ct)/PI+.5
      SET COLOR MIX(10) g,g,g 
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

CALL FFT(bmp)

REM 【周波数スペクトラム】

SET VIEWPORT .5+win,1-win , win,.5-win
SET WINDOW -dn2, dn2, -dn2, dn2
FOR y= -dn2 TO dn2-1
   FOR x = -dn2 TO dn2-1
      LET g=ATN(ABS(bmp(MOD(x+dn,dn)+1,MOD(y+dn,dn)+1))*fct)/PI*2
      SET COLOR MIX(10) g,g,g 
      PLOT POINTS: x,y
   NEXT x
NEXT y
DRAW axes(dn/8,dn/8)

END

REM  [2D FFT/IFFT]
REM   元データ x(,) → FFT/IFFT → 変換後データ x(,)
REM   ユニタリ変換になるようFFT/IFFT両方を1/N
REM   Fk = 1/N * Σ[n=1...N] (fn * W^( (k-1)*(n-1)) )   :FFT
REM   fk = 1/N * Σ[n=1...N] (Fn * W^(-(k-1)*(n-1)) )  :IFFT
REM   W = exp(-��(-1) * 2π/N)
REM   配列添字は 1 から出発

REM 【2D FFT】
REM θ = -2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB FFT(x(,))
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = -2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM 【2D IFFT】
REM θ = 2π/N
REM W = exp(-��(-1)θ)

EXTERNAL SUB IFFT(x(,))
OPTION ARITHMETIC complex
DECLARE EXTERNAL SUB FFTMAIN
LET sx = SIZE(x,1)
DIM d(sx)
LET th = 2 * PI / sx
FOR t=1 TO sx
   FOR s=1 TO sx
      LET  d(s)=x(s,t)
   NEXT s
   CALL FFTMAIN(d, th)
   FOR s=1 TO sx
      LET  x(s,t)=d(s)
   NEXT s
NEXT t
MAT x = (1/sx) * x
FOR s=1 TO sx
   FOR t=1 TO sx
      LET  d(t)=x(s,t)
   NEXT t
   CALL FFTMAIN(d, th)
   FOR t=1 TO sx
      LET  x(s,t)=d(t)
   NEXT t
NEXT s
END SUB

REM [バタフライ演算]
REM 配列サイズを半分にしながら再帰呼び出しする
REM wk= exp( (k-1)θ )
REM x(k)    --- ×1  -- x0(k)
REM         × 
REM x(k+N/2)--- ×wk -- x1(k)

EXTERNAL SUB FFTMAIN(x(), th)
OPTION ARITHMETIC complex
LET sx = SIZE(x)
IF sx>1 THEN
   LET sx2 = sx/2
   DIM x0(sx2), x1(sx2)
   FOR k=1 TO sx2
      LET x0(k) = x(k) + x(k+sx2)
      LET wk = EXP( complex(0, th * (k-1)) )
      LET x1(k) = wk * ( x(k) - x(k+sx2) )
   NEXT k
   LET th0 = 2 * th
   CALL FFTMAIN(x0, th0)
   LET th1 = 2 * th
   CALL FFTMAIN(x1, th1)
   FOR k=1 TO sx2
      LET x(2*k-1) = x0(k)
      LET x(2*k) = x1(k)
   NEXT k
END IF
END SUB

CSVからBMPへ

調子に乗って、今度はCSVからBMPへ戻すプログラム。
3枚のCSV化した画像データを
EXCELシート上で「目に見える形で」フィルター処理をして、
それを再び10進BASICでBMPに戻すわけだ。


10進BASICプログラム→[csv2bmp.bas]

REM [csv2bmp.bas]
REM RGBに分解された3つのCSVファイルからBMP画像データを作る
REM
REM office2003のEXCELシート互換にするため
REM BMPの横幅は256pxまで、縦幅は65536pxまで
REM
REM [使用法]
REM (1)このプログラムと同じフォルダに加工後のRGBデータ
REM       ○○_rr_out.csv,○○_gg_out.csv,○○_bb_out.csv
REM ファイルを置く
REM (2)読み込むCSVファイル名はf$="○○"で指定する
REM (3)実行すると3つのCSVファイルを読み込み
REM       ○○_out.bmp
REM としてBMPファイルが同じフォルダに出力される

OPTION ARITHMETIC NATIVE
SET COLOR MODE "NATIVE"
LET  f$="neko"
!比較のためあらかじめ原画像も表示しておく
GLOAD f$ & ".bmp"
ASK PIXEL SIZE (0,0; 1,1) bw,bh
DIM bc(bw,bh), rr(bw,bh), gg(bw,bh), bb(bw,bh), bf(bw,bh)
DIM dr(bw), dg(bw), db(bw)
ASK PIXEL ARRAY (0,1) bc
SET BITMAP SIZE bw,bh

OPEN #1:NAME f$ & "_rr_out.csv", RECTYPE INTERNAL
OPEN #2:NAME f$ & "_gg_out.csv", RECTYPE INTERNAL
OPEN #3:NAME f$ & "_bb_out.csv", RECTYPE INTERNAL
FOR y=1 TO bh
   MAT READ #1: dr
   MAT READ #2: dg
   MAT READ #3: db
   FOR x=1 TO bw
      LET  rr(x,y)=MAX(MIN(dr(x),255),0)
      LET  gg(x,y)=MAX(MIN(dg(x),255),0)
      LET  bb(x,y)=MAX(MIN(db(x),255),0)
   NEXT x
NEXT y
CLOSE #3
CLOSE #2
CLOSE #1

MAT gg=256*gg
MAT bb=65536*bb
MAT bf=rr+gg
MAT bf=bc+bb

!一瞬だけ変換後のBMPを表示してBMPファイルとして保存
MAT PLOT CELLS, IN 0,1; 1,0 :bf
GSAVE f$ & "_out.bmp"

!その後、変換前と後を並べて表示
SET BITMAP SIZE bw*2,bh*2
MAT PLOT CELLS, IN 0,1; 1/2,1/2 :bc
MAT PLOT CELLS, IN 0,1/2; 1/2,0 :bf

END

BMPからCSVへ

今年度卒業のゼミ生、卒論ではFourier変換の利用をテーマに進めている。
とりかかりが遅かったため、かなりピンチな状況だ。


何人かがデジタル画像処理をテーマにしていることもあって、
画像処理をブラックボックス化してしまわないで
教育的に目に見える形で処理させたい。
そこで、BMP画像をEXCELで処理できるように
BMP画像からRGBに分解してCSV形式で出力するプログラムを
(半分は趣味を兼ねて)10進BASICで作ったので、
ここに公開しておきます。


10進BASICプログラム→[bmp2csv.bas]

REM [bmp2csv.bas]
REM BMP画像データをRGBに分解し3つのCSVファイルとして出力する
REM
REM office2003のEXCELシート互換にするため
REM BMPの横幅は256pxまで、縦幅は65536pxまで
REM
REM [使用法]
REM (1)このプログラムと同じフォルダにBMP画像を置く
REM (2)読み込む画像ファイル名はf$="○○"で指定する
REM (3)実行するとCSV形式でRGB情報が
REM       ○○_rr.csv,○○_gg.csv,○○_bb.csv
REM として同じフォルダに出力される

OPTION ARITHMETIC NATIVE
SET COLOR MODE "NATIVE"
LET  f$="neko"
GLOAD f$ & ".bmp"
ASK PIXEL SIZE (0,0; 1,1) bw,bh
DIM bc(bw,bh), rr(bw,bh), gg(bw,bh), bb(bw,bh), bf(bw,bh), dr(bw), dg(bw), db(bw)
ASK PIXEL ARRAY (0,1) bc
SET BITMAP SIZE bw*4,bh

FOR y=1 TO bh
   FOR x=1 TO bw
      LET  c=bc(x,y)
      LET  rr(x,y)=MOD(c,256)
      LET  c=INT(c/256)
      LET  gg(x,y)=MOD(c,256)
      LET  bb(x,y)=INT(c/256)
   NEXT x
NEXT y

OPEN #1:NAME f$ & "_rr.csv", RECTYPE INTERNAL
ERASE #1
OPEN #2:NAME f$ & "_gg.csv", RECTYPE INTERNAL
ERASE #2
OPEN #3:NAME f$ & "_bb.csv", RECTYPE INTERNAL
ERASE #3
FOR y=1 TO bh
   FOR x=1 TO bw
      LET  dr(x)=rr(x,y)
      LET  dg(x)=gg(x,y)
      LET  db(x)=bb(x,y)
   NEXT x
   MAT WRITE #1: dr
   MAT WRITE #2: dg
   MAT WRITE #3: db
NEXT y
CLOSE #3
CLOSE #2
CLOSE #1

! 確認のためRGBそれぞれで画像を表示してみる

MAT gg=256*gg
MAT bb=65536*bb

MAT PLOT CELLS, IN 1/4,1; 2/4,0 :rr
MAT PLOT CELLS, IN 2/4,1; 3/4,0 :gg
MAT PLOT CELLS, IN 3/4,1;   1,0 :bb
END