数値解析法及び演習 第十二回

画像解析
  1. 画像の表現
  2. 4バイトデータの入出力
  3. 1バイトデータの入出力
  4. 複数ファイルのオープン
  5. 全地球の温度上昇トレンド
  6. 今日の宿題

[1] 画像の表現


[2] 4バイトデータの入出力

今回用いる白黒画像は多数の1バイトデータからなります.1バイトのデータ(8個の0と1)をそのまま書き込んだファイルをバイナリファイルと呼びます.2の8乗=256個の整数を表現できます.一方で今まで扱ってきた,人間が読めるファイルは「テキストファイル」といって区別されます.シフトJIS(macのデフォルト文字コード)では,半角文字で1バイト,全角文字で2バイトのメモリを使います.256は3桁ですから,半角の数字で表現した場合3バイト必要になります.したがって,テキストでファイルに書き出した場合とバイナリで書き出した場合とで,必要な記憶容量が3倍も異なることになります.桁数が多くなればもっと差が大きくなります.

今回あつかうデータはバイナリファイルのため,バイナリファイルの読み書きを行う必要があります.次のプログラムは,0.1から1.0までの10個の単精度実数を生成し,それをバイナリ形式でファイル10.datに書き込むプログラムです.

      implicit none
      integer i
      real x
      open(1,file='10.dat',access='direct',form='unformatted',recl=4)

      do i=1,10
         x=0.1*i
         write(1,rec=i) x
      enddo

      end

ここで,form='unformatted'でバイナリ形式を指定し,recl=4で一つのデータ(レコード)の大きさ(今の場合単精度実数=4バイト=32ビット)を指定します.access='direct'は,直接参照ファイルであることを指定します.ファイル中の場所を指定することで,特定の数値に直接アクセスします.今までは'sequential'となっていて,ファイルの先頭から順番にアクセスすることが暗に指定されていました.またwrite文では,rec=iでレコード番号を指定します.このプログラムを実行し,生成された10.datのサイズが40バイトになっていることを確認しましょう.

先ほど生成したデータを読み込んで画面に表示するプログラムは以下になります.

      implicit none
      integer i
      real x
      open(1,file='10.dat',access='direct',form='unformatted',recl=4)

      do i=1,10
       	 read(1,rec=i) x
	 write(*,*) x
      enddo

      end

実行して結果を確認しましょう.


[3] 1バイトデータの入出力

今回扱う画像は白黒256階調なので,各ピクセルは1バイトのデータで表現されます.FORTRANでは符号なし1バイトデータを扱う適当な形式が用意されていないので,character文で宣言される文字型変数に数値を格納します.

それぞれの文字には文字コードと呼ばれる数字が対応しています.そこで,0から255までの整数を文字に一度変換し,その文字をファイルに書き込むと,0から255までの整数(本来は文字コード)がファイルに書き出されます.

0から255までの整数をバイナリデータとしてファイル256.datに書き込み,そのデータを読み出して0から255までの整数を画面に表示するプログラムは以下になります.

      implicit none
      character n
      integer i

      open(1,file='256.dat',access='direct',form='unformatted',recl=1)

      do i=0,255
         n=char(i)
         write(1,rec=i+1) n
      enddo

      do i=1,256
         read(1,rec=i) n
         write(*,*) i,ichar(n)
      enddo

      end

ポイントは,256階調の整数はchar()で文字型変数 n に変換してファイルに書き込む.読み出す際は文字型変数 n にまずは代入し,それをichar()関数で256の整数に変換する です.

練習問題1

このデータは,1980年-2000年における全地球の表面気温のデータです.360x180=64800個の1バイトデータ(0-255)が並んだファイルが21個あります.0が-15度,255が25度に対応します.それぞれのデータは地球上のある経度(0-359),緯度(0-179)における温度です.注意してほしいのは,二次元的にデータが並んでいるわけではなく,一次元的にデータは並んでいることです.データの頭は北極(緯度0)であり,まずは経度方向に360個のデータが並びます.361番目のデータは緯度1,経度0となります.360個のデータの組みが180個並んでいる,という構造になってます.

