numpy (python) でボクセル値を高速に取得する方法を教えてください

Question

いま２つの配列 coordinate と voxel があるとします。
voxelは三次元配列 (512,512,512) でそれぞれにボクセル値を持っています。
coordinateは二次元配列 (512,512) で今から取得したいボクセル値の座標が三次元のリストとして格納されています。

このcoordinateは、いわゆる斜め掛けの座標を保持しているため、例えばvoxel[:,:,512] のようにして任意の座標軸のボクセル値を一度に取得できないと考えています。

私は以下のようなコードを書いて二次元配列 imgAry (512, 512) に任意の座標のボクセル値を格納させましたが、より高速なコードがあるなら教えて頂きたいです。よろしくお願いします。

import numpy as np
imgAry = np.array( [voxel[ j[0], j[1], j[2] ] for j in coordinate ] ).reshape(512, 512)

Comments

[kirara0048]

`coordinate`は、「リストを要素にもつndarray」ということですか？
`np.concatenate(coordinate).shape`は`(262144, 3)`になりますか？

`coordinate`の要素が`[x, y, z]`であるとき、それを`voxel[x, y, z]`の値に置換する、という認識であってますか？

[fuji_san]

ご質問ありがとうございます。
>> `coordinate`は、「リストを要素にもつndarray」ということですか？
おっしゃる通りです。
>> `np.concatenate(coordinate).shape`は`(262144, 3)`になりますか？
上図にミスがありました。coordinate.shape が (262144, 3)になります。
よって np.concatenate(coordinate).shape は (786432,) になります。
`coordinate`の要素が`[x, y, z]`であるとき、それを`voxel[x, y, z]`の値に置換する、という認識であってますか？
おっしゃる通りです。

[bsdfan]

多次元のndarrayを使うのではなく、リストを要素に持つndarrayを使う理由がありますか？
速度を重視するならndarrayだけでやる方が良さそうですが。

Answer 1

いまいち状況がよくわかっていませんが……

#
voxel[coordinate[:, 0], coordinate[:, 1], coordinate[:, 2]]

coordinateが (262144, 3)形状の整数データ型配列なのであれば、上記で (262144,)形状のndarrayが出力され、以下で求める形になるはずです。

voxel[coordinate[:, 0], coordinate[:, 1], coordinate[:, 2]].reshape(512, 512)

Answer 2

512の例を出すのは訳が分からなくなるので4に縮小します。各boxelの値は、RGBを想定しています。
512のときは、N=512にしてください。

まずサンプルデータを作ります。
voxelです。

>>> import numpy as np
>>> import random
>>>
>>> N = 4
>>> voxel = np.zeros(N*N*N, dtype=object)
>>> for i in range(N*N*N):
...     voxel[i] = [random.randint(0,255) for i in range(3)]
...
>>> voxel = voxel.reshape(N,N,N)
>>> print(voxel)
[[[list([149, 18, 49]) list([220, 1, 186]) list([168, 111, 239])
   list([105, 140, 107])]
  [list([110, 181, 29]) list([203, 197, 147]) list([185, 157, 219])
   list([133, 123, 114])]
  [list([198, 89, 152]) list([131, 62, 104]) list([12, 185, 234])
   list([88, 82, 84])]
  [list([138, 215, 119]) list([136, 117, 121]) list([158, 98, 121])
   list([13, 181, 209])]]

 [[list([172, 0, 153]) list([132, 39, 155]) list([105, 190, 101])
   list([240, 192, 199])]
  [list([38, 254, 24]) list([9, 69, 173]) list([211, 89, 210])
   list([223, 159, 43])]
  [list([42, 92, 125]) list([101, 194, 24]) list([134, 215, 16])
   list([54, 176, 51])]
  [list([168, 116, 18]) list([94, 81, 49]) list([192, 41, 45])
   list([137, 182, 178])]]

 [[list([31, 68, 194]) list([174, 1, 104]) list([236, 63, 160])
   list([19, 20, 44])]
  [list([238, 60, 184]) list([172, 106, 34]) list([121, 173, 227])
   list([198, 145, 202])]
  [list([183, 225, 163]) list([93, 101, 74]) list([166, 39, 45])
   list([69, 169, 247])]
  [list([195, 47, 140]) list([205, 195, 218]) list([142, 240, 180])
   list([230, 50, 121])]]

 [[list([174, 5, 140]) list([157, 144, 74]) list([206, 86, 223])
   list([121, 40, 95])]
  [list([158, 68, 102]) list([3, 42, 87]) list([61, 76, 57])
   list([2, 77, 103])]
  [list([216, 241, 127]) list([175, 45, 225]) list([11, 153, 174])
   list([199, 98, 154])]
  [list([192, 174, 217]) list([12, 36, 70]) list([72, 164, 22])
   list([18, 121, 128])]]]

