CopernicusでRicohTHETAのFisheye画像をEquirectangular画像に変換する

Houdini Advent Calendar 2024 7日目 の記事です

はじめに

記事の概要

以下の記事でnumPyとOpenCVを使って記述されている画像変換を、簡略化しつつHoudiniのCopernicusで動くように実装し直した記録になります。

引用は、すべてこの記事からのものです。

モチベーション

チャレンジした理由としては、自前のHDRIを手軽に生成してみたいと思ったからです。
CGSLAB LLC.さんによるCCWORLDの記事に触発されて始めたものでした。

記事の通りにRicohTHETAのRaw画像をdaktableでHDRI合成すると、上段のような画像が生成されます。
ダイナミックレンジは確保できたものの、このままだとCGで扱いにくい画像形式(DualFisheye)です。

この画像をCGで扱いやすくするためには、下段の画像形式(Equirectangular)に変換する必要があります。
CGSLAB LLC.さんの記事ではBlenderで行っていますが、それをHoudiniで、それもCopernicusでやってみたいと思いました。
Houdiniで実装できれば、PDGで簡単に量産できるのも夢があるなと思いました。

しかし、とりあえず実装までできた矢先に
撮影時にHDRIまで生成してくれちゃうプラグインがあることを知り、PDGまで実装して量産体制を整えるモチベーションがなくなってしまいました。
というわけで、道半ばですが学習記録として残しておこうと思います。

免責

あくまで簡易的な実装です

今回挑戦する画像変換はあくまで簡易的なもので、レンズプロファイルを使った完璧な結果を得るものではありません。

OpenCL部分のご利用は自己責任でおなしゃす

OSごとクラッシュすることもあるというOpenCL、初めて書いたので何が起こるかわかりません。

用語解説

• Dualfisheye形式
  ２個の180度半天球映像を横に並べたフォーマットです。Theta Sからはこの形式で出力されます。
• Equirectangular形式
  日本語では正距円筒図法といわれるもので、360度の映像が世界地図のメルカトル図法のように展開（正確にはちょっと違う）された形式です。

動作環境

Houdini Apprentice 20.5.278
NVIDIA GeForce RTX2070
Ricoh THEATA Z1

仮組みまで

全部AIに任せられないか試す

ソースコードをAIに投げて、"これをOpenCLに書き換えてください" とやってみました。
残念ながらうまくいきませんでした。
numPyは触ったことありましたが、OpenCVもOpenCL全くの初見です。
AI側もリリースされたばかりのHoudiniの新機能であるCopernicusのことをよく知りません。
この状態でちゃんと動くものを作れるわけないですよね。(2024年時点では)

AIにコード解説をさせる

ソースコードを全部貼り付けて、ブロックごとに解説してもらいました。

AIによる解説(クリックで展開)

ソースコード+解説

# ライブラリインポート
import sys
import numpy as np
import cv2

# エラーハンドリング (マップファイルが存在しない場合の再生成)
try:
    # 既存のマップファイルがあれば読み込む (高速化)
    xmap = np.load("xmap.npy")
    ymap = np.load("ymap.npy")
