anime gif 画像出力ルーチン OO 化

type bound procedure で OO 風

少しバグがあって ビット列長が 255byte の倍数になっている時、余分な 00 が出力されていました。

それが丁度、intel compiler の associate の関わる bug と思われるものと重なってちょっと苦しみました。associate のポインタがらみのバグも中々取れてません。

図画用の変数はスカラーの allocatable として宣言して、コンストラクタで割り付けと初期化を行います。fig = t_gif('mandel.gif', nx, ny, n_global_color, global_color_table)

FINALIZE を使って、ファイルのクローズやメモリ解放処理を行おうとしたのですが、一方で関数でコンストラクタを作ったために、代入・自動割り付け後のコンストラクタの作った右辺のオブジェクトが解放される時にファイナライズ処理が一度なされてしまい、ファイルがクローズされるので、FINALIZE はやめて、明示的に終了処理をすることにしました。気長に解決法を考えたいと思います。

(H30.5.9)
二次元配列を一次元にするのに、Fortran2008 で拡張された pointer assignment を用いました。(MFE §15.6) これによってテンポラリ配列への代入の必要がなくなりました。

(H30.5.11)
２色の時の扱いが間違っていたので修正。

出力結果

アニメでサイケ調のマンデルブロ集合の収束を眺める

２色

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Scientific Software Design: The Object-Oriented Way

ソース・プログラム

gfortran には findloc が無いようなので、代替品をば。あと forall(integer:: ...) にも対応していないようです。

        pure integer function my_findloc(dict, k) result(ires)
          integer, intent(in) :: dict(:), k
          do ires = 1, size(dict)
            if (dict(ires) == k) exit
          end do  
          if (ires > size(dict)) ires = 0 ! not found
        end function my_findloc

