めたへいログ

GIS系エンジニアの備忘録

IDLの積演算で行われる型変換について勘違いしていた件

公私ともにIDLという言語を使用している。
これは、CやFortranに影響を受けたプログラミング言語で、さっと使えるため、ちょっとした計算や処理はこれで書くことが多い。

www.harrisgeospatial.co.jp

今回、これでプログラムを書いていた際におかしな現象を確認したので、備忘録として残しておく。

大容量のファイルを読み込みたい

ある画像データを考える。最近だと広い範囲をセンチメートルオーダの高解像度のドローンで取るようなことも珍しくないので、ファイル容量は近年増加傾向にある。
こういった、時には数十GBにも達するような大容量のファイルを一度に読み込んで変数に設定する場合、処理に大きな負荷をかけることになるので好ましくない。

data = read_tiff("/path/to/bigSizeTiffData.tif") ; メモリに負荷がかかる。。

このような場合にはREAD_BINARYなどの関数を使用して、特定の範囲だけにアクセスすることが推奨される。
この関数のキーワードではデータタイプや読み込み開始位置を指定することができるので、ファイルへの柔軟なアクセスが行え、大きなファイルの処理の際にはよく使われている。

www.l3harrisgeospatial.com

たとえば、ある3次元のデータを考える。

  • XYZのそれぞれの次元が2000個の分解能を持つ
  • ある点(Xi,Yi,Zi)の点にFLOAT型で定義された浮動小数の値が定義されている

この場合のファイルサイズをバイト単位で計算すると、2000*2000*2000*4(FLOATは4バイトのサイズを占める) = 32,000,000,000で32GBの容量になる。

これを一気に読み込むのは現実的でない。
もし処理したい範囲がXが100の地点から10ポイント、Yが100の地点から10ポイント、Zが100の地点から10ポイントの都合1000個分だけの範囲に収まるなら、以下のような書き方ができる。
(※実際にはどのような順序で形でデータが格納されているかまで考慮する必要があるが、ここではあまり厳密には触れないものとする)

pnt = ULONG64(100)*100*100*4

readData = READ_BINARY($
  "/path/to/bigSize3dData.data", $
  ENDIAN="big", $
  DATA_TYPE=4, $ ; floatを示す番号
  DATA_START= pnt, $ ; 読み込み開始地点
  DATA_DIMS=[10,10,10]) ; 指定の場所から読んで[10,10,10]で取得
  help, readData ; readDataが[10,10,10]の配列と確認できる

今回何があったか。

上記のような方法で読み込もうとするときに、pntの計算でおかしなことが起きた。
読み込み位置の起点となるxiyiziを指定して読み込み開始位置を計算する際に、以下のような関数を用意していたが、この戻り値がどうもおかしい。

function get_xyz_point, x, y, z
    compile_opt idl2

    xmax = 2000
    ymax = 2000
    zmax = 2000
    floatCoff = LONG64(4)

    pnt = 0
    pnt += z * xmax * ymax * floatCoff
    pnt += y * ymax * floatCoff
    pnt += x * floatCoff

    return, pnt

end

途中まではうまく動いていたように思ったのだが、以下のように設定した際の戻り値が奇妙になっていた。

pnt = get_xyz_point(500, 500, 800)
print, pnt 
; answer > -4,379,869,184 であり、意図したものでない。。

デバッグで気づいたのだが、以下の計算結果がおかしくなっていた。

pnt += z * xmax * ymax * floatCoff 

どういった理由かというと、floatCoffでLONG64型を指定しているのでまず全てがLONG64型にスケールされて計算されるかと思っていたのだが、違ったらしい。
正しくは、先頭から順番に積を取っていき、その都度一番大きな型にスケールされるというのがIDLでの動作だったようだ。

IDL> help, z
Z               INT       =      800  ; 引数のzはint型
IDL> help, xmax
XMAX            LONG      =         2000 ; 内部で定義したxmaxはlong型
IDL> help, ymax
YMAX            LONG      =         2000 ; 内部で定義したymaxはlong型

IDL> help, z * xmax
<Expression>    LONG      =      1600000 
; z* xmaxまでは int * longなのでlong型

IDL> help, z * xmax * ymax
<Expression>    LONG      =  -1094967296 
; ここで int * long * long で long型となってしまう。
; 意図した 800*2000*2000 = 3,200,000,000 はlongの範囲外

IDL> help, z * xmax * ymax * floatCoff
<Expression>    LONG64    =            -4379869184 
; 上記のlongの結果はすでにマイナス。
; このもとで long * long64 でlong64型にスケールされたが
; 前段階の z * xmax * ymax の時点で既に意図した結果にはなっていない。。

内部で定義する値についても明示的にLONG64で定義すればこの問題は避けることができた。
大容量のファイルを扱う場合には、常に型を意識しないと怖い目にあうと肝に銘じておく。

最終的な修正プログラムは以下の通り。

function get_xyz_point, x, y, z
    compile_opt idl2

    xmax = LONG64(2000)
    ymax = LONG64(2000)
    zmax = LONG64(2000)
    floatCoff = LONG64(4)

    pnt = 0
    pnt += z * xmax * ymax * floatCoff
    pnt += y * ymax * floatCoff
    pnt += x * floatCoff

    return, pnt

end

pnt = get_xyz_point(500, 500, 800)
print, pnt ; answer > 12,804,002,000 という意図したものになっている。