次にcoordinateを作ります。

>>> coordinate = np.zeros(N*N, dtype=object)
>>> for i in range(N*N):
...     coordinate[i] = [random.randint(0,N-1) for i in range(3)]
...
>>> coordinate = coordinate.reshape(N,N)
>>> print(coordinate)
[[list([1, 1, 1]) list([2, 0, 0]) list([0, 2, 3]) list([0, 2, 1])]
 [list([1, 1, 3]) list([2, 1, 2]) list([3, 3, 3]) list([0, 3, 2])]
 [list([0, 0, 1]) list([0, 0, 1]) list([1, 0, 3]) list([1, 2, 0])]
 [list([0, 1, 3]) list([0, 0, 1]) list([2, 3, 2]) list([1, 3, 2])]]

最後にimgAryを生成します。

>>> coordinate_a = np.array(list(coordinate.reshape(-1))).reshape(N,N,3)
>>> coordinate_s = (coordinate_a[:,:,0]*(N*N)+coordinate_a[:,:,1]*N +coordinate_a[:,:,2]).reshape(-1)
>>> imgAry = voxel.reshape(-1)[coordinate_s].reshape(N, N)
>>> print(imgAry)
[[list([9, 69, 173]) list([31, 68, 194]) list([88, 82, 84])
  list([131, 62, 104])]
 [list([223, 159, 43]) list([121, 173, 227]) list([18, 121, 128])
  list([158, 98, 121])]
 [list([220, 1, 186]) list([220, 1, 186]) list([240, 192, 199])
  list([42, 92, 125])]
 [list([133, 123, 114]) list([220, 1, 186]) list([142, 240, 180])
  list([192, 41, 45])]]

このやり方は、fuji_sanさんの方法よりはかなり速いと思います。
これ以上高速にしたければ、リストのndarrayではなく、高次元のndarrayを使うのが良いのですが、元データの作成や、imgAryの使い方によるところもあるので、なんとも言えないですね。

Comments

[fuji_san]

ご回答ありがとうございます。
ご提示の方法をColaboratoryで試してみましたが速度が２倍遅くなりました。実際に必要としているN=512ではご教授頂いた方法では難しいようです。

以下、試したコードです。ボクセル値はRGBではなくGrey値なので一部コードを改変しています。

```python
import numpy as np
import random

N = 512
voxel = np.zeros(N*N*N, dtype=object)
for i in range(N*N*N):
voxel[i] = random.randint(0,255)
voxel = voxel.reshape(N,N,N)

coordinate = np.zeros(N*N, dtype=object)
for i in range(N*N):
coordinate[i] = [random.randint(0,N-1) for i in range(3)]

def test1():
coordinate = coordinate.reshape(N,N)
coordinate_a = np.array( list(coordinate.reshape(-1))).reshape(N,N,3)
coordinate_s = (coordinate_a[:, :, 0]*(N*N) + coordinate_a[:, :, 1]*N + coordinate_a[:,:,2]).reshape(-1)
imgAry = voxel.reshape(-1)[coordinate_s].reshape(N, N)

def test2():
imgAry = np.array( [voxel[ j[0], j[1], j[2] ] for j in coordinate] ).reshape(N,N)

%timeit -n 1 test1()
%timeit -n 1 test2()

'''
1 loop, best of 5: 496 ms per loop
1 loop, best of 5: 233 ms per loop
'''
```