# er2ed.R
#
# RICOH THETA で得られる正距円筒図法(Equirectangular projection)の
# 写真画像のうち地平線上の範囲を、等距離射影方式(Equidistance 
# projection)の画像に変換するRスクリプト。
#
#  2021-08-28  TAKENAKA, Akio
#
# 
#  使用法
#
#  本スクリプトはjpeg パッケージを読み込んで使用する。
#  インストールされていなければ、Rguiのメニューバーの「パッケージ」
#  からインストールしておく。
# 
#  変換したい画像ファイル(jpeg形式）をひとつのディレクトリにまとめて入れる。
#  Rgui を起動し、上記のディレクトリを作業ディレクトリに設定する。
# 
#  source('equirec2equidist.R') などとして本スクリプトを実行すると、
#  作業ディレクトリ内のjpeg ファイルを自動的に読み込んで変換する。
#  ただしjpegファイルでも縦横比が 1:2 ではないファイルはTHETA が作成したものでは
#  ないと判断してスキップする。
#  作業ディレクトリ内に eqdという名前のサブディレクトリが作成され、
#  この中に  変換後の画像ファイルが保存される。 作成されるファイルの名前は、
#  元のファイル名の最初に"eqd_" を付加したもの（これも変更可）。
#
#  作成される等距離射影の画像の大きさは、スクリプト中の radius で設定している
#  (初期設定は 500ピクセル）。 radius に設定する値を変えれば、大きさを変更
#  することができる。
#  また、地平線に対応する円周の外側に余白ができる。この大きさは、左右はmargin.x、
#  上下はmargin.y で設定する（初期設定はそれぞれ100ピクセルと50ピクセル）。
#  
#  作成される等距離射影の画像の周囲には灰色の点線が描かれている。
#  これが不要なら、draw.horizon をFALSEに設定する（初期設定はTRUE)。
#
#  上記で変更可となっている設定は、スクリプト中の
#  「### ↓諸設定。いずれも編集可 ↓ ##」
#  と書かれた行のあとにある。
#
#  参考文献：
#  Sky view factor measurement by using a spherical camera
#  Tsuyoshi HONJO, Tzu-Ping LIN, Yuhwan SEO (2019)
#  Journal of Agricultural Meteorology  Vol. 75 Issue 2 Pages 59-66
#   DOI https://doi.org/10.2480/agrmet.D-18-00027
#
# -------------------------------
# 著作権表示
# 
# Copyright (c) 2021 by TAKENAKA, Akio
#
# This software is released under the MIT License.
# http://opensource.org/licenses/mit-license.php
#
# このスクリプトは MIT ライセンスに従って提供されています。
# 
# 以下に定める条件に従い、本ソフトウェアおよび関連文書のファイル（以下「ソフトウェア」）の
# 複製を取得するすべての人に対し、ソフトウェアを無制限に扱うことを無償で許可します。これには、
# ソフトウェアの複製を使用、複写、変更、結合、掲載、頒布、サブライセンス、および/または販売
# する権利、およびソフトウェアを提供する相手に同じことを許可する権利も無制限に含まれます。
#
# 上記の著作権表示および本許諾表示を、ソフトウェアのすべての複製または重要な部分に記載する
# ものとします。
#
# ソフトウェアは「現状のまま」で、明示であるか暗黙であるかを問わず、何らの保証もなく提供され
# ます。ここでいう保証とは、商品性、特定の目的への適合性、および権利非侵害についての保証も
# 含みますが、それに限定されるものではありません。 作者または著作権者は、契約行為、不法行為、
# またはそれ以外であろうと、ソフトウェアに起因または関連し、あるいはソフトウェアの使用または
# その他の扱いによって生じる一切の請求、損害、その他の義務について何らの責任も負わないものと
# します。
# -------------------------------

# jpeg パッケージを読み込む。あらかじめインストールしておくこと。
library(jpeg)  

### ↓諸設定。いずれも編集可 ↓ ##

# 作成する等距離射影画像の半径(画素数）
radius = 500

# 作成する画像で、投影された円の周囲に設けるマージンのサイズ（画素数）
margin.x <- 100
margin.y <- 50

# 作成した変換画像ファイルを収めるディレクトリ。
# 作業ディレクトリの直下に作成される。
subdir <- 'eqd'

# 作成する等距離射影の画像ファイルのファイル名の最初に付加する文字列。
# 空文字列を設定すれば、元ファイルと同名で作成する。
eqd.prefix <- 'eqd_'

# 地平線を点線で描く。不要なら FALSE を設定
draw.horizon <- TRUE

###  ↑ここまで初設定↑  ###