except IOError:
    # マップファイルが存在しない場合、新規生成
    print("マップファイルを生成します...")

    # 出力画像サイズの設定
    COLS = 1280
    ROWS = 720

    # 変換マップ初期化 (float32型で)
    xmap = np.zeros((ROWS, COLS), np.float32)
    ymap = np.zeros((ROWS, COLS), np.float32)

    # 出力画像の幅と高さを取得
    DST_X = float(COLS)
    DST_Y = DST_X / 2

    # カメラ1の中心 (画像幅の1/4)
    SRC_CX1 = DST_X / 4
    # カメラ2の中心 (画像幅 - カメラ1の中心)
    SRC_CX2 = DST_X - SRC_CX1

    # カメラの中心からの距離 (画像幅の1/4 * 0.884)
    SRC_CY = DST_X / 4
    SRC_R = 0.884 * DST_X / 4

    # カメラ1, 2の半径 (それぞれ補正係数をかける)
    SRC_RX = SRC_R * 1.00
    SRC_RY = SRC_R * 1.00

    # ループ処理で各ピクセルに対する変換座標を計算
    for y in range(COLS // 2):
        for x in range(COLS):
            # 出力画像の極座標 (θ, φ) を計算
            ph1 = np.pi * x / DST_Y
            th1 = np.pi * y / DST_Y

            # カメラ1の球面座標 (x1, y1, z1) を計算
            x1 = np.sin(th1) * np.cos(ph1)
            y1 = np.sin(th1) * np.sin(ph1)
            z1 = np.cos(th1)

            # 逆球面変換でカメラ2の極座標 (ph2, th2) を求める
            ph2 = np.arccos(-x1)
            th2 = (1 if y1 >= 0 else -1) * np.arccos(-z1 / np.sqrt(y1 * y1 + z1 * z1))

            # カメラ2の領域判定 (画像の半分を超えるか)
            if ph2 < np.pi / 2:
                r0 = ph2 / (np.pi / 2)  # 正規化された半径
                xmap[y, x] = SRC_RX * r0 * np.cos(th2) + SRC_CX1  # x座標の変換
                ymap[y, x] = SRC_RY * r0 * np.sin(th2) + SRC_CY  # y座標の変換
            else:
                r0 = (np.pi - ph2) / (np.pi / 2)  # 正規化された半径
                xmap[y, x] = SRC_RX * r0 * np.cos(np.pi - th2) + SRC_CX2  # x座標の変換
                ymap[y, x] = SRC_RY * r0 * np.sin(np.pi - th2) + SRC_CY  # y座標の変換

    # 生成したマップファイルを保存 (次回以降の高速化のため)
    np.save("xmap.npy", xmap)
    np.save("ymap.npy", ymap)


# デュアル魚眼画像を全方位投影に変換する関数
def convert_dualfisheye_to_equirectangular(frame):
    # OpenCV の remap 関数を使って変換 (双線形補間、境界は一定値)
    return cv2.remap(frame, xmap, ymap, cv2.INTER_LINEAR, cv2.BORDER_CONSTANT)

# スクリプト実行時の処理
if __name__ == '__main__':
    if len(sys.argv) < 3:
        print('使い方: Python3 dualfisheye.py <入力ファイル> <出力ファイル>')
        sys.exit(1)

    # 入力動画ファイルを開く
    cap = cv2.VideoCapture(sys.argv[1])
    if not cap.isOpened():
        print('ファイルが開けません')
        sys.exit(1)

    # 動画のフレームレート、高さ、幅を取得
    fps = cap.get(cv2.CAP_PROP_FPS)
    height = cap.get(cv2.CAP_PROP_FRAME_HEIGHT)
    width = cap.get(cv2.CAP_PROP_FRAME_WIDTH)

    # 出力動画ファイルの設定 (MPEG-4形式)
    fourcc = cv2.VideoWriter_fourcc('m', 'p', '4', 'v')
    out = cv2.VideoWriter(sys.argv[2], fourcc, fps, (int(width), int(height)))

    # フレーム処理のループ
    count = 0
    while cap.isOpened():
        # フレームを読み込む
        ret, frame = cap.read()
        if not ret:
            break

        # 魚眼画像から全方位画像に変換
        frame = convert_dualfisheye_to_equirectangular(frame)

        # 変換後のフレームを出力動画に書き込む
        out.write(frame)

        # 進捗表示 (30フレームごとに.)を出力
        if count % 30 == 0:
            print('.', end="")
            sys.stdout.flush()
        count += 1

    # 処理終了
    print()
    cap.release()
    cv2.destroyAllWindows()

Qiitaの記事も合わせ読みして、大まかに以下のような処理の流れを理解しました。

• 画像の情報をインポート、必要な変数を設定
• 座標変換"だけ"を出力画像のピクセルごとにループ処理
• 座標変換情報を元に変換先となる入力画像のピクセル情報を取得

大まかに処理の流れを掴んだら、Houdini上でどうやって楽に実装するか考えました。

• 画像のI/Oなどは必要ないので簡略化できそう
• FisheyeからEquirectangularへの変換は三角関数を用いた座標変換を行う必要があり、コーディングが必要そう
• "Single"Fisheyeならコーディングが楽になりそう
• DualFisheye→"Single"Fisheyeはノードで簡単にできそう
• ループ処理の考え方は変える必要がありそう

これらの気づきを得られたところで、実際にネットワークを組んでみました。

COPネットワーク仮組み

• インプットはDualFisheyeのEXR画像
• プレビューが暗いのでBright挟む
• Cropで右半分と左半分に分割
• OpenCLでEquirectangularに変換(予定)
• ContactSheetで再合成

OpenCLノードごとリファレンスコピー

↑のOpenCLノードですが、左右の画像で同じ処理をするつもりです。
処理を見直すたびにコードをコピペししたくないので、OpenCLノードごとリファレンスコピーします。
これで左右でコードが同期しました。
なんて快適なんでしょうか。
これがワンアクションで(Ctrl+Shift+Alt+ドラッグ)できてしまうの、本当に最高です。
Houdiniは文化だ。

OpenCL書く準備

基礎の基礎をインプット

• 堀川先生のyoutube いつも助けられています、ありがとうございます。
• SideFX公式マニュアル VEXユーザー視点で書かれていてわかりやすいです。

このあたりを流し見して、OpenCLノードのことを少し学びました。

• 書式はVEXとかなり似てる(どちらもCベース言語だからそれはそう)
• ピクセル毎に並列計算するらしい
• ピクセル情報は、@xres,@yres,@ix,@iyでアクセスできる
• I/OをSignatureタブで指定しておいてコードの中では@hogeでアクセスする

このぐらいは抑えておきました。

処理速度を上げる工夫

OpenCLとはいえフルレゾで処理するとCookに2-3秒かかっていたので、
ソース画像を読み込んだあとにResampleノードでレゾを落とし、デバッグ時のレスポンスを上げておきました。

デバッグ方法

値のデバッグは、単純に出力画像に乗せてInspectする方法を取りました。

float4 debug_color = ( float4 )( hoge, hoge, hoge, 1.0f );
@dst.set(debug_color);

InspectはCompositeWindow上で右クリックすると出てきます。

実際にソースコードをOpenCLに書き換える

座標変換情報の格納とループ処理

細かい座標変換部分はさておき、まずは何をやってるのかを見ています。

ソースコード+AI解説

# エラーハンドリング (マップファイルが存在しない場合の再生成)
try:
    # 既存のマップファイルがあれば読み込む (高速化)
    xmap = np.load("xmap.npy")
    ymap = np.load("ymap.npy")
except IOError:
    # マップファイルが存在しない場合、新規生成
    print("マップファイルを生成します...")
    
    ----省略----
    
    for y in range(COLS // 2):
        for x in range(COLS):
    
            ----省略----
            
            xmap[y, x] = ----省略----
            ymap[y, x] = ----省略----
    
    # 生成したマップファイルを保存 (次回以降の高速化のため)
    np.save("xmap.npy", xmap)
    np.save("ymap.npy", ymap)

def convert_dualfisheye_to_equirectangular(frmae):
    return cv2.remap(frame, xmap, ymap, cv2.INTER_LINEAR, cv2.BORDER_CONSTANT)

座標変換情報を一フレーム目で格納して、後フレームでは再利用することで処理負荷を抑えているようです。

DualfisheyeからEquirectangularへの変換するには、三角関数を駆使してEquirectangular上のピクセルがDualfisheye上のどのピクセルに相当するかを計算する座標変換処理を行うのですが、高速化のためにあらかじめ座標変換行列を作成しておき、OpenCVのremap関数を使って１フレーム分を一気に変換する方法を使います。

今回の目的としては、とりあえずCopernicusで変換できさえすれば御の字と考えました。
なので座標変換情報は、配列として格納はせずにそのままピクセル毎の並列計算でやってしまうことにしました。

最初の変数設定

方針が決まったところで、早速変数設定から書き換えていきます。

ソースコード+AI解説

# 出力画像サイズの設定
COLS = 1280
ROWS = 720

# 変換マップ初期化 (float32型で)
xmap = np.zeros((ROWS, COLS), np.float32)
ymap = np.zeros((ROWS, COLS), np.float32)

# 出力画像の幅と高さを取得
DST_X = float(COLS)
DST_Y = DST_X / 2

# カメラ1の中心 (画像幅の1/4)
SRC_CX1 = DST_X / 4
# カメラ2の中心 (画像幅 - カメラ1の中心)
SRC_CX2 = DST_X - SRC_CX1

# カメラの中心からの距離 (画像幅の1/4 * 0.884)
SRC_CY = DST_X / 4
SRC_R = 0.884 * DST_X / 4

# カメラ1, 2の半径 (それぞれ補正係数をかける)
SRC_RX = SRC_R * 1.00
SRC_RY = SRC_R * 1.00

出力画像の解像度やアスペクト比は、変換後にいくらでも変更できるので入力画像と同じままにします。
よって、COLS,ROWSは@xres,@yresに書き換えられます。
xmap,ymapはソースコードでは配列になっていますが、ピクセルごとに並列計算するためのインデックスの@ix,@iyがあれば今回は十分です。
SingleFisheyeの中心座標は、@xres,@yresを1/2するだけなのでわざわざ変数にしなくても良いですね。
SingleFisheyeにしたので、半径の係数は1/4ではなく1/2になります。
補正係数も特に必要なさそうなのでスルーします。必要になってもノードできそうですしね。

OpenCL

float xres = ( float )( @xres );
float yres = ( float )( @yres );
float ix = ( float )( @ix );
float iy = ( float )( @iy );
float fish_RADIUS = 0.884f * xres / 2.0f;

だいぶシンプルになりました。

座標変換その1

ソースコード+AI解説

# 出力画像の極座標 (θ, φ) を計算
ph1 = np.pi * x / DST_Y
th1 = np.pi * y / DST_Y

# カメラ1の球面座標 (x1, y1, z1) を計算
x1 = np.sin(th1) * np.cos(ph1)
y1 = np.sin(th1) * np.sin(ph1)
z1 = np.cos(th1)

出力となるEquirectangular画像の平面座標系を極座標系で評価しなおして、仮想球体の三次元座標にマッピングしてるようです。
この記事によれば、原点や軸の設定によって式が変化するようですが、一旦そのまま書き換えてしまいます。

OpenCL

float PI = ( float )( 3.141592653589793);
float eqr_PHI = PI * ix / xres;
float eqr_THETA = PI * iy / yres;
float sphere_x = sin( eqr_THETA ) * cos( eqr_PHI );
float sphere_y = sin( eqr_THETA ) * sin( eqr_PHI );
float sphere_z = cos( eqr_THETA );

元の変数名が抽象的だったので接頭語を付けました。

• eqr：出力となるEquirectangular画像を評価するもの
• sphere：仮想球体を評価するもの

PIには予約語があるかもしれないですが、呼び出し方を知らないのでとりあえず適当に定義しときました。

座標変換その2

ソースコード+AI解説

# 逆球面変換でカメラ2の極座標 (ph2, th2) を求める
ph2 = np.arccos(-x1)
th2 = (1 if y1 >= 0 else -1) * np.arccos(-z1 / np.sqrt(y1 * y1 + z1 * z1))

その1で作った仮想球体の三次元座標に対応する、Fisheye座標系における緯度と経度のようなものを求めているようですね。
この記事によれば、今回素材撮影に使用したTHETA Z1は等距離射影方式で画像が記録されているそうです。
等距離射影方式は球体をただ正投影したものに比べると、端の画像のピクセルの歪みが少ないようです。
極座標系が2種類出てきて混乱しますが、あまり深く考えずに書き換えを進めます。

また、仮想球体のy座標で緯度の正負を反転しているようです。

OpenCL

float fish_PHI = acos( -sphere_x );
float fish_THETA  = -acos( -sphere_z / sqrt( sphere_y * sphere_y + sphere_z * sphere_z ));
if ( shpere_y > 0 )
    {
    fish_THETA = -fish_THETA;
    }

前項と同じく接頭語を付けました。

• fish：入力となるFisheye画像を評価するもの

座標変換その3

ソースコード+AI解説

# カメラ2の領域判定 (画像の半分を超えるか)
if ph2 < np.pi / 2:
    r0 = ph2 / (np.pi / 2)  # 正規化された半径
    xmap[y, x] = SRC_RX * r0 * np.cos(th2) + SRC_CX1  # x座標の変換
    ymap[y, x] = SRC_RY * r0 * np.sin(th2) + SRC_CY  # y座標の変換
else:
    r0 = (np.pi - ph2) / (np.pi / 2)  # 正規化された半径
    xmap[y, x] = SRC_RX * r0 * np.cos(np.pi - th2) + SRC_CX2  # x座標の変換
    ymap[y, x] = SRC_RY * r0 * np.sin(np.pi - th2) + SRC_CY  # y座標の変換

ここが最後の座標変換で、Fisheye画像のどのピクセルをサンプルするべきかを求めているようです。
仮想球体の緯度ph2で条件分岐していて一瞬イメージできなくて戸惑いますが、分岐の差分を見るとdualfisheyeの左右のことだとわかります。
"single"fisheyeである今回の実装では必要ないので、とりあえず片っぽだけ使います。

OpenCL

float r0 = (fish_PHI) / PI * 2;
float fish_dest_x = ( float ) fish_RADIUS * r0 * cos( eqr_dest_THETA ) + xres / 2.0f;
float fish_dest_y = ( float ) fish_RADIUS * r0 * sin( eqr_dest_THETA ) + yres / 2.0f;

変数定義の時に確認したように"Single"Fisheyeの中心点はただの画像の中心なので、xres/2.0f,yres/2.0fと書き換えてあります。

ピクセル補間と出力

用意した座標を元に、Fisheye画像からピクセルをサンプルします。

ソースコード+AI解説

# デュアル魚眼画像を全方位投影に変換する関数
def convert_dualfisheye_to_equirectangular(frame):
    # OpenCV の remap 関数を使って変換 (双線形補間、境界は一定値)
    return cv2.remap(frame, xmap, ymap, cv2.INTER_LINEAR, cv2.BORDER_CONSTANT)

OpenCVのremap関数については元の記事に解説がありました。

remap関数はOpenCVに予め定義されている関数です。カメラのレンズの歪み補正などのためによく使われる関数だそうです。map1とmap2は、変換後の座標(xd, yd)をもとに変換前の座標(xs, ys)値を出力するテーブルです。逆変換のテーブルを用意しておくことがポイントです。
dst = cv.remap( src, map1, map2, interpolation)
src: 入力画像のフレームデータ
map1: 変換後の座標(x_d, y_d)を与えると変換前のx座標値(x_s)を出力するテーブル
map2: 変換後の座標(x_d, y_d)を与えると変換前のy座標値(y_s)を出力するテーブル
interpolation: 補完方法
dst: 出力画像のフレームデータ

補間しつつピクセルをサンプルする関数ということがわかりました。
OpenCLで似たような関数を探します。
SideFXのマニュアルで"Interpolate"とページ内検索したら、以下の関数が見つかりました。

@name.bufferSample(xy)
Bilinearly interpolates buffer elements at floating point buffer coordinates.
浮動小数点バッファ座標でバッファ要素を双線形補間します。

補間方法も同じなので多分これを使えば同等の処理になるでしょう。
ピクセル毎に処理するので、入力は配列ではなく座標変換その3で作成した座標にすればいいですね。

OpenCL

float2 fish_dest = ( float2 )( fish_x, fish_y );
float4 dest_color = ( float4 ) @src.bufferSample( fish_dest );
@dst.set( dest_color );

最後に出力画像として出すアトリビュート？@dstに取得したピクセル情報をsetして終わりです。

つまづいたポイント

ここまでの実装をつなげて動かしてみましたが、いい結果は得られませんでした。

今までのコードをつなげたもの

OpenCL

float xres = ( float )( @xres );
float yres = ( float )( @yres );
float ix = ( float )( @ix );
float iy = ( float )( @iy );
float fish_RADIUS = 0.94f * xres / 2.0f;
float PI = ( float )( 3.141592653589793);
float eqr_PHI = PI * ix / xres;
float eqr_THETA = PI * iy / yres;
float sphere_x = sin( eqr_THETA ) * cos( eqr_PHI );
float sphere_y = sin( eqr_THETA ) * sin( eqr_PHI );
float sphere_z = cos( eqr_THETA );

float fish_PHI = acos( -sphere_x );
float fish_THETA  = -acos( -sphere_z / sqrt( sphere_y * sphere_y + sphere_z * sphere_z ));
if ( sphere_y > 0 )
    {
    fish_THETA = -fish_THETA;
    }

float r0 = (fish_PHI) / PI * 2;
float fish_dest_x = ( float ) fish_RADIUS * r0 * cos( eqr_dest_THETA ) + xres / 2.0f;
float fish_dest_y = ( float ) fish_RADIUS * r0 * sin( eqr_dest_THETA ) + yres / 2.0f;

float2 fish_dest = ( float2 )( fish_dest_x, fish_dest_y );
float4 dest_color = ( float4 ) @src.bufferSample( fish_dest );
@dst.set( dest_color );

理想の結果が得られるように、一つずつエラーを潰していきます。

画像の向きが合わない

出力結果は90度回転しているように見えます。
タイプミスを探したり入念にデバッグしても、コードの書き換え自体にミスはありませんでした。
しかし、座標変換の部分のどこかで予期しないことが起きているのは確かです。
落ち着いて元の記事をよく見ると、元のDualFisheye画像は向きが違うことに気が付きました。

カメラの保存時の挙動が違うのか、空の向きが違いますね。
なるほど、座標系がずれてしまうわけですね。

座標変換の最後でFisheye画像を評価するXとYを入れ替えてみました。

修正前

float fish_x = ( float ) fish_RADIUS * r0 * cos( fish_THETA ) + xres / 2.0f;
float fish_y = ( float ) fish_RADIUS * r0 * sin( fish_THETA ) + yres / 2.0f;

↓

修正後

float fish_y = ( float ) fish_RADIUS * r0 * cos( fish_THETA ) + xres / 2.0f;
float fish_x = ( float ) fish_RADIUS * r0 * sin( fish_THETA ) + yres / 2.0f;

だいぶ近づきましたが、まだエラーがあるようです。

水平方向に回転している

ヘリポート腰掛けている人(僕)に注目すると、左にずれてしまっていることがわかります。
角度にすると、ちょうど90度に見えます。
Dual,SingleFisheye双方のEquirectangular変換を解説しているページによれば、DualFishEyeと比較してSingleFisheyeにはπ/2のオフセットが入っているので90度回せばよさそうです。
オフセットは出力ピクセルの経度の部分にかかっているので、真似してみます。

修正前

float fish_PHI = PI * ix / xres;

↓

修正後

float fish_PHI = PI * ix / xres - PI / 2;

ついに繋がりました！ 感動の瞬間です。
ここまでたどり着くのに、何気に1週間ぐらいかかっていました。

サイズ調整

ReampleノードをBypassしてフルレゾで確認してみると、つなぎ目が目立って見えました。
つなぎ目付近に写り込んでいるものから特徴的な場所を探して、唯一残した調整パラメータFisheyeRADIUSの最適な値を探りました。

つなぎはだいぶ緩和したかなと思います。
今回は実装するまでが目的なので、このぐらいのクオリティで良しとします。

nanピクセル回避

縦に白いピクセルがあるにも気が付きました。

Inspectしてみるとnanでした。

どうやらゼロ除算が発生してnanピクセルができてしまっているようです。
例外処理書いて回避します。

nanピクセル回避

    // If the point is very close to the horizon (|sphere_y| is nearly zero), 
    // handle the special case for the horizon.
     if ( fabs(sphere_y) < 1e-3 )
        {
        if ( sphere_z > 0 )
            {
            fish_dest_THETA = PI;
            }
        if ( sphere_z <= 0 )
            {
            fish_dest_THETA =  0;
            }
        }

結果

というわけで、どうにかまともに動くようになりました。

初めての道具で作ったものが動いたときの感動はひとしおですね。
OpenCLがめちゃくちゃ速いのも非常に驚きました。
記録のために一応貼っておきます。誰か(もしくは未来の自分)の助けになれば幸いです。

COPネットワーク

ほとんど仮組みと変わってないです。

Resampleのあとに入っているOpenCLは、Brightノードの代わり↓です。

OpenCL

#bind layer src? val=0
#bind layer !&dst

@KERNEL
{
    @dst.set(@src*1000);
}

Fisheye to Equirectangular

Fisheye to Equirectangular

#bind layer src? val=0
#bind layer !&dst

@KERNEL
{
    // Debug colors for visualization
    float4 debug_color0 = ( float4 )( 1.0f,0.0f,0.0f,1.0f );

    // Get input image dimensions
    float xres = ( float )( @xres );
    float yres = ( float )( @yres );
    float ix = ( float )( @ix );
    float iy = ( float )( @iy );
    float4 debug_color2 = ( float4 )( ix / xres, iy / yres, 0.0f, 1.0f );
    
    // Define PI for calculations
    float PI = ( float )( 3.141592653589793);
    
    // Define the radius of the fisheye projection
    float fish_RADIUS = 0.94f * xres / 2.0f;
    
    // Calculate polar coordinates on the fisheye sphere
    float eqr_PHI = PI * ix / xres + PI /2;
    float eqr_THETA = PI * iy / yres;
    float4 debug_color3 = ( float4 )( eqr_PHI, eqr_THETA, 0.0f, 0.0f );

    // Convert spherical coordinates to equirectangular coordinates
    float sphere_x = sin( eqr_THETA ) * cos( eqr_PHI );
    float sphere_y = sin( eqr_THETA ) * sin( eqr_PHI );
    float sphere_z = cos( eqr_THETA );

    
    // Convert equirectangular coordinates back to polar coordinates
    float fish_PHI = acos( -sphere_x );
    float fish_THETA  = -acos( -sphere_z / sqrt( sphere_y * sphere_y + sphere_z * sphere_z ));

    // Handle edge cases
    // - If the point is below the horizon (sphere_y > 0), flip the theta angle
    if ( sphere_y > 0 )
        {
        fish_THETA = -fish_THETA;
        }

    // If the point is very close to the horizon (|sphere_y| is nearly zero), 
    // handle the special case for the horizon.
     if ( fabs(sphere_y) < 1e-3 )
        {
        if ( sphere_z > 0 )
            {
            fish_THETA = PI;
            }
        if ( sphere_z <= 0 )
            {
            fish_THETA =  0;
            }
        }
        
    float4 debug_color5 = ( float4 )( fish_PHI, fish_THETA, 0.0f, 0.0f );
    
    // Convert polar coordinates to 2D image coordinates
    float r0 = fish_PHI / PI * 2;
    float fish_y = ( float ) fish_RADIUS * r0 * cos( fish_THETA ) + xres / 2.0f;
    float fish_x = -( float ) fish_RADIUS * r0 * sin( fish_THETA ) + yres / 2.0f;
    float2 fish_dst = ( float2 )( fish_x, fish_y );
    float4 debug_color6 = ( float4 )( fish_x, fish_y, 0.0f, 0.0f );

    // Sample the source image at the calculated destination coordinates
    float4 dest_color = ( float4 ) @src.bufferSample( fish_dst );
    
    @dst.set( dest_color );

    // Enable debug colors
    //@dst.set(debug_color0); // single color
    //@dst.set(debug_color1); // 
    //@dst.set(debug_color2); // UV color
    //@dst.set(debug_color3); // ( eqr_PHI, eqr_THETA, 0.0f, 0.0f )
    //@dst.set(debug_color4); // ( sphere_x, sphere_y, sphere_z, 0.0f )
    //@dst.set(debug_color5); // ( fish_PHI, fish_THETA, 0.0f, 0.0f )
    //@dst.set(debug_color6); // ( fish_x, fish_y, 0.0f, 0.0f )
}

余談

実は最初はVEXで実装した

レガシーCOP(2.0)上でVEXを使ったピクセル処理は経験があったので、最初はVEXで実装を試しました。
文法はわかっていて、数式を書き換えるだけだったので結構
すぐに実装できました。
しかし、絶望的にCookに時間がかかったのですぐに諦めました。やっぱりGPUってすごいですね。

適当にいじって泥沼化した

つまづいた時、記事に書いたような正解?にすぐたどり着いたわけではありませんでした。
コードの内容を深く理解しないままに、出力結果を見ながらガチャガチャと符号を変えたりオフセットし続けました。
このガチャガチャに3-4日分のアフター5が消えました。
結果的には正解っぽい挙動を見つけて、それを裏付ける資料はあとから(記事にまとめながら)見つけて理解を深めていきました。
いつもこんな感じでめちゃ効率悪いのを何とかしたいんですが、この辺がアーティスト上がり?の僕の限界なのかなと思ったり。

～ですね

記事の文体は元同僚のYさんを参考にしました。
彼の日報を読むのが好きでした。諸々リスペクトを込めて。