計算負荷の偏り
C/C++プログラミングで，OpenMPのお世話になっています。頭が上がりません 
画像処理をはじめ，配列を扱うプログラムでは，配列の各要素を並列計算することがよくありますよね。
例えば，2次元画像の場合，こんな感じです。
// 画像サイズ
int img_shape[2] = { 512, 512 };
// 画素数
int num_pixel = img_shape[0] * img_shape[1];
// 2次元画像
std::vector<float> img( num_pixel );
// バイナリファイルを読み込む（架空の関数です）
read_from_file( "image.f256x256.bin", img );
// 画素毎の並列処理
#pragma omp parallel for
for ( int n = 0; n < num_pixel; n++ ) {
// 多次元インデックス
int index[2] = { n%img_shape[0], n / img_shape[0] };
// 各画素の処理
// ...
}
ここで，値が正の画素のみを処理対象とする場合を考えます。
並列計算の部分はこんな感じになります。
#pragma omp parallel for
for ( int n = 0; n < num_pixel; n++ ) {
// 多次元インデックス
int index[2] = { n % img_shape[0], n / img_shape[0] };
// 各画素の処理
if ( img[n] > 0.f )
{
// ...
}
}
この何気ない実装には，実はとても大きな問題があります。
それは，処理対象の画素が一部の領域に偏っている場合に，並列計算の効率が落ちることです。
例えば，下図のように偏っている場合に，並列スレッド数を4にすると，スレッド0,1,2,3のうち，スレッド1とスレッド2だけに負荷が集中し，スレッド0とスレッド3は開店休業状態になり，演算能力を持て余してしまいます 
この記事では 「スレッド負荷を分散して本来の演算能力を引き出す方法」 を紹介します。
対処方法１
一番ナイーブな対処方法は，画素番号のシーケンシャルな配列index_map={0,1,2,...,num_pixel-1}
をシャッフルして，画素番号のマッピング（読み替え） を行うことです。
#pragma omp parallel for
for ( int k = 0; k < num_pixel; k++ ) {
// マッピング
int n = index_map[k];
// 多次元インデックス
int index[2] = { n % img_shape[0], n / img_shape[0] };
// 各画素の処理
if ( img[n] > 0.f )
{
// ...
}
}
これで計算負荷は分散できますが，この実装には，別の大きな問題があります。
それは，画素数と同じサイズのマッピング配列が必要であること，つまり，Memory-Efficient でないことです。
例えば，非破壊検査CTでは，2048 * 2048 * 2048といった巨大な3次元画像を扱う場合があり，その場合，マッピング配列(8バイト要素)のサイズは，なんと約70GBになります。流石にキツイ。
ということで，この対処方法は，巨大な画像を扱う場合，実用的ではありません。
対処方法２
シンプルでMemory-Efficientな対処方法があります。
軸毎にシャッフルされたマッピング配列を準備すればいいのです。
例えば，2次元画像であれば，横方向・縦方向のそれぞれで，マッピング配列index_map_x， index_map_yを準備すれば済むのです。
この方法では，マッピング配列のサイズは，各軸のサイズになります。2048 * 2048 * 2048の3次元画像であっても，たったの24KBです。余裕！
#pragma omp parallel for
for ( int k = 0; k < num_pixel; k++ ) {
// 多次元インデックス
int index[2] = { k % img_shape[0], k / img_shape[0] };
// マッピング
index[0] = index_map_x[index[0]];
index[1] = index_map_y[index[1]];
// 線形インデックス
int n = index[0] + index[1] * img_shape[0];
// 各画素の処理
if ( img[n] > 0.f )
{
// ...
}
}
この方法でも負荷分散でき，Memory-Efficientなので実用的でもあります。
この方法をふと思い付きましたが，極めてフツーな考え方なので，既に誰かがやってるはずです。
対処方法３
わざわざインデックスのマッピング配列を準備しなくても，マッピング関数を準備すれば，それで済みます。
例えば，配列サイズが15，並列スレッド数が4の場合は，このような「等間隔」なマッピングが考えられます。
{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14} --> {0,4,8,12,1,5,9,13,2,6,10,14,3,7,11}
つまり，f(0)=0, f(1)=4, ..., f(14)=11 を満たすマッピング関数f(･)が必要です。
以下は，(恐らく)任意の配列サイズNと任意の並列スレッド数Tに対応したマッピング関数です。
// マッピング関数
int equipitched_map( int index, int N, int T )
{
return
( N > T )
? ( ( index * T ) - ( ( index * T ) / N ) * ( N - ( N % T == 0 ) ) )
: index;
}
このマッピング関数を使うと，次のような実装に帰着します。
線形インデックスをマッピングする場合（対処方法１の改良版）
#pragma omp parallel for
for ( int k = 0; k < num_pixel; k++ ) {
// マッピング
int n = equipitched_map( k, num_pixel, omp_get_max_threads() );
// 多次元インデックス
int index[2] = { n % img_shape[0], n / img_shape[0] };
// 各画素の処理
if ( img[n] > 0.f )
{
// ...
}
}
多次元インデックスをマッピングする場合（対処方法２の改良版）
#pragma omp parallel for
for ( int k = 0; k < num_pixel; k++ ) {
// 多次元インデックス
int index[2] = { k % img_shape[0], k / img_shape[0] };
// マッピング
index[0] = equipitched_map( index[0], img_shape[0], omp_get_max_threads() );
index[1] = equipitched_map( index[1], img_shape[1], omp_get_max_threads() );
// 線形インデックス
int n = index[0] + index[1] * img_shape[0];
// 各画素の処理
if ( img[n] > 0.f )
{
// ...
}
}
最後に
なぜこんなことを考えるに至ったかというと，MySimuに3次元画像を距離変換するプログラムを実装していて，「なんか計算遅いな」となって，「前景の画素しか処理対象にならないからか？」と思って調べてみた，というような感じです。
下図が実例（画像サイズ：128x128x63）で，インデックスのマッピングを実装したら，計算時間は20秒から半減しました 
以上です