本体

    module m_gif_types
      use, intrinsic :: iso_fortran_env
      implicit none    
      private
      public :: t_gif_header, t_image_block, t_graphic_control_extension, t_application_extension 
      ! constructors 
      interface t_gif_header
        procedure :: init_gif_header
      end interface 

      interface t_image_block
        procedure :: init_image_block
      end interface 
      
      interface t_graphic_control_extension
        procedure :: init_graphic_control_extension 
      end interface
      
      interface t_application_extension
        procedure :: init_application_extension
      end interface  
      !
      type :: t_gif_header
        sequence  
        character(3)   :: signature               = 'GIF'
        character(3)   :: version                 = '89a' ! '87a'
        integer(int16) :: width     
        integer(int16) :: height    
        integer(int8)  :: pck                     = int(B'10001010', int8) !1:Global Col. 3:Color Res. 1:Sort Flg. 3:Size of Global Col.
        integer(int8)  :: background_color_index  = 0
        integer(int8)  :: pixel_aspect_ratio      = 0
      end type  

      type :: t_image_block
        sequence
        integer(int8)  :: image_separator         = int(Z'2C', int8)
        integer(int16) :: image_left_position     = 00
        integer(int16) :: image_top_position      = 00
        integer(int16) :: image_width  
        integer(int16) :: image_height 
        integer(int8)  :: pck                     = int(B'00000000', int8)
        ! local color table
        !integer(int8) :: LZW_minimum_code_size 
        !integer(int8) :: block_size  
        !integer(int8) :: image_data (:)
        !integer(int8) :: block_terminator        = int(Z'00', int8)  
      end type t_image_block

      type :: t_graphic_control_extension
        sequence  
        integer(int8)  :: extention_introducer    = int(Z'21', int8)
        integer(int8)  :: graphic_control_label   = int(Z'F9', int8)
        integer(int8)  :: lock_size               = int(Z'04', int8)
        integer(int8)  :: pck                     = int(Z'00', int8)
        integer(int16) :: delay_time              = 50 ! msec
        integer(int8)  :: transparent_color_index = 0
        integer(int8)  :: block_terminator        = int(Z'00', int8)
      end type t_graphic_control_extension
      
      type :: t_application_extension ! for animation gif
        sequence  
        integer(int8)  :: extension_intrducer     = int(Z'21', int8)
        integer(int8)  :: extension_label         = int(Z'FF', int8)
        integer(int8)  :: block_size_01           = int(Z'0B', int8)
        character(len = 8) :: application_identifier          = 'NETSCAPE'
        character(len = 3) :: application_authentication_code = '2.0'
        integer(int8)  :: block_size_02           = int(Z'03', int8) ! 0:block terminator
        integer(int8)  :: n                       = int(Z'01', int8) !application_data
        integer(int16) :: nloop                   = 0                ! 0:unlimited loop 
        integer(int8)  :: block_terminator        = int(Z'00', int8)
      end type t_application_extension   
    contains
      ! constructors 
      pure type(t_gif_header) function init_gif_header(nx, ny, ngcol) result(res)
        integer, intent(in) :: nx, ny, ngcol
        res%width  = nx
        res%height = ny
        res%pck = iand(res%pck, int(B'11111000', int8)) + int(mod(ngcol, 8), int8) !least 3bits:size of global color:2**(ngcol + 1) 2..256
      end function init_gif_header
    
      pure type(t_image_block) function init_image_block(nx, ny) result(res)
        integer, intent(in) :: nx, ny
        res%image_width  = nx
        res%image_height = ny
      end function init_image_block
    
      pure type(t_graphic_control_extension) function init_graphic_control_extension(it) result(res)
        integer, intent(in) :: it
        res%delay_time = int(it, int16) ! msec 
      end function  init_graphic_control_extension
    
      pure type(t_application_extension) function init_application_extension(nloop) result(res)
        integer, intent(in) :: nloop
        res%nloop = int(nloop, int16)
      end function init_application_extension  
    end module m_gif_types 

    module m_gif_lzw
      use, intrinsic :: iso_fortran_env
      use :: m_gif_types
      implicit none
      private
      public :: encoder
      type :: t_enc
        integer :: kbits = 0, id = 0
      end type    
    contains
      subroutine enc(irgb, nbits, code)               ! gif LZW 
        integer, intent(in) :: irgb(:), nbits
        type (t_enc), intent(out) :: code(:)
        integer :: dict(0:2**12 - 1) ! 0:4095
        integer :: kbits, kd, m0                      ! <== m0 state variable
        integer :: ip, ic, ienc
        call clear_dict()                             ! clear dictionary 
        ic = 1
        code(ic) = t_enc(kbits, 2**nbits)             ! clear_code 
        m0 = irgb(1)
        do ip = 2, size(irgb)
          call subenc(ienc, irgb(ip)) 
          if (ienc /= -1) then                        ! new word 
            ic = ic + 1
            code(ic) = t_enc(kbits, ienc)
          end if  
          if (kd == 2**12 - 1) then                   ! dictionary full
            ic = ic + 1  
            code(ic) = t_enc(kbits, 2**nbits)         ! clear_code
            call clear_dict()                         ! clear dictionary 
          end if
        end do
        ic = ic + 1
        code(ic) = t_enc(kbits, m0)
        ic = ic + 1 
        code(ic) = t_enc(kbits, dict(2**nbits + 1))   ! end_code  
      contains
        subroutine subenc(ienc, m1)
          integer, intent(out) :: ienc
          integer, intent(in ) :: m1
          integer :: k, id
          k = ishft(m0 + 1, 16) + m1                  ! m0 + 1 to avoid 00...0 degeneracy 
          id = findloc(dict, k, dim = 1)              ! dictionary is 0-based 
          if (id == 0) then                           ! not found in the dictionary
            kd = kd + 1
            dict(kd) = k
            if (kd == 2**kbits + 1) kbits = kbits + 1
            ienc = m0
            m0 = m1
          else                                        ! found in the dictionary
            ienc = -1
            m0 = id - 1                               ! because dictionary is 0-based, m0 must be shifted by 1
          end if    
        end subroutine subenc

        subroutine clear_dict()
          kbits = nbits + 1
          kd = 2**nbits + 1 ! dictinary defined
          dict = 0
          forall(integer :: i = 0:kd) dict(i) = i    ! dict(nbit):clear code; dict(nbit+1):end code 
        end subroutine clear_dict
      end subroutine enc
      
      subroutine encoder(irgb, nbits, icode)
        integer, intent(in ) :: irgb(:), nbits
        integer(int8), allocatable, intent(out) :: icode(:)
        type(t_enc), allocatable :: code(:)
        integer :: i, j, k, ib, ic, nb  
        allocate(code(size(irgb) + 100))              ! ? large enough ?
        call enc(irgb, nbits, code)                   ! compress color code to LZW code
        nb = ceiling(sum(code(:)%kbits) / 8.0)        !   required bytes  
        allocate(icode(nb))
        icode = 0                                     ! pack LZW code to bit string    
        k = 0
        do i = 1, size(code)                          
          do j = 1, code(i)%kbits   
            if (btest(code(i)%id, j - 1)) then 
              ic = k / 8 + 1 
              ib = mod(k, 8)
              icode(ic) = ibset(icode(ic), ib)
            end if
            k = k + 1
          end do  
        end do            
      end subroutine encoder 
    end module m_gif_lzw
        
    module m_anime_gif
      use, intrinsic :: iso_fortran_env
      use :: m_gif_types
      use :: m_gif_lzw
      implicit none
      private 
      public :: t_gif
      interface t_gif 
        procedure :: init_gif
      end interface
      type :: t_gif
        integer :: iw, nbits
        integer, allocatable :: irgb2(:, :)
      contains
        procedure :: wr_gif_img, gif_close
      end type t_gif
    contains
      impure type(t_gif) function init_gif(fn, nx, ny, n_global_color, global_color_table, loop) result(this)
        character(*), intent(in) :: fn
        integer, intent(in) :: nx, ny, n_global_color, global_color_table(:)
        integer, intent(in), optional :: loop
        type(t_gif_header) :: head
        type(t_application_extension) :: anime
        integer :: nloop = 0
        if (present(loop)) nloop = loop
        associate(iw => this%iw, nbits => this%nbits)
          nbits = max(2, n_global_color + 1) ! 1..8 -> 2,2,3,4,5,6,7,8
          allocate(this%irgb2(nx, ny))
          head  = t_gif_header(nx, ny, n_global_color)
          anime = t_application_extension(nloop)
          open(newunit = iw, file = fn, access = 'stream') 
          write(iw) head      
          write(iw) int(global_color_table, int8)
          write(iw) anime
        end associate  
      end function init_gif

      subroutine wr_gif_img(this, itime)
        class(t_gif), target, intent(in) :: this
        integer, intent(in) :: itime
        type(t_graphic_control_extension), allocatable :: ext
        type(t_image_block), allocatable :: img
        integer(int8), allocatable :: icode(:)
        integer :: i, nbyte  
        integer, pointer :: irgb(:)
        associate(iw => this%iw, nbits => this%nbits, irgb2 => this%irgb2)
          ext = t_graphic_control_extension(itime)
          img = t_image_block(size(irgb2, 1), size(irgb2, 2))
          irgb(1:size(irgb2)) => irgb2      ! 1D => 2D 
          call encoder(irgb, nbits, icode)  ! color code to GIF LZW code
          write(iw) ext                     ! write bit strings to file
          write(iw) img     
          write(iw) int(nbits, int8)        ! LZW_minimum_code_size = 2,2,3,4,5,6,7,8bits 
          do i = 1, size(icode), 255         
            nbyte = min(i + 255, size(icode)) - i
            write(iw) achar(nbyte)             
            write(iw) icode(i:i+nbyte-1)
          end do  
          write(iw) achar(00)               ! block_terminator
        end associate
      end subroutine wr_gif_img
      
      subroutine gif_close(this)
        class(t_gif), intent(in out) :: this
        associate(iw => this%iw)
          write(iw) achar(int(Z'3B'))                       
          close(iw)
        end associate  
      end subroutine gif_close
    end module m_anime_gif
    
    program mandel
      use :: m_anime_gif
      implicit none
      integer, parameter :: nx = 500, ny = 500, maxiter = 254
      real   , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5
      complex, allocatable :: c(:, :), z(:, :)                 
      type(t_gif), allocatable :: fig
      allocate(c(nx, ny), z(nx, ny))                     ! make 2D-mesh: c(:, :)   
      forall(integer :: ix = 1:nx, iy = 1:ny) &
          c(ix, iy) = cmplx((x1 - x0) / (nx - 1) * (ix - 1), (y1 - y0) / (ny - 1) * (iy - 1)) + cmplx(x0, y0)
      block                                              ! set global color table
        integer :: n_global_color                        ! 0..7 ! 2..256 colors (2^(n_global_color + 1)) 
        integer, allocatable :: global_color_table(:)
        n_global_color = 7
        allocate(global_color_table(3 * 2**(n_global_color + 1)))
        global_color_table = 0
        forall(integer :: i = 0:254) global_color_table(3*i+mod(i, 3)+1) = 255 ! global_color_table(3*i+1:3*i+3) = i 
        fig = t_gif('mandel.gif', nx, ny, n_global_color, global_color_table)
      end block
      block                                              !  Mandelbrot Set
        integer :: iter
        fig%irgb2 = 0
        do iter = 0, maxiter 
          where (abs(z) <= 2.0) 
            z = z**2 + c
            fig%irgb2 = fig%irgb2 + 1
          end where 
          if (iter < 10 .or. mod(iter, 10) == 1) call fig%wr_gif_img(max(80 - iter, 10))
        end do
        call fig%wr_gif_img(500)
      end block
      call fig%gif_close()
      deallocate(fig)
    end program mandel

（H30.5.9/10 微改善, 11 コメント修正)