ダウンロードしたら

  mv ~/Downloads/data/* ./
とすることで現在のディレクトリの中にデータファイルが移動します.lsで確認してください.

display -depth 8 -size 360x180 gray:conv1980.raw を実行すると,1980年における温度分布が画像として表示されます.黒が低温,白が高温です.displayコマンドが使えない場合はこちら

次に,以下のプログラムを実行し,数値(0-255)が画面に表示されることを確認してください.

      implicit none
      integer i
      character n

      open(1,file='conv1980.raw',access='direct',form='unformatted',recl=1)

      do i=1,64800
         read(1,rec=i) n
         write(*,*) i,ichar(n)
      enddo

      end

    

[3] 複数ファイルのオープン

今回は,1980年から2000年までの21年分のデータを解析し,その温度上昇のトレンドが地域ごとにどう分布しているかを調べます.

そのためには,ファイルを21個開かないといけません.正直にプログラムに記述してもいいのですが効率が悪いので,文字処理のテクニックを用いてファイルを操作します.次のプログラムを実行すると,21個のファイルを順番にオープンしてそのファイル名が画面に表示されます.

      implicit none
      integer i,num,nyear,j
      character name*12,year*4,n

      nyear=21

      do i=1,nyear

         num=1979+i
         write(year,'(I4)') num  ! 整数numを文字型に変換しyearに代入
         name='conv'//year//'.raw' ! conv, year, .rawからファイル名を生成


         open(10+i,file=name,access='direct',form='unformatted',recl=1)
         ! ここでnameには上で作られた文字列が代入される
      
         do j = 1, 64800
           read(10+i,rec=j) n
           write(*,*) num,j,ichar(n)
         enddo

      enddo

      end


write(year,'(I4)') num で,4桁の整数numが文字となって文字型変数yearに格納されます.


[4] 全地球の温度上昇トレンド

さてここまでくると,全地球の温度変化を解析できるようになりました.以下の手順を踏んで,地球上のどこで温暖化がより進行しているのか解析しましょう.

練習問題2

  1. 21個のファイルからそれぞれ64800個のデータを読み込み,0-255の整数から温度に変換し,結果を配列temp(21,64800)に格納する.
  2. 64800個すべての地点において最小二乗法を用いることにより温度変化の傾き(℃/年)を算出する.
  3. 2で算出した数値を256階調の整数に変換する.0から0.1(℃/年)までを,0から255までの整数に対応させる.0より小さい値には0,0.1よりも大きい値には255を対応させる.こうすることで,温度がより上昇している地域は白っぽく見えることになる.
  4. 64800個の0から255までの整数をcharacter変数を用いることでバイナリファイルplus.rawに書き込む.
  5. 同様に,0から-0.1(℃/年)までを,0から255までの整数に対応させる.0より大きい値には0,-0.1よりも小さい値には255を対応させる.こうすることで,温度がより下降している地域が白っぽく見えることになる.こちらの結果は,minus.rawに書き込む.
  6. プログラムの大まかな構造は以下のようになるでしょう.
    
        implicit none
    !    宣言文など
        
        
        open(1,file='plus.raw',access='direct',form='unformatted',recl=1) ! 書き込むファイルのオープン
    
        do i=1,nyear
          
             num=1979+i
             write(year,'(I4)') num
             name='conv'//year//'.raw'
             write(*,*) name
          
             open(10+i,file=name,access='direct',form='unformatted',recl=1) !各年のファイルのオープン
    
          do j = 1, 64800
    
        
     ! 10+i番のファイルから数値(0-255)を64800個読み込み,温度に変換し配列T(i,j)に格納する.読み出し方は
              read(10+i,rec=j) n
        T(i,j) = ichar(n)  ! 配列T(i,j)に温度データ(-30度-30度)が代入される
        
          enddo
          
        enddo
    
    
        do j = 1, 64800
           do i = 1, nyear
    
        
    ! この2重ループで,最小二乗法により64800個の傾きを算出する
        ! xが年(1980-2000) yが温度T(年i,場所の番号j)となっている.
        
           enddo
        enddo
    
        do j = 1, 64800 ! 64800個の傾きの書き出し
    
        !    ここで温度上昇の傾きから0-255の整数に変換する
        !    プラスの場合とマイナスの場合と二通り実行する必要あり
    
        n = char(変換後の整数がここに入る)
        write(1,rec=i) n ! ファイルへの書き込み
        
       enddo
    
      
  7. display -depth 8 -size 360x180 gray:plus.raw (およびminus.raw) を実行する.

画像がだせたら,以下の点を観察してください.

練習問題3

  1. 温度が上昇(下降)している地域はどこか?
  2. 陸域と海洋でどちらがムラが多いか?それはなぜか?

[5] 今日の宿題

プログラムおよび練習問題3の答えを送ること. 〆切は1/18.


日程表へ戻る <<