s-taka.org

高橋 賢

Lang:

全天球カメラを用いた天空視界の定量的評価

category: gnss
tags: image rtklib

はじめに

衛星測位や無線通信の屋外実験を行うとき、地形や建物などによる天空の遮蔽（しゃへい）を数値的に表現できると、その環境を説明しやすくなります。全天球カメラにて撮影した画像ファイルから、Pythonコードにより、空を見渡せる割合である天空視界率（SVF: sky-view factor）を求めてみました。SVFは樹冠開空度とも呼ばれます。

全天球カメラ

全天球カメラとは、すべての景色を1ショットで取り込むことのできるカメラです。これは、複数の広角レンズを備えた撮像素子により同時に撮影し、これらの映像つなぎ合わせを行って、すべての景色を1枚の写真に収めます。

つなぎ合わされた画像は、標準的なJPEG画像ファイルにて出力されますので、一般のパソコンやスマートフォンでも見ることができます。しかし、この映像は、広角レンズにて撮影されたものなので、そのままでは湾曲して表示されます。この画像ファイルは、カメラ添付アプリによりVR表示できますが、地図でお馴染みのメルカトル図法と同等の表現形式ですので、多くの方々が様々な画像形式への変換にチャレンジされています。

ここでは、広角レンズを2枚備えた、Ricoh Theta Vを用います。このカメラは、姿勢センサにより画像の天地補正（zenith correction）ができるので、また、スマートフォンアプリによる遠隔シャッタが利用できるので、とても便利です。例えば、広島市立大学RTK基準局の3Dパノラマ写真では、このカメラにて撮影した画像をJavaScriptでレンダリングしています。

天空視界率の計算

ここでは、画像輝度の明るい部分を空とみなして、画像中の空の割合を求めます。この輝度しきい値は、大津の二値化法にて求められます。

天空視界率（SVF）は、(1) 天地補正を前提に水平から天頂までの上半球の切り取り、(2) 正距円筒図法（equidistant cylindrical projection）から正積図法（equidistant and equal-area projection）への変換、(3) グレイスケール化、(4) 大津の二値化法による輝度しきい値決定、(5) 二値化の実行、(6) しきい値を超えたピクセル割合の計算、の手順にて求めます。この画像処理のために、OpenCVを、またPythonモジュールnumpy、pillow、scikit-imageを用いています。

pip install pillow scikit-image

ここで作成したPythonコードskyview.pyは次のとおりです。-fオプションで魚眼画像の出力を、-bオプションで二値化画像の出力を、-tオプションで二値化しきい値調整を、それぞれ行うようになっています。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# This file is for sky view factor calculation
# pip install pillow scikit-image

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import sys
import argparse
from scipy.ndimage import map_coordinates
from skimage.filters import threshold_otsu
from skimage.color import rgb2gray

