bem_importers.f90 Source File

!> OBJメッシュを走査・解析し、内部 `mesh_type` へ変換するインポートモジュール。
module bem_importers
  use bem_kinds, only: dp, i32
  use bem_types, only: mesh_type
  use bem_mesh, only: init_mesh
  use, intrinsic :: iso_fortran_env, only: iostat_end, iostat_eor
  implicit none
contains

  !> OBJファイルを2パス(件数取得→実データ読込)で読み込み、メッシュを構築する。
  !! @param[in] path 読み込むOBJファイルのパス。
  !! @param[out] mesh OBJから構築した三角形メッシュ。
  subroutine load_obj_mesh(path, mesh)
    character(len=*), intent(in) :: path
    type(mesh_type), intent(out) :: mesh
    integer(i32) :: nvert, ntri
    real(dp), allocatable :: vertices(:, :)
    integer(i32), allocatable :: faces(:, :)

    call scan_obj(path, nvert, ntri)
    if (nvert == 0 .or. ntri == 0) error stop "OBJ has no vertices/faces"
    allocate (vertices(3, nvert), faces(3, ntri))
    call parse_obj(path, vertices, faces)
    call build_mesh_from_indexed(vertices, faces, mesh)
  end subroutine load_obj_mesh

  !> 頂点配列と三角形インデックス配列から面頂点配列(v0,v1,v2)を展開して初期化する。
  !! @param[in] vertices 頂点座標配列 `vertices(3,nvert)`。
  !! @param[in] faces 三角形の頂点インデックス配列 `faces(3,ntri)`（1始まり）。
  !! @param[out] mesh 展開後に初期化した三角形メッシュ。
  subroutine build_mesh_from_indexed(vertices, faces, mesh)
    real(dp), intent(in) :: vertices(:, :)
    integer(i32), intent(in) :: faces(:, :)
    type(mesh_type), intent(out) :: mesh
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :)
    integer(i32) :: i, ntri

    if (size(vertices, 1) /= 3 .or. size(faces, 1) /= 3) then
      error stop "vertices/faces shape mismatch"
    end if

    ntri = size(faces, 2)
    allocate (v0(3, ntri), v1(3, ntri), v2(3, ntri))
    do i = 1, ntri
      v0(:, i) = vertices(:, faces(1, i))
      v1(:, i) = vertices(:, faces(2, i))
      v2(:, i) = vertices(:, faces(3, i))
    end do
    call init_mesh(mesh, v0, v1, v2)
  end subroutine build_mesh_from_indexed

  !> OBJを行走査して頂点数と三角形分割後の面数を事前計数する。
  !! @param[in] path 走査対象のOBJファイルパス。
  !! @param[out] nvert 頂点行 `v` の総数。
  !! @param[out] ntri 面行を扇形分割した後の三角形総数。
  subroutine scan_obj(path, nvert, ntri)
    character(len=*), intent(in) :: path
    integer(i32), intent(out) :: nvert, ntri
    character(len=:), allocatable :: line
    integer :: u, ios, ntok

    nvert = 0
    ntri = 0
    open (newunit=u, file=path, status='old', action='read', iostat=ios)
    if (ios /= 0) error stop "failed to open OBJ"

    do
      call read_obj_line(u, line, ios)
      if (ios /= 0) exit
      if (is_vertex_line(line)) nvert = nvert + 1
      if (is_face_line(line)) then
        ntok = count_face_tokens(line)
        if (ntok >= 3) ntri = ntri + (ntok - 2)
      end if
    end do
    close (u)
  end subroutine scan_obj

  !> OBJの頂点/面行を解析し、負インデックス対応で配列へ格納する。
  !! @param[in] path 読み込み対象のOBJファイルパス。
  !! @param[out] vertices 解析した頂点座標配列 `vertices(3,nvert)`。
  !! @param[out] faces 扇形分割後の三角形インデックス配列 `faces(3,ntri)`。
  subroutine parse_obj(path, vertices, faces)
    character(len=*), intent(in) :: path
    real(dp), intent(out) :: vertices(:, :)
    integer(i32), intent(out) :: faces(:, :)
    character(len=:), allocatable :: line
    integer :: u, ios, i
    integer(i32) :: iv, itri
    integer(i32), allocatable :: idx(:)

    iv = 0
    itri = 0
    open (newunit=u, file=path, status='old', action='read', iostat=ios)
    if (ios /= 0) error stop "failed to open OBJ"

    do
      call read_obj_line(u, line, ios)
      if (ios /= 0) exit
      if (is_vertex_line(line)) then
        iv = iv + 1
        call parse_vertex_line(line, vertices(:, iv))
      else if (is_face_line(line)) then
        call parse_face_line(line, iv, idx)
        do i = 2, size(idx) - 1
          itri = itri + 1
          faces(:, itri) = [idx(1), idx(i), idx(i + 1)]
        end do
      end if
    end do
    close (u)
  end subroutine parse_obj

  !> OBJの1行を長さ制限なしで読み、行末の CR を除く。
  !! 固定長バッファで読むと長い面行が短い有効な面へ切り詰められるため、
  !! 非前進読込で行末まで連結する。
  !! @param[in] u 読み込み中の装置番号。
  !! @param[out] line 読み込んだ1行（行末の改行と CR を除く）。
  !! @param[out] ios 0 なら1行を読んだ。ファイル終端または読込失敗では非0。
  subroutine read_obj_line(u, line, ios)
    integer, intent(in) :: u
    character(len=:), allocatable, intent(out) :: line
    integer, intent(out) :: ios
    character(len=1024) :: chunk
    integer :: nread

    line = ''
    do
      read (u, '(a)', advance='no', size=nread, iostat=ios) chunk
      line = line//chunk(1:nread)
      if (ios == iostat_eor) then
        ios = 0
        exit
      end if
      if (ios /= 0) then
        ! 改行のない最終行は、その行を返してから次の呼出しで終端を返す。
        if (ios == iostat_end .and. len(line) > 0) ios = 0
        exit
      end if
    end do
    call strip_cr(line)
  end subroutine read_obj_line

  !> 与えられた行がOBJ頂点行(`v `)かを判定する。
  !! @param[in] line 判定対象の1行文字列。
  !! @return is_vertex_line 関数の戻り値。
  logical pure function is_vertex_line(line)
    character(len=*), intent(in) :: line
    character(:), allocatable :: adj

    adj = adjustl(line)
    is_vertex_line = (len_trim(adj) > 1) .and. (adj(1:2) == 'v ')
  end function is_vertex_line

  !> 与えられた行がOBJ面行(`f `)かを判定する。
  !! @param[in] line 判定対象の1行文字列。
  !! @return is_face_line 関数の戻り値。
  logical pure function is_face_line(line)
    character(len=*), intent(in) :: line
    character(:), allocatable :: adj

    adj = adjustl(line)
    is_face_line = len_trim(adj) > 1 .and. adj(1:2) == 'f '
  end function is_face_line

  !> 面行に含まれる頂点トークン数を数え、扇形分割時の三角形数算出に使う。
  !! @param[in] line `f` で始まるOBJ面行。
  !! @return count_face_tokens 関数の戻り値。
  integer(i32) pure function count_face_tokens(line)
    character(len=*), intent(in) :: line
    character(len=:), allocatable :: s
    integer :: pos, first, last

    s = trim(adjustl(line))
    pos = 3
    count_face_tokens = 0
    do
      call next_face_token(s, pos, first, last)
      if (first > last) exit
      count_face_tokens = count_face_tokens + 1
    end do
  end function count_face_tokens

  !> 面行の `pos` 以降にある次の空白区切りトークンの範囲を返す。
  !! 計数パスと解析パスが同じ規則でトークンを数えるための共通 tokenizer。
  !! @param[in] s `f` で始まる左詰めの面行。
  !! @param[inout] pos 走査開始位置。戻り値はトークン直後の位置。
  !! @param[out] first トークン先頭位置。トークンがなければ `first > last`。
  !! @param[out] last トークン末尾位置。
  pure subroutine next_face_token(s, pos, first, last)
    character(len=*), intent(in) :: s
    integer, intent(inout) :: pos
    integer, intent(out) :: first, last
    integer :: n

    n = len_trim(s)
    do while (pos <= n)
      if (s(pos:pos) /= ' ') exit
      pos = pos + 1
    end do
    first = pos
    last = pos - 1
    do while (pos <= n)
      if (s(pos:pos) == ' ') exit
      last = pos
      pos = pos + 1
    end do
  end subroutine next_face_token

  !> `v x y z` 形式の頂点行を3次元座標へ変換する。
  !! @param[in] line OBJ頂点行（例: `v 0.0 1.0 2.0`）。
  !! @param[out] p 解析した頂点座標 `(x,y,z)`。
  subroutine parse_vertex_line(line, p)
    character(len=*), intent(in) :: line
    real(dp), intent(out) :: p(3)
    character(len=:), allocatable :: s
    s = trim(adjustl(line))
    read (s(3:), *) p(1), p(2), p(3)
  end subroutine parse_vertex_line

  !> `f` 行の頂点参照を抽出し、`v/vt/vn` 形式から頂点インデックスのみを取り出す。
  !! @param[in] line OBJ面行（`f i j k ...`）。
  !! @param[in] nvert 負インデックス解決に使う現在有効な頂点数。
  !! @param[out] idx 抽出した頂点インデックス配列（トークン数と同じ長さ）。
  subroutine parse_face_line(line, nvert, idx)
    character(len=*), intent(in) :: line
    integer(i32), intent(in) :: nvert
    integer(i32), allocatable, intent(out) :: idx(:)
    character(len=:), allocatable :: s
    integer :: pos, first, last, slash, k
    integer(i32) :: vi

    s = trim(adjustl(line))
    allocate (idx(count_face_tokens(s)))
    pos = 3
    do k = 1, size(idx)
      call next_face_token(s, pos, first, last)
      slash = index(s(first:last), '/')
      if (slash > 0) last = first + slash - 2
      read (s(first:last), *) vi
      if (vi < 0) vi = nvert + vi + 1
      if (vi <= 0 .or. vi > nvert) error stop "OBJ face index out of range"
      idx(k) = vi
    end do
  end subroutine parse_face_line

  !> 行末の CR 文字 (char 13) をスペースに置換し、CRLF 改行の OBJ に対応する。
  !! @param[inout] line 処理対象の行文字列。
  subroutine strip_cr(line)
    character(len=*), intent(inout) :: line
    integer :: i
    do i = len_trim(line), 1, -1
      if (ichar(line(i:i)) == 13) then
        line(i:i) = ' '
      else
        exit
      end if
    end do
  end subroutine strip_cr

end module bem_importers