# 作成したを用意する。
# すでに存在するならなにもしない（オブションで警告の表示を抑制）。
dir.create(subdir, showWarnings = FALSE)

# 作業ディレクトリ直下のファイルのリスト
files.here <- list.files('./')  

# 拡張子が .jpg ないしは jpeg のものを取り出す（大文字でも可）
files.jpg <- files.here[grep("*.(jpg|jpeg)$", files.here, ignore.case = TRUE)]


for (src.file in files.jpg) { # 各jpg ファイルを処理

    # 画像の読み込み
    # 正常に読み込めれば配列、失敗すればtry-error クラスのオブジェクトが返される。

    source.img <- try(readJPEG(src.file))  

    if (class(source.img)  == "try-error") {
        next  # 読み込みに失敗したファイルはスキップする
    }

    pxl.y <- nrow(source.img)  # 元画像の縦（天頂角）方向の画素数（180度相当）
    pxl.x <- ncol(source.img)  # 元画像の横（方位角）方向の画素数（360度相当）

    # 縦横の比率が 1:2 でなければこの画像はスキップする
    # 正距円筒射影の画像なら 1: 2 (180度:360度)のはず
    if (pxl.y * 2 != pxl.x) {
        cat(src.file, "skipped.  !! H:W is not 1:2 !! \n")
        next
    }

    # 作成する等距離射影画像の情報を持たせる 3次元配列(x, y と RGB)
    equidistant.img <-array(0, dim = c(2 * (radius + margin.y), 2 * (radius + margin.x), 3))

    # 元画像の上辺は天頂角 0度、下端は180度。90度までが解析対象。
    # 元画像の左辺は方位角 0度、右に向かって360度まで方位角が大きくなる。
    #
    # 作成する画像の画素ごとに、対応する元画像の画素を探してコピーする
    # 等距離射影の画像の中心を (0, 0) とする座標で計算。

    for (x in -radius : radius) {  # 横軸方向、(中心-半径)から(中心+半径)まで
        for (y in -radius : radius) {  # 縦軸方向、(中心-半径)から(中心+半径)まで

            # 天頂からの距離
            d.zenith <- sqrt(x ^2 + y ^ 2)
            
            # 地平線以下の点ならスキップ
            if (d.zenith > radius) {
                next
            }
            
            # その点の天頂角を計算
            zenith.angle <- pi / 2 * (d.zenith / radius)

            # その点の方位角を計算
            # atan の値域を象限（xの正負、yの正負の組み合わせ）により調整する必要がある。
            # 正距円筒射影の画像で、北（方位角0度）→東→南（方位角180度）と左から右に視線を
            # 動かすことは、 空を見上げた等距離射影の画像では反時計回りとなること、
            # またコンピュータ画面上ではx軸（横軸）のプラス方向は左から右と通常通りだが、
            # y軸（縦軸）のプラス方向は上から下となることに注意。

            if (y != 0) { # atan(x/y) で分母が 0 でない場合
                azimuth <- atan(x / y)  # x < 0, y < 0
                if (y > 0) {
                    azimuth <- azimuth + pi
                }
                else if (x > 0) {
                    azimuth <- azimuth + 2 * pi
                }
            }
            else { # atan(x/y) で分母が 0 の場合の例外処理
                if (x > 0) {
                    azimuth = pi * (3 / 2)
                }
                else {
                    azimuth = pi / 2
                }
            }

            # 計算した天頂角と方位角に対応する正距円筒射影の画像上の位置
            x.source <- round(azimuth / (2 * pi) * pxl.x) + 1
            y.source <- round(zenith.angle / pi  * pxl.y) + 1

            # 元画像からはみ出さないように。念のため。
            if (x.source > pxl.x || y.source > pxl.y) {
                next
            }

            # 変換先の座標（原点の移動）
            x.projected <- x + radius + margin.x
            y.projected <- y + radius + margin.y

            # 作成する画像に元画像の画素の内容をコピー
            equidistant.img[y.projected, x.projected, ] <- source.img[y.source, x.source, ]
        }
    }

    # 地平線を点線で書き込む。
    if (draw.horizon == TRUE) {
        n.dots <- 720      # 円周上の点の数
        dot.color <- c(0.5, 0.5, 0.5) # 色の指定 (R, G, B)
        for (i in 1:n.dots) {
            deg <- i * 360 / n.dots
            x = round(radius * sin(deg * pi / 180)) + radius + margin.x
            y = round(radius * cos(deg * pi / 180)) + radius + margin.y
            equidistant.img[y, x, ] <- dot.color # 灰色の点
        }
    }

    new.file.name <- sprintf("./%s/%s%s", subdir, eqd.prefix, src.file);
    writeJPEG(equidistant.img, new.file.name)
}