def remap_equidistant_to_fisheye_array(img, fov_deg=180):
    # Calculate the new dimensions for the square fisheye image
    # Assuming the field of view is a full hemisphere
    fov = np.radians(fov_deg)
    h, w = img.shape[:2]
    w_new = int(w * fov / np.pi)
    h_new = w_new

    # Define the new image with the same channel count as the original
    fisheye_img = np.zeros((h_new, w_new, img.shape[2]), dtype=img.dtype)

    # Coordinates of the new image
    xx_new, yy_new = np.meshgrid(np.arange(w_new), np.arange(h_new))

    # Re-center the coordinates at the center of the fisheye image
    xx_new = xx_new - w_new / 2.0
    yy_new = yy_new - w_new / 2.0

    # Calculate the r, theta coordinates
    r = np.sqrt(xx_new**2 + yy_new**2)
    theta = np.arctan2(-yy_new, xx_new) + np.pi

    # Calculate the corresponding coordinates in the source image
    xx = (w * theta / (2 * np.pi)).astype(np.float32)
    yy = (h * r / r.max()).astype(np.float32)

    # Map the polar coordinates back to the cylindrical coordinates
    map_x = xx.clip(0, w-1)
    map_y = yy.clip(0, h-1)

    # Use scipy's map_coordinates for efficient remapping
    for i in range(img.shape[2]):
        fisheye_img[:, :, i] = map_coordinates(img[:, :, i], [map_y.flatten(), map_x.flatten()], order=1).reshape((h_new, w_new))

    return Image.fromarray(fisheye_img)


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Calculate the sky view factor of a fisheye image')
    parser.add_argument('images', type=str, nargs='+', help='Path to the original image')
    parser.add_argument('-f', '--fisheye', action='store_true' , help='switch to save fisheye image')
    parser.add_argument('-b', '--binary', action='store_true', help='switch to save binary image')
    parser.add_argument('-t', '--threshold', type=float, default=0.0, help='threshold offset')
    args = parser.parse_args()
    for image in args.images:  # For each image in the command line arguments
        # Path to the original image
        img_original = np.array(Image.open(image))

        # Cut the original image to only keep the sky
        img_sky = img_original[:img_original.shape[0] // 2, :]

        # Convert the sky part to fisheye
        fisheye_sky_image = remap_equidistant_to_fisheye_array(img_sky)

        # Save the fisheye sky image
        if args.fisheye:
            path_image_fisheye = image.split('.')[0] + '_fish.jpg'
            fisheye_sky_image.save(path_image_fisheye)

        # Convert the image to grayscale for thresholding
        gray_fisheye_sky_image = rgb2gray(fisheye_sky_image)

        # Find the Otsu threshold value
        threshold_value = threshold_otsu(gray_fisheye_sky_image)
        threshold_value += args.threshold

        # Apply the threshold to binarize the image
        binary_fisheye_sky_image = gray_fisheye_sky_image > threshold_value

        # Calculate the proportion of sky and clouds
        # Assuming that sky is white after thresholding
        sky_ratio = np.sum(binary_fisheye_sky_image) / binary_fisheye_sky_image.size

        if args.binary:
            path_image_sky = image.split('.')[0] + '_bin.jpg'
            # Convert the binary image back to uint8 for saving
            binary_fisheye_sky_image = (binary_fisheye_sky_image * 255).astype(np.uint8)
            binary_fisheye_sky_image = Image.fromarray(binary_fisheye_sky_image)

            # Save the binary image
            binary_fisheye_sky_image.save(path_image_sky)

        print(f'{image} sky ratio: {sky_ratio * 100:.1f}%')

全天球カメラで撮影した次の写真を例に、SVFを求めます。

240318-144315.jpg Original image of an example of sky-view factor calculation

多くの画像を整理するため、私は、カメラ内部時刻を設定の上、jheadにてJPEGファイルのEXIF（イグジフ、Exchangeable image file format）情報をもとにファイル名を日時を含むものに変更しています（jhead -nf%y%m%d-%H%M%S -autorot *.jpg）。私は、仕事もプライベートも、この方法で画像を整理しています。Macならばbrew install jheadで、Debian Linuxならばapt install jhead -yで、jheadを簡単にインストールできます。

この画像ファイル240318-144315.jpgを処理します。これは、2024-03-18 14:43:15 JSTに撮影したものです。

$ ./skyview.py -f -b 240318-144315.jpg
240318-144315.jpg sky ratio: 52.4%

SVFは52.4%と評価されました。ここで出力された魚眼画像240318-144315_fish.jpgと二値化画像240318-144315_bin.jpgは次のとおりです。

240318-144315_fish.jpg fish-eye image for sky-view factor calculation 240318-144315_bin.jpg binary image for sky-view factor calculation

二値化画像において、明らかな空の部分が遮蔽されているとして扱われています。そこで、空の部分を空と認識されるように、しきい値の調整と2値化画像の表示を繰り返しました。最終的に、この画像の二値化しきい値として-0.30を選択し、このSVFを92.2%と評価しました。

$ ./skyview.py -b -t -0.30 240318-144315.jpg
240318-144315.jpg sky ratio: 92.2%

240318-144315_bin.jpg adjusted binary image for sky-view factor calculation

このしきい値では、建物の一部が空とみなされてしまいました。しきい値調整のみの方法では、この程度が限界だと思っています。

天空視界率の観測

全天球カメラとスマートフォンを外に持ち出し、職場の周囲道路の天空視界率を観測します。あらかじめ、自動車に乗せたGPS受信機により大学周辺道路を測位しながら走行し、RTKLIBのrtkplot.exeにより走行コースをプロットして、印刷しておきました。

全天球カメラを走行コース上に設置し、スマートフォンアプリにて遠隔撮影を行って、その時刻を走行コース上に記録します。アナログな方法です。

experiment trajectory

その後、全天球カメラから画像ファイルをPCにダウンロードし、上述のjheadを用いた方法にてファイル名に日時を含むようリネームしました。そして、次のようなシェルスクリプトを作成し、すべての画像ファイルに対して上述のskyview.pyを実行し、魚眼画像と二値画像を見比べ、しきい値を調整しました。かなりアナログな方法です。

./skyview.py -b -f -t -0.30 240318-144315.jpg
./skyview.py -b -f -t -0.30 240318-144447.jpg
./skyview.py -b -f -t -0.30 240318-144804.jpg
./skyview.py -b -f -t -0.30 240318-144957.jpg
./skyview.py -b -f -t -0.00 240318-145410.jpg
./skyview.py -b -f -t -0.20 240318-145529.jpg
./skyview.py -b -f -t -0.10 240318-145731.jpg
./skyview.py -b -f -t -0.10 240318-150032.jpg
./skyview.py -b -f -t  0.00 240318-150231.jpg
./skyview.py -b -f -t -0.25 240318-150350.jpg
./skyview.py -b -f -t -0.25 240318-150716.jpg
./skyview.py -b -f -t  0.00 240318-151024.jpg
./skyview.py -b -f -t -0.15 240318-151251.jpg
./skyview.py -b -f -t  0.13 240318-151456.jpg
./skyview.py -b -f -t -0.30 240318-151751.jpg

obtaining sky-view factor thresholds

このようにして、この走行コース上のSVFを求めました。

240318-144315.jpg sky ratio: 92.2%
240318-144447.jpg sky ratio: 88.3%
240318-144804.jpg sky ratio: 91.4%
240318-144957.jpg sky ratio: 94.5%
240318-145410.jpg sky ratio: 64.9%
240318-145529.jpg sky ratio: 76.0%
240318-145731.jpg sky ratio: 74.7%
240318-150032.jpg sky ratio: 72.8%
240318-150231.jpg sky ratio: 63.9%
240318-150350.jpg sky ratio: 88.4%
240318-150716.jpg sky ratio: 85.8%
240318-151024.jpg sky ratio: 53.2%
240318-151251.jpg sky ratio: 78.7%
240318-151456.jpg sky ratio: 69.3%
240318-151751.jpg sky ratio: 85.0%

そして、ファイル名の時刻情報をもとに、SVF値を走行コースに書き込みました。

これらのJPG画像ファイルのEXIF情報にスマートフォン座標が自動記録されることを期待しましたが、そのような記録はありませんでした。残念。

$ jhead 240318-144315.jpg
File name    : 240318-144315.jpg
File size    : 4950603 bytes
File date    : 2024:03:18 21:07:33
Camera make  : RICOH
Camera model : RICOH THETA V
Date/Time    : 2024:03:18 14:43:15
Resolution   : 5376 x 2688
Flash used   : No
Focal length :  1.3mm
Exposure time: 0.0008 s  (1/1250)
Aperture     : f/2.0
ISO equiv.   : 64
Whitebalance : Auto
Metering Mode: pattern
Exposure     : program (auto)
JPEG Quality : 95
======= IPTC data: =======
City          :
Record vers.  : 2
Time Created  : 144315+0900
Caption       :
DateCreated   : 20240318

今後は、自動車のルーフに全天球カメラとGPS受信機を設置し、その動画から走行コースとSVFを自動計算できるようにしたいと思っています。また、Unicore UM982受信機などの進行方向（ヘディング）情報の取れるGPS受信機と画像処理により、全天球映像の正面が常に真北になるように画像記録したいです。現在の方法は、とてもアナログな方法なので。

まとめ

全天球カメラの映像から天空視界率（SVF）を求めました。大津の二値化法を用いて空の割合を求めてきましたが、画像を見ながらのしきい値の調整は必要になりました。撮影時刻情報をもとに撮影場所を同定して、撮影場所とSVFを対応づけました。手作業の部分が多いので、今後は自動化を考えてみたいと思っています。