|
HPC/並列プログラミングポータルでは、HPC（High Performance Computing）プログラミングや並列プログラミングに関する情報を集積・発信しています。 |
［記事一覧を見る］
インテルのCoreシリーズCPUが登場して以来、マルチコアCPUは爆発的に普及している。現在では比較的低価格なCPUでも複数のコアを搭載するようになり、現在販売されているPCのほとんどはマルチコアに対応しているといっても過言ではない。しかし、まだマルチコアCPUの性能を生かす、並列処理を行っているアプリケーションは多くない。
並列処理は、一般には実装が難しい、という印象があるようだ。確かにスレッドを駆使して処理を並列化する場合、スレッドの管理やスレッド毎の連携など、考慮しなければならないことが増え面倒ではある。しかし、プログラムを並列化するための言語規格「OpenMP」や、C++用の並列化ライブラリ「Threading Building Blocks」といった並列化支援技術を利用することで、プログラムの並列化へのハードルは大幅に低くなる。また、インテルの開発製品「Parallel Studio」に含まれる「Parallel Composer」には、独自の並列化支援機能を備えたC/C++コンパイラーが含まれている。このような技術や製品を利用すれば、手軽に並列処理を実装することが可能だ。
本記事では、このような技術を利用した並列処理の実装方法を、簡単な画像処理アプリケーションを例に説明していこう。
今回サンプルに使用するアプリケーションは、画像に対して「メディアンフィルタ」と呼ばれるフィルタ処理を行うものだ（ソースコード全文は記事末にリスト1として掲載）。メディアンフィルタとは、ノイズ除去等を目的に利用される画像処理フィルタで、あるピクセルに対し、そのピクセルの周辺にあるピクセルの値の中央値をそのピクセルの値とするものである（図1）。これにより、図2のように単発的に発生するノイズ（ごましおノイズとも呼ばれる）を軽減する効果がある。
今回作成したアプリケーションでは、指定したカラーJPEG画像に対してメディアンフィルタを適用し、JPEG形式で保存するコマンドラインアプリケーションだ。対象画像ファイルおよびフィルタを適用後の画像ファイルはコマンドライン引数で指定する。処理の流れは図3のようになる（ソースコード全文は記事末にリスト1として掲載）。
さて、並列処理を実装するにあたり、もっともプリミティブな方法がスレッドを利用した実装である。Windowsの場合、新しいスレッドを作成する方法にはいくつかあるが、Cでスレッドを作成する場合は_beginthreadex関数を使用するのが基本となる。
_beginthreadex関数は「process.h」内で定義されており、そのプロトタイプは下記のようになっている。
uintptr_t _beginthreadex( void *security, unsigned stack_size, unsigned ( *start_address )( void * ), void *arglist, unsigned initflag, unsigned *thrdaddr );
詳細についてはライブラリリファレンス等を参照してほしいが、基本的にはstart_address引数に新たなスレッドで実行したい関数を、arglistに関数に引数として与える変数もしくは配列のポインタを与えて_beginthreadex()を呼び出せば良い。それ以外の引数については、特に問題がなければ0もしくはNULLで構わない。
さて、今回のプログラムではループを使って各画素に対するメディアン処理を繰り返し行っている。このように繰り返し回数が多いループを含むプログラムの場合、このループを分割して並列処理させるのが定石である（図4）。
まず、並列化したい個所を別の関数として分離し、この関数を別スレッドで実行することで並列処理を行うわけだ。この関数を別スレッドで実行する場合、関数には1つの引数しか与えられないので、複数の引数を与える場合はあらかじめそれを構造体として宣言しておき、その構造体のポインタを引数として与える。
今回のサンプルプログラムで並列実行したい個所は次の部分である。この部分は最大で4重のループとなっているが、このうちもっとも外側のループ（「for(y = 1; y height - 1; y++)」の部分）を分割することを考えよう。つまり、このループを「for(y = 1; y (height - 1)/2; y++)」と「for(y = (height - 1)/2; y height - 1; y++)」という2つのループに分割し、それぞれを並列に処理する、という流れである（並列化後のソースコード全文は記事末にリスト2として掲載）。
edit
/* フィルタを適用 */
for(y = 1; y height - 1; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3 * x + i] = buf_in-array[y+dy][3 * (x+dx) + i];
}
}
}
この処理内ではいくつかの変数が使われているが、この処理を別の関数にまとめる際に引数として与える必要があるものはbuf_tmp、buf_in、buf_outの3つである。そこで、この3つの変数に、ループを開始する値（y_start）および終了する値（y_end）を加えた5つの変数を関数に与える配列として宣言する。そのほかの変数についてはすべてこのループ内で初期化され、また出力にも影響しないので関数内で確保を行えば良い。
typedef struct {
image_buffer* buf_in;
image_buffer* buf_out;
image_buffer* buf_tmp;
int y_start;
int y_end;
} thread_param;
この構造体のポインタを引数として受け取り、並列に処理を実行する関数は次のようになる。なお、_beginthreadexの引数に与える関数は__stdcall形式でなければならない点に注意しよう。
void __stdcall work_thread(thread_param* param) {
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
const image_buffer* buf_in = param-buf_in;
const image_buffer* buf_out = param-buf_out;
const image_buffer* buf_tmp = param-buf_tmp;
const int y_start = param-y_start;
const int y_end = param-y_end;
const int width = param-buf_in-width;
/* フィルタを適用 */
for(y = y_start; y y_end; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3*x + i] = buf_in-array[y+dy][3*(x+dx) + i];
}
}
}
}
実際にスレッドを生成し、関数を呼び出す個所は次のようになる。
/* 与える引数を設定 */
for (i = 0; i 2; i++) {
param[i].buf_in = buf_in;
param[i].buf_out = buf_out;
param[i].buf_tmp = buf_tmp;
}
param[0].y_start = 1;
param[0].y_end = (height - 1) / 2;
param[1].y_start = (height - 1) / 2;
param[1].y_end = height - 1;
/* スレッド生成 */
thread = (HANDLE)_beginthreadex(NULL, 0, (void*)work_thread, (param[0]), 0, NULL);
work_thread((param[1]));
/* スレッドの終了を待つ */
WaitForSingleObject(thread, INFINITE);
CloseHandle(thread);
ここでは_beginthreadexでy=1からy=(height-1)/2までの処理を行う別スレッド（スレーブスレッド）を起動させた後、もともとあったスレッド（マスタースレッド）でy=(height-1)/2からy=height-1までの処理を実行するように関数を呼び出している。マスタースレッドでの処理が終了した後は、マスタースレッドはスレーブスレッドが終了するまで待機する（WaitForSingleObject関数）。スレーブスレッドの処理が終了したら、スレッドハンドルを閉じて処理を完了させる。
このようにして並列化を行ったプログラムを表2のような条件でコンパイルし実行したところ、並列化を行うことで約2倍近い高速化の効果が得られた（表1）。
| プログラム | 実行時間 |
|---|---|
| シングルスレッド版（med_serial.c）並列化前 | 5553ミリ秒 |
| マルチスレッドによる並列化版（med_thread.c） | 2870ミリ秒 |
| 要素 | スペック |
|---|---|
| CPU | Core 2 Duo E6550（2.33GHz） |
| メモリ | 2GB |
| OS | Windows Vista Business SP1 |
| コンパイラ | Visual C++ 2008 Professional Edition（SP1） |
| テストに使用した画像 | 4000×3000ピクセルのJPEG画像 |
今回並列化を行ったプログラムでは、このように高速化に成功しているように見えるものの、実はいくつかの問題を含んでいる。まず、今回の実装では最大で2スレッドしか利用しないため、3コア以上を持つシステムではそのパフォーマンスをフルに引き出すことはできない。もし3スレッド以上を同時に実行させたい場合、そのためのコードを別途追加する必要がある。また、それぞれのスレッドに割り当てる仕事量は固定されているため、スレッドに割り当てられた処理が終了したら、そのスレッドはアイドル状態となってしまう。そのほか実装上の問題として、マルチスレッドで処理する部分を別の関数として分離させる必要があるほか、_beginthreadex()関数では実行する関数の引数としてポインタを1つだけしか与えられないため、引数の処理が面倒となる点も挙げられる。
このような問題は、OpenMPを利用することで大幅に改善が可能だ。OpenMPはC/C++およびFortranで利用できる、並列処理を記述するための言語規格である。OpenMPによる並列化では、並列化を行う個所にプラグマとして命令を埋め込み、対応コンパイラでコンパイルすることで並列化を行う。現在ではVisual C++やGCCなど、多くのコンパイラがこのOpenMPをサポートしており、またWindowsやLinux、各種UNIXからスーパーコンピュータまで、さまざまなプラットフォームで利用できるのも特徴だ。
OpenMPは1997年に発表された標準規格であり、1997年にFortran向けの「OpenMP fortran 1.0」が、1998年にC/C++向けの「OpenMP C/C++ 1.0」がリリースされた。2009年時点の最新版は「OpenMP 3.0」である。OpenMP規格の策定にはインテルやAMDといったCPUメーカーや富士通、ヒューレット・パッカード、IBM、NEC、サン・マイクロシステムズ、マイクロソフトといったハードウェアベンダーやOSベンダーが参加しており、現在ではインテル コンパイラーやVisual C++、GCC、各社が提供するSolarisやAIXなどの商用UNIX向けコンパイラーなどでなんらかのOpenMPサポートが行われている。
C/C++でOpenMPを利用する場合「omp.h」をインクルードし、並列化を行いたい個所に「#pragma omp ＜ディレクティブ＞」というプラグマを挿入することで、コンパイラに並列化を行う指示を与える。OpenMPには多数のディレクティブが用意されているが、そのなかで多く利用されるであろうものを表3にまとめておいた。そのほかの詳細については別途文献などを参照してほしい。また、OpenMPを使用するプログラムをコンパイルする際は、Visual C++の場合/openmpオプションを、インテル コンパイラの場合/Qopenmpオプションを、GCCの場合-fopenmpオプションを付ける必要がある。
| ディレクティブ | 意味 |
|---|---|
| 並列処理ブロックの生成 | |
| parallel | 並列化するブロックを開始する。このプラグマの直後の「{」から「}」までが並列化を行うブロックとなる |
| for | このプラグマの直後のforループを並列化する。「schedule」オプションでどのようにループを並列化するかも指定可能 |
| sections | このプラグマに続くブロック（{}で囲まれた個所）内の複数の「section」を並列実行する |
| section | 並列実行するブロックを指定する。「sections」と組み合わせて使用する |
| single | 並列処理の際、続くブロックを1つのスレッドだけで実行する |
| 変数の扱い | |
| threadprivate（変数） | それぞれのスレッドで独立して利用したいグローバル変数を指定する |
| private（変数） | それぞれのスレッドで独立して利用したいローカル変数を指定する |
| 同期に関する構文 | |
| master | 続くブロックをマスタースレッドだけで実行する |
| critical | 続くブロックを同時に1つのスレッドだけで実行する |
| barrier | 並列実行しているすべてのスレッドがこのプラグマ部分に到着するまで、このプラグマ部分で待機する |
たとえばOpenMPを利用してループを並列化するには、次のようにすればよい。
#pragma omp parallel
#pragma omp for
for (i = 0; i count; i++) {
/* 並列実行したい処理をここに記述する */
}
このように記述された処理を実行する際、OpenMPライブラリは実行環境で同時に実行できるスレッド数などに応じて適切にスレッドを作成し、処理を各スレッドに振り分けて並列に実行させる。もちろん、ユーザーが利用するスレッド数を明示的に指定することも可能だ。このとき、ループの制御変数「i」は各スレッドで独立したものとなる。
また、下記は複数の処理を同時に実行させる場合の例である。
#pragma omp parallel
#pragma omp sections
{
#pragma omp section
/* 並列に実行したい処理A */
#pragma omp section
/* 並列に実行したい処理B */
#pragma omp section
/* 並列に実行したい処理C */
}
/* すべての処理が完了後に実行する処理 */
この例の場合、処理A、B、Cがそれぞれ別のスレッドで並列に実行され、すべての処理が完了後に#pragma omp parallel以降の処理が実行される。
さて、今度はOpenMPを利用し、先のシングルスレッド版メディアンフィルタプログラムを並列化してみよう。OpenMPを使用したループの並列化は非常に簡単で、ループの前に2行のコードを追加するだけである。
ここで、「#pragma omp parallel private(x, i, j, tmp_array, dx, dy)」という行ではスレッドごとに独立して利用すべき変数を指定しており、また「#pragma omp for」という行では続くforループを並列実行せよ、という指示を与えている。
#pragma omp parallel private(x, i, j, tmp_array, dx, dy)
#pragma omp for
for(y = 1; y height - 1; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3 * x + i] = buf_in-array[y+dy][3 * (x+dx) + i];
}
}
}
このように、OpenMPではスレッドを利用した場合と比べ、少ない変更点で並列化を実装できていることが分かる（表4）。処理時間はマルチスレッドを利用したものと比べてわずかに多いが、無視できる範囲であろう。
| プログラム | 実行時間 |
|---|---|
| シングルスレッド版（med_serial.c） | 5553ミリ秒 |
| OpenMP版（med_omp.c） | 2995ミリ秒 |
OpenMPはC++でも利用できるものの、ループやブロックを対象として並列化を行うため、C++で利用する場合には扱いにくい場合もある。そこで、C++で並列処理を実装する場合に検討したいのがThreading Building Blocks（TBB）である。
TBBはC++のテンプレート機構を活用した並列処理実装のためのテンプレートライブラリで、処理を分割して実行するアルゴリズムや並列処理を行うデータを格納するためのコンテナ、メモリ管理機構、排他・同期処理機構、タスクスケジューラなどから構成されている。TBBはループを置き換えるのではなく、処理もしくはデータを分割して実行するという観点から設計されているのも特徴である。
TBBはインテルが開発し、登場した当初は商用製品としてリリースされていたが、2007年にオープンソース化され（ライセンスはGPL2）、現在は無償で入手および利用が可能だ。また、インテル コンパイラーやParallel Studioなどには非フリーのソフトウェアからでも利用できる商用版が付属している。
TBBを用いて並列処理を実装する基本的な考え方を示したものが図5である。
TBBを利用する場合、まず分割して処理するデータにアクセスするためのイテレータを定義し、そのイテレータと実行する処理となる関数オブジェクトをTBBが用意する並列処理アルゴリズムに与えて分割実行させる、という形になる。データを分割するだけでなく、処理を分割してパイプライン処理させることも可能だ。
データにアクセスするためのイテレータはTBBのテンプレートクラスとしていくつかが用意されているほか、利用するデータに応じてオリジナルのイテレータを実装することもできる。並列処理アルゴリズムとしては、表5のようなものに加え、パイプライン処理を実装するためのクラス「pipeline」が用意されている。
| アルゴリズム（関数）名 | 説明 |
|---|---|
| parallel_for() | 指定された範囲に対して処理を並列に実行する |
| parallel_reduce() | 指定された範囲に対して処理を実行し、それぞれの実行結果に対して和や積、最大/最小といった集計操作を行うリダクション処理を実行する |
| parallel_scan() | 指定された範囲に対して並列スキャン処理を実行する |
| parallel_while() | 実行前にはループ回数が分からないループ処理を並列に実行する |
| parallel_sort() | 並列ソートを実行する |
たとえばTBBを使ってループを並列化する場合、次のようにする。
/* 必要なヘッダーファイルをインポート */
#include "tbb/task_scheduler_init.h"
#include "tbb/parallel_for.h"
#include "tbb/blocked_range.h"
/* tbb名前空間をインポート */
using namespace tbb;
：
：
class ParallelProc {
public:
/*
* operator(): 並列実行する処理を定義する関数オブジェクト
* const blocked_rangesize_t r: 処理する範囲
*/
void operator() (const blocked_rangeint r) const {
int i;
for (i = r.begin(); i != r.end(); i++) {
/* 並列化したい処理をここに記述
}
}
}
void hogehoge() {
/* 1スレッドで実行する処理 /*
：
：
/* 1スレッドで実行する処理ここまで /*
/* TBBのスケジューラを初期化 */
task_scheduler_init init;
/* 並列処理を行うデータ範囲を用意 */
blocked_range2dint range(0, n, 1000)
/* 処理を行う関数オブジェクトを作成 */
ParallelProc proc_obj;
/* データ範囲と関数オブジェクトを与えて並列実行 */
parallel_for(range, proc_obj);
/* 1スレッドで実行する処理 /*
：
：
}
TBBを利用した並列化は一見OpenMPによる並列化よりも複雑に見えるが、実行する処理をクラスとして記述するため汎用性が高いのが特徴である。また、ここでは1次元のループを並列化する例を挙げているが、2次元、3次元のループについてもその構造を保ったままに並列化が可能だ。
それでは、TBBを利用してメディアンフィルタプログラムを並列化してみよう。まず、必要なヘッダファイルをインクルードする。TBBは機能ごとにヘッダファイルも分割されており、必要なものを個別にインクルードする形になる。今回はタスクスケジューラおよび「parallel_for」アルゴリズム、そしてデータ範囲を指定するイテレータは2次元の範囲を指定できる「blocked_range2d」を使用するため、それぞれに対応するヘッダファイルをインクルードする。また、TBBは「tbb」という名前空間を使用するので、そちらを省略して利用できるように指定しておく。
#include "tbb/task_scheduler_init.h" #include "tbb/parallel_for.h" #include "tbb/blocked_range2d.h" /* tbb名前空間をインポート */ using namespace tbb;
次に、並列処理を行う関数を関数オブジェクトとして定義する。ここでは、「ApplyMedFilter」というクラスとして定義している。関数に与える引数はクラスのパブリック変数として定義し、コンストラクタで初期化する形とした。そのほか、処理に必要とする関数もこのクラスのprivateメンバ関数として定義しておく。
class ApplyMedFilter {
public:
const image_buffer* buf_in;
const image_buffer* buf_out;
const image_buffer* buf_tmp;
ApplyMedFilter(const image_buffer* in,
const image_buffer* out,
const image_buffer* tmp) {
buf_in = in;
buf_out = out;
buf_tmp = tmp;
}
private:
/* 使用する作業関数を定義 */
：
：
public:
/*
* operator(): 並列実行する処理
*
* const blocked_range2dsize_t r: 処理する範囲
*/
void operator() (const blocked_range2dint r) const {
int x, y, i, j, dx, dy;
JSAMPLE* tmp_array[9];
for (y = r.rows().begin(); y != r.rows().end(); y++) {
for (x = r.cols().begin(); x != r.cols().end(); x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, x, y);
dx = get_med_position_x(tmp_array, x, y, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3*x + i] = buf_in-array[y+dy][3*(x+dx) + i];
}
}
}
}
};
処理を並列実行させる個所は次のようになる。まずタスクスケジューラを初期化し、関数オブジェクトおよび処理範囲を設定するイテレータを割り当てる。ここで使用しているblocked_range2dは2次元の範囲を指定できるイテレータで、ここではX=1からX=width-1、Y=1からY=height-1までの範囲を指定している。
最後に使用する処理アルゴリズム（ここではparallel_for）に処理を定義した関数オブジェクトと処理する範囲を指定したイテレータを与えて実行すれば良い。
task_scheduler_init init;
ApplyMedFilter filter(buf_in, buf_out, buf_tmp);
/*
* blocked_range2dで処理範囲を指定する。
* 第1引数および第2引数で1次元目の開始値及び終了値を、
* 第4引数及び第5引数で2次元目の開始値および終了値を指定する。
* 第3引数および第6引数では1次元目および2次元目の分割の粒度を指定する
*/
blocked_range2dint range(1, height-1, 10000, 1, width-1, 10000);
parallel_for(range, filter);
なお、このプログラムを実行した際にかかった時間は下記のとおりだ（表6）。
| プログラム | 実行時間 |
|---|---|
| シングルスレッド版（med_serial.c） | 5553ミリ秒 |
| TBB版（med_tbb.cpp） | 3182ミリ秒 |
TBBを利用する場合も、環境に応じてスレッドが生成され、マルチコアCPUでより効率的な処理が可能である。また、並列化の粒度を柔軟に制御できるのも特徴である。
コンパイラによっては、ループを並列に実行するなどの並列化を自動的に行う「自動並列化」機能を備えているものがある。自動並列化を備えたコンパイラとして代表的なのが、Parallel Composerに含まれているほか、単体製品としても販売されている「インテル C++ コンパイラー」だ。
自動並列化を利用することで、ユーザーがコードを変更することなしに、プログラムを並列化することができる。ただし、自動並列化機能ではあらかじめ実行する回数が決まっているループしか並列化できず、OpenMPやTBBを使って人力で並列化を行ったコードと比べるとパフォーマンスはやはり劣ることが多い。とはいえ、何も手間をかけずにコードを並列化できるため、使わない手はない。OpenMPやTBBと併用し、自動並列化機能では並列化できないようなループについてはOpenMPやTBBを利用する、といった使い方が効果的だろう。
また、インテル コンパイラーはループをCPUがより高速に実行できる形に置き換える機能や、CPUが備えるSSEなどの高速演算機能をより活用するようなコードを出力する機能も備えている。これにより、Visual C++でコンパイルするよりも高速なプログラムを作成することが可能だ。
たとえば先で紹介してきたメディアンフィルタプログラムについて、Visual C++でコンパイルしたものとParallel Composerに含まれるインテル C++ コンパイラーでコンパイルしたもの、そしてインテル C++ コンパイラーで自動並列化を有効にしてコンパイルしたものを比較したのが表7である。この結果からも、単純にインテル C++ コンパイラーを利用するだけでパフォーマンスが向上していることが分かる。
| 並列化手法 | 所要時間 | ||
|---|---|---|---|
| Visual C++ 2008 | インテル C++ コンパイラー 11.1（Parallel Composer） | C++ コンパイラー 11.1＋自動並列化 | |
| シングルスレッド版（med_serial.c） | 5598 | 3882 | 3802 |
| マルチスレッド版（med_thread.c） | 2926 | 2064 | 2008 |
| OpenMP版（med_omp.c） | 3033 | 2037 | 2017 |
| TBB版（med_tbb.cpp） | 3186 | 2210 | 2099 |
※単位はすべてミリ秒
以上、マルチスレッド、OpenMP、TBBを使用した並列処理の実装について紹介してきたが、それぞれに特徴があることが理解いただけただろうか。マルチスレッドによる実装はもっとも柔軟ではあるものの、実装の手間やスケーラビリティに欠ける一方、OpenMPによる実装は手軽で、かつスケーラビリティも確保できる。また、TBBはOpenMPよりも柔軟だが、マルチスレッドの場合ほどは実装が複雑にならない、ちょうど中間に位置するものとも言えるだろう。
今回は多数繰り返し実行されるループを並列化する、という形で並列処理を実装しているが、もちろんアプリケーションによってどのように並列化を実装すべきか、というのは異なる。そのため、状況に応じて、適切な並列化技術を使い分けると良いだろう。
/*
* -*- coding: utf-16 -*-
* med_serial.c
*/
#include stdio.h
#include stdlib.h
#include windows.h
#include mmsystem.h
#include "jpeg-6b/jpeglib.h"
/* 画像データを保存するバッファ */
typedef struct {
int width;
int height;
int components;
JSAMPARRAY array;
} image_buffer;
int exit_with_usage()
{
char usage[] = "usage: test input_file output_file\n";
fprintf(stderr, usage);
return -1;
}
image_buffer* alloc_image_buffer(int width, int height, int components)
{
int i, j;
image_buffer* buf;
buf = (image_buffer*)malloc(sizeof(image_buffer));
if (buf == NULL) {
return NULL;
}
buf-width = width;
buf-height = height;
buf-components = components;
buf-array = (JSAMPARRAY)malloc(sizeof(JSAMPROW)*height);
if (buf == NULL) {
free(buf);
return NULL;
}
/* 連続したメモリ領域を割り当てる */
buf-array[0] = (JSAMPROW)calloc(sizeof(JSAMPLE), width*components * height);
if (buf-array[0] == NULL) {
free(buf-array[0]);
free(buf-array);
free(buf);
return NULL;
}
for (i = 1; i height; i++) {
buf-array[i] = buf-array[0] + width*components * i;
}
return buf;
}
void free_image_buffer(image_buffer* buf)
{
int i;
free(buf-array[0]);
free(buf-array);
free(buf);
}
int write_jpg(image_buffer* buf, const char* filename)
{
struct jpeg_compress_struct cinfo;
struct jpeg_error_mgr jerr;
int quality = 90;
FILE* fp;
/* initialize */
cinfo.err = jpeg_std_error(jerr);
jpeg_create_compress(cinfo);
/* open and set output file */
fp = fopen(filename, "wb");
if (fp == NULL) {
fprintf(stderr, "cannot save file: %s.\n", filename);
return -1;
}
jpeg_stdio_dest(cinfo, fp);
/* initialize image attributes */
cinfo.image_width = buf-width;
cinfo.image_height = buf-height;
cinfo.input_components = buf-components;
cinfo.in_color_space = JCS_RGB;
jpeg_set_defaults( cinfo );
/* set quality*/
jpeg_set_quality(cinfo, quality, TRUE);
/* do compression */
jpeg_start_compress( cinfo, TRUE );
jpeg_write_scanlines(cinfo, buf-array, buf-height);
/* finish */
jpeg_finish_compress(cinfo);
jpeg_destroy_compress(cinfo);
fclose(fp);
return 1;
}
image_buffer* read_jpg(const char* filename)
{
/* Allocate and initialize a JPEG decompression object */
struct jpeg_decompress_struct cinfo;
struct jpeg_error_mgr jerr;
FILE* fp;
image_buffer* buf = NULL;
cinfo.err = jpeg_std_error(jerr);
jpeg_create_decompress(cinfo);
/* Specify the source of the compressed data (eg, a file) */
fp = fopen(filename, "rb");
if (fp == NULL) {
fprintf(stderr, "cannot open file: %s.\n", filename);
return NULL;
}
jpeg_stdio_src(cinfo, fp);
/* Call jpeg_read_header() to obtain image info */
jpeg_read_header(cinfo, TRUE);
jpeg_start_decompress(cinfo);
/* allocate buffer */
buf = alloc_image_buffer(cinfo.output_width,
cinfo.output_height,
cinfo.out_color_components);
if (buf == NULL) {
fclose(fp);
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
return NULL;
}
/* copy jpg image's data to memory buffer */
while(cinfo.output_scanline cinfo.output_height) {
jpeg_read_scanlines(cinfo,
buf-array + cinfo.output_scanline,
cinfo.output_height - cinfo.output_scanline);
}
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
fclose(fp);
return buf;
}
void monochrome(image_buffer* buf_in, image_buffer* buf_out)
{
int x, y;
for (y = 0; y buf_in-height; y++) {
for (x = 0; x buf_in-width; x++) {
buf_out-array[y][x] =
0.298912*buf_in-array[y][3*x + 0]
+ 0.586611*buf_in-array[y][3*x + 1]
+ 0.114478*buf_in-array[y][3*x + 2];
}
}
}
int compar(const void* ptr1, const void* ptr2) {
JSAMPLE** p1 = (JSAMPLE**)ptr1;
JSAMPLE** p2 = (JSAMPLE**)ptr2;
return (int)(**p1) - (int)(**p2);
}
/*
* フィルタの探索対象9ピクセルをソートする
*
* JSAMPLE* values[]:
* JSAMPARRAY中の対象のJSAMPLEのアドレスが入った配列
*
* int count: 配列の要素数
*
* 戻り値: 中間値となるJSAMPLEのインデックス
*/
void med_sort(JSAMPLE* values[], int count)
{
qsort(values, count, sizeof(JSAMPLE*), compar);
}
/* ソート結果から中央値を持つピクセルのy位置を計算 */
int get_med_position_y(JSAMPLE* tmp_array[], JSAMPLE* ptr_base) {
int dif;
const JSAMPLE* ptr_med = tmp_array[4];
dif = ptr_med - ptr_base;
if (dif -1) {
return -1;
} else if (1 dif) {
return 1;
} else {
return 0;
}
}
/* ソート結果から中央値を持つピクセルのx位置を計算 */
int get_med_position_x(JSAMPLE* tmp_array[], JSAMPLE* ptr_base, int width, int pos_y) {
const JSAMPLE* ptr_med = tmp_array[4];
return ptr_med - (pos_y * width) - ptr_base;
}
void med_filter(image_buffer* buf_in, image_buffer* buf_out, image_buffer* buf_tmp)
{
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
JSAMPLE* ptr;
const int width = buf_in-width;
const int height = buf_in-height;
fprintf(stderr, "create monochrome image...\n");
/* 作業用バッファにモノクロ化した画像をセット */
monochrome(buf_in, buf_tmp);
/* 簡略化のため、1行目および最終行、1列目および最終列にはフィルタを適用しない */
y = 0; /* 1行目 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
y = height - 1; /* 最終行 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
/* 1列目および最終列 */
for (y = 0; y height; y++) {
for (i = 0; i 3; i++) {
buf_out-array[y][0 + i] = buf_in-array[y][0 + i];
buf_out-array[y][(width - 1) * 3 + i] = buf_in-array[y][(width - 1) * 3 + i];
}
}
fprintf(stderr, "apply filter...\n");
/* フィルタを適用 */
for(y = 1; y height - 1; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3 * x + i] = buf_in-array[y+dy][3 * (x+dx) + i];
}
}
}
/* 終了 */
fprintf(stderr, "done.\n");
}
int main(int argc, char* argv[])
{
image_buffer* buf = NULL;
image_buffer* buf_out = NULL;
image_buffer* buf_tmp = NULL;
DWORD start;
DWORD end;
/* check arguments */
if (argc 3) {
return exit_with_usage();
}
fprintf(stderr, "reading jpeg...\n");
buf = read_jpg(argv[1]);
if (buf == NULL) {
fprintf(stderr, "file read error occured. exit.\n");
return -1;
}
/*
fprintf(stderr, "writing ppm...\n");
if (write_ppm(buf, argv[2]) 0)
return -1;
*/
/* proc */
fprintf(stderr, "apply filter...\n");
buf_out = alloc_image_buffer(buf-width, buf-height, buf-components);
buf_tmp = alloc_image_buffer(buf-width, buf-height, 1);
if ( !buf_out || !buf_tmp) {
free_image_buffer(buf_out);
free_image_buffer(buf_tmp);
free_image_buffer(buf);
return -1;
}
start = timeGetTime();
med_filter(buf, buf_out, buf_tmp);
end = timeGetTime();
/* write */
fprintf(stderr, "writing jpg...\n");
if (write_jpg(buf_out, argv[2]) 0) {
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
return -1;
}
/* finish */
fprintf(stderr, "finishing...\n");
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
fprintf(stderr, "done.\n");
printf("%s, time: %d ms.\n", argv[0], end - start);
return 0;
}
/*
* -*- coding: utf-16 -*-
* test.c
* test code for paralel programming.
* by hylom, 2009.
*/
#include stdio.h
#include stdlib.h
#include process.h
#include windows.h
#include mmsystem.h
#include "jpeg-6b/jpeglib.h"
/*
* usage: test input_file output_file
*/
typedef struct {
int width;
int height;
int components;
JSAMPARRAY array;
} image_buffer;
typedef struct {
image_buffer* buf_in;
image_buffer* buf_out;
image_buffer* buf_tmp;
int y_start;
int y_end;
} thread_param;
/* exit_with_usage: 用法を出力する */
int exit_with_usage()
{
char usage[] = "usage: test input_file output_file\n";
fprintf(stderr, usage);
return -1;
}
/*
* alloc_image_buffer: 作業用の画像バッファを確保
*
* int width: 画像の幅
* int height: 画像の高さ
* int components: モノクロなら1、カラーなら3
*
* 戻り値：確保した画像バッファのポインタ
*/
image_buffer* alloc_image_buffer(int width, int height, int components)
{
int i, j;
image_buffer* buf;
buf = (image_buffer*)malloc(sizeof(image_buffer));
if (buf == NULL) {
return NULL;
}
buf-width = width;
buf-height = height;
buf-components = components;
/*
* jpeglibでは各列のデータを格納する配列（JSAMPROW）の
* 配列として画像を扱う
*/
buf-array = (JSAMPARRAY)malloc(sizeof(JSAMPROW)*height);
if (buf == NULL) {
free(buf);
return NULL;
}
/*
* 各列を独立したメモリ領域として割り当てるのではなく、
* 連続したメモリ領域として割り当てる
*/
buf-array[0] = (JSAMPROW)calloc(sizeof(JSAMPLE), width*components * height);
if (buf-array[0] == NULL) {
free(buf-array[0]);
free(buf-array);
free(buf);
return NULL;
}
/*
* 連続して割り当てたメモリ領域から、各列に相当するアドレスを
* JSAMPARRAYに代入
*/
for (i = 1; i height; i++) {
buf-array[i] = buf-array[0] + width*components * i;
}
return buf;
}
/* free_image_buffer: 割り当てた画像バッファを解放 */
void free_image_buffer(image_buffer* buf)
{
int i;
free(buf-array[0]);
free(buf-array);
free(buf);
}
/*
* write_jpg: 画像バッファの内容をJPEG画像としてファイルに出力
*
* image_bugger* buf: 出力する画像が含まれた画像バッファ
* const char* filename: 出力ファイルのパス名
*
* 戻り値：成功なら1、失敗ならそれ以外
*/
int write_jpg(image_buffer* buf, const char* filename)
{
struct jpeg_compress_struct cinfo;
struct jpeg_error_mgr jerr;
int quality = 90;
FILE* fp;
/* initialize */
cinfo.err = jpeg_std_error(jerr);
jpeg_create_compress(cinfo);
/* open and set output file */
fp = fopen(filename, "wb");
if (fp == NULL) {
fprintf(stderr, "cannot save file: %s.\n", filename);
return -1;
}
jpeg_stdio_dest(cinfo, fp);
/* initialize image attributes */
cinfo.image_width = buf-width;
cinfo.image_height = buf-height;
cinfo.input_components = buf-components;
cinfo.in_color_space = JCS_RGB;
jpeg_set_defaults( cinfo );
/* set quality*/
jpeg_set_quality(cinfo, quality, TRUE);
/* do compression */
jpeg_start_compress( cinfo, TRUE );
jpeg_write_scanlines(cinfo, buf-array, buf-height);
/* finish */
jpeg_finish_compress(cinfo);
jpeg_destroy_compress(cinfo);
fclose(fp);
return 1;
}
/*
* read_jpg: JPEG画像ファイルを読み込み画像バッファに格納
*
* const char* filename: 読み込むファイルのパス名
*
* 戻り値：成功なら画像バッファのポインタ、失敗ならNULL
*
* 注意点：ここで返された画像バッファは呼び出し元の責任において解放すること
*/
image_buffer* read_jpg(const char* filename)
{
/* Allocate and initialize a JPEG decompression object */
struct jpeg_decompress_struct cinfo;
struct jpeg_error_mgr jerr;
FILE* fp;
image_buffer* buf = NULL;
cinfo.err = jpeg_std_error(jerr);
jpeg_create_decompress(cinfo);
/* Specify the source of the compressed data (eg, a file) */
fp = fopen(filename, "rb");
if (fp == NULL) {
fprintf(stderr, "cannot open file: %s.\n", filename);
return NULL;
}
jpeg_stdio_src(cinfo, fp);
/* Call jpeg_read_header() to obtain image info */
jpeg_read_header(cinfo, TRUE);
jpeg_start_decompress(cinfo);
/* allocate buffer */
buf = alloc_image_buffer(cinfo.output_width,
cinfo.output_height,
cinfo.out_color_components);
if (buf == NULL) {
fclose(fp);
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
return NULL;
}
/* copy jpg image's data to memory buffer */
while(cinfo.output_scanline cinfo.output_height) {
jpeg_read_scanlines(cinfo,
buf-array + cinfo.output_scanline,
cinfo.output_height - cinfo.output_scanline);
}
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
fclose(fp);
return buf;
}
/*
* monochrome: 画像をモノクロ化
*
* image_buffer* buf_in: モノクロ化する画像を含む画像バッファ
* image_buffer* buf_out: モノクロ後の画像を含む画像バッファ。
*
* 注意点：buf_outはcomponent=1で割り当てられている必要がある
*/
void monochrome(image_buffer* buf_in, image_buffer* buf_out)
{
int x, y;
for (y = 0; y buf_in-height; y++) {
for (x = 0; x buf_in-width; x++) {
/* RGB各色に重みをかけた平均を取る*/
buf_out-array[y][x] =
0.298912*buf_in-array[y][3*x + 0]
+ 0.586611*buf_in-array[y][3*x + 1]
+ 0.114478*buf_in-array[y][3*x + 2];
}
}
}
/*
* compar: qsortで使用する比較関数
*
* ptr1ポインタが指すポインタの指す変数の値と、
* ptr2ポインタが指すポインタの指す変数の値を比較する
*/
int compar(const void* ptr1, const void* ptr2) {
JSAMPLE** p1 = (JSAMPLE**)ptr1;
JSAMPLE** p2 = (JSAMPLE**)ptr2;
return (int)(**p1) - (int)(**p2);
}
/*
* med_sort:フィルタの探索対象9ピクセルをソートする
*
* JSAMPLE* values[]:
* JSAMPARRAY中の対象のJSAMPLEのアドレスが入った配列
*
* int count: 配列の要素数
*/
void med_sort(JSAMPLE* values[], int count)
{
qsort(values, count, sizeof(JSAMPLE*), compar);
}
/* ソート結果から中央値を持つピクセルのy位置を計算 */
int get_med_position_y(JSAMPLE* tmp_array[], JSAMPLE* ptr_base) {
int dif;
const JSAMPLE* ptr_med = tmp_array[4];
dif = ptr_med - ptr_base;
if (dif -1) {
return -1;
} else if (1 dif) {
return 1;
} else {
return 0;
}
}
/* ソート結果から中央値を持つピクセルのx位置を計算 */
int get_med_position_x(JSAMPLE* tmp_array[], JSAMPLE* ptr_base, int width, int pos_y) {
const JSAMPLE* ptr_med = tmp_array[4];
return ptr_med - (pos_y * width) - ptr_base;
}
void __stdcall work_thread(thread_param* param) {
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
const image_buffer* buf_in = param-buf_in;
const image_buffer* buf_out = param-buf_out;
const image_buffer* buf_tmp = param-buf_tmp;
const int y_start = param-y_start;
const int y_end = param-y_end;
const int width = param-buf_in-width;
/* フィルタを適用 */
for(y = y_start; y y_end; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3*x + i] = buf_in-array[y+dy][3*(x+dx) + i];
}
}
}
}
void med_filter(image_buffer* buf_in, image_buffer* buf_out, image_buffer* buf_tmp)
{
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
const int width = buf_in-width;
const int height = buf_in-height;
thread_param param[2];
HANDLE thread;
fprintf(stderr, "create monochrome image...\n");
/* 作業用バッファにモノクロ化した画像をセット */
monochrome(buf_in, buf_tmp);
/* 簡略化のため、1行目および最終行、1列目および最終列にはフィルタを適用しない */
y = 0; /* 1行目 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
y = height - 1; /* 最終行 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
/* 1列目および最終列 */
for (y = 0; y height; y++) {
for (i = 0; i 3; i++) {
buf_out-array[y][0 + i] = buf_in-array[y][0 + i];
buf_out-array[y][3*(width - 1) + i] = buf_in-array[y][3*(width - 1) + i];
}
}
fprintf(stderr, "apply filter...\n");
/* フィルタを適用 */
for (i = 0; i 2; i++) {
param[i].buf_in = buf_in;
param[i].buf_out = buf_out;
param[i].buf_tmp = buf_tmp;
}
param[0].y_start = 1;
param[0].y_end = (height - 1) / 2;
param[1].y_start = (height - 1) / 2 + 1;
param[1].y_end = height - 1;
/* スレッド生成 */
thread = (HANDLE)_beginthreadex(NULL, 0, (void*)work_thread, (param[0]), 0, NULL);
work_thread((param[1]));
/* スレッドの終了を待つ */
WaitForSingleObject(thread, INFINITE);
CloseHandle(thread);
/* 終了 */
fprintf(stderr, "done.\n");
}
int main(int argc, char* argv[])
{
image_buffer* buf = NULL;
image_buffer* buf_out = NULL;
image_buffer* buf_tmp = NULL;
DWORD start;
DWORD end;
/* check arguments */
if (argc 3) {
return exit_with_usage();
}
fprintf(stderr, "reading jpeg...\n");
buf = read_jpg(argv[1]);
if (buf == NULL) {
fprintf(stderr, "file read error occured. exit.\n");
return -1;
}
/*
fprintf(stderr, "writing ppm...\n");
if (write_ppm(buf, argv[2]) 0)
return -1;
*/
/* proc */
fprintf(stderr, "apply filter...\n");
buf_out = alloc_image_buffer(buf-width, buf-height, buf-components);
buf_tmp = alloc_image_buffer(buf-width, buf-height, 1);
if ( !buf_out || !buf_tmp) {
free_image_buffer(buf_out);
free_image_buffer(buf_tmp);
free_image_buffer(buf);
return -1;
}
start = timeGetTime();
med_filter(buf, buf_out, buf_tmp);
end = timeGetTime();
/* write */
fprintf(stderr, "writing jpg...\n");
if (write_jpg(buf_out, argv[2]) 0) {
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
return -1;
}
/* finish */
fprintf(stderr, "finishing...\n");
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
fprintf(stderr, "done.\n");
printf("%s, time: %d ms.\n", argv[0], end - start);
return 0;
}
/*
* -*- coding: utf-16 -*-
* test.c
* test code for paralel programming.
* by hylom, 2009.
*/
#include stdio.h
#include stdlib.h
#include windows.h
#include mmsystem.h
#include "jpeg-6b/jpeglib.h"
/*
* usage: test input_file output_file
*/
typedef struct {
int width;
int height;
int components;
JSAMPARRAY array;
} image_buffer;
int exit_with_usage()
{
char usage[] = "usage: test input_file output_file\n";
fprintf(stderr, usage);
return -1;
}
image_buffer* alloc_image_buffer(int width, int height, int components)
{
int i, j;
image_buffer* buf;
buf = (image_buffer*)malloc(sizeof(image_buffer));
if (buf == NULL) {
return NULL;
}
buf-width = width;
buf-height = height;
buf-components = components;
buf-array = (JSAMPARRAY)malloc(sizeof(JSAMPROW)*height);
if (buf == NULL) {
free(buf);
return NULL;
}
/* 連続したメモリ領域を割り当てる */
buf-array[0] = (JSAMPROW)calloc(sizeof(JSAMPLE), width*components * height);
if (buf-array[0] == NULL) {
free(buf-array[0]);
free(buf-array);
free(buf);
return NULL;
}
for (i = 1; i height; i++) {
buf-array[i] = buf-array[0] + width*components * i;
}
return buf;
}
void free_image_buffer(image_buffer* buf)
{
int i;
free(buf-array[0]);
free(buf-array);
free(buf);
}
int write_jpg(image_buffer* buf, const char* filename)
{
struct jpeg_compress_struct cinfo;
struct jpeg_error_mgr jerr;
int quality = 90;
FILE* fp;
/* initialize */
cinfo.err = jpeg_std_error(jerr);
jpeg_create_compress(cinfo);
/* open and set output file */
fp = fopen(filename, "wb");
if (fp == NULL) {
fprintf(stderr, "cannot save file: %s.\n", filename);
return -1;
}
jpeg_stdio_dest(cinfo, fp);
/* initialize image attributes */
cinfo.image_width = buf-width;
cinfo.image_height = buf-height;
cinfo.input_components = buf-components;
cinfo.in_color_space = JCS_RGB;
jpeg_set_defaults( cinfo );
/* set quality*/
jpeg_set_quality(cinfo, quality, TRUE);
/* do compression */
jpeg_start_compress( cinfo, TRUE );
jpeg_write_scanlines(cinfo, buf-array, buf-height);
/* finish */
jpeg_finish_compress(cinfo);
jpeg_destroy_compress(cinfo);
fclose(fp);
return 1;
}
image_buffer* read_jpg(const char* filename)
{
/* Allocate and initialize a JPEG decompression object */
struct jpeg_decompress_struct cinfo;
struct jpeg_error_mgr jerr;
FILE* fp;
image_buffer* buf = NULL;
cinfo.err = jpeg_std_error(jerr);
jpeg_create_decompress(cinfo);
/* Specify the source of the compressed data (eg, a file) */
fp = fopen(filename, "rb");
if (fp == NULL) {
fprintf(stderr, "cannot open file: %s.\n", filename);
return NULL;
}
jpeg_stdio_src(cinfo, fp);
/* Call jpeg_read_header() to obtain image info */
jpeg_read_header(cinfo, TRUE);
jpeg_start_decompress(cinfo);
/* allocate buffer */
buf = alloc_image_buffer(cinfo.output_width,
cinfo.output_height,
cinfo.out_color_components);
if (buf == NULL) {
fclose(fp);
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
return NULL;
}
/* copy jpg image's data to memory buffer */
while(cinfo.output_scanline cinfo.output_height) {
jpeg_read_scanlines(cinfo,
buf-array + cinfo.output_scanline,
cinfo.output_height - cinfo.output_scanline);
}
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
fclose(fp);
return buf;
}
void monochrome(image_buffer* buf_in, image_buffer* buf_out)
{
int x, y;
for (y = 0; y buf_in-height; y++) {
for (x = 0; x buf_in-width; x++) {
buf_out-array[y][x] =
0.298912*buf_in-array[y][3*x + 0]
+ 0.586611*buf_in-array[y][3*x + 1]
+ 0.114478*buf_in-array[y][3*x + 2];
}
}
}
int compar(const void* ptr1, const void* ptr2) {
JSAMPLE** p1 = (JSAMPLE**)ptr1;
JSAMPLE** p2 = (JSAMPLE**)ptr2;
return (int)(**p1) - (int)(**p2);
}
/*
* フィルタの探索対象9ピクセルをソートする
*
* JSAMPLE* values[]:
* JSAMPARRAY中の対象のJSAMPLEのアドレスが入った配列
*
* int count: 配列の要素数
*
* 戻り値: 中間値となるJSAMPLEのインデックス
*/
void med_sort(JSAMPLE* values[], int count)
{
qsort(values, count, sizeof(JSAMPLE*), compar);
}
/* ソート結果から中央値を持つピクセルのy位置を計算 */
int get_med_position_y(JSAMPLE* tmp_array[], JSAMPLE* ptr_base) {
int dif;
const JSAMPLE* ptr_med = tmp_array[4];
dif = ptr_med - ptr_base;
if (dif -1) {
return -1;
} else if (1 dif) {
return 1;
} else {
return 0;
}
}
/* ソート結果から中央値を持つピクセルのx位置を計算 */
int get_med_position_x(JSAMPLE* tmp_array[], JSAMPLE* ptr_base, int width, int pos_y) {
const JSAMPLE* ptr_med = tmp_array[4];
return ptr_med - (pos_y * width) - ptr_base;
}
void med_filter(image_buffer* buf_in, image_buffer* buf_out, image_buffer* buf_tmp)
{
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
JSAMPLE* ptr;
const int width = buf_in-width;
const int height = buf_in-height;
fprintf(stderr, "create monochrome image...\n");
/* 作業用バッファにモノクロ化した画像をセット */
monochrome(buf_in, buf_tmp);
/* 簡略化のため、1行目および最終行、1列目および最終列にはフィルタを適用しない */
y = 0; /* 1行目 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
y = height - 1; /* 最終行 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
/* 1列目および最終列 */
for (y = 0; y height; y++) {
for (i = 0; i 3; i++) {
buf_out-array[y][0 + i] = buf_in-array[y][0 + i];
buf_out-array[y][(width - 1) * 3 + i] = buf_in-array[y][(width - 1) * 3 + i];
}
}
fprintf(stderr, "apply filter...\n");
/* フィルタを適用 */
#pragma omp parallel private(x, i, j, tmp_array, dx, dy)
#pragma omp for
for(y = 1; y height - 1; y++) {
for (x = 1; x width - 1; x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, (buf_tmp-array[y][x]));
dx = get_med_position_x(tmp_array, (buf_tmp-array[y][x]), width, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3 * x + i] = buf_in-array[y+dy][3 * (x+dx) + i];
}
}
}
/* 終了 */
fprintf(stderr, "done.\n");
}
int main(int argc, char* argv[])
{
image_buffer* buf = NULL;
image_buffer* buf_out = NULL;
image_buffer* buf_tmp = NULL;
DWORD start, end;
/* check arguments */
if (argc 3) {
return exit_with_usage();
}
fprintf(stderr, "reading jpeg...\n");
buf = read_jpg(argv[1]);
if (buf == NULL) {
fprintf(stderr, "file read error occured. exit.\n");
return -1;
}
/*
fprintf(stderr, "writing ppm...\n");
if (write_ppm(buf, argv[2]) 0)
return -1;
*/
/* proc */
fprintf(stderr, "apply filter...\n");
buf_out = alloc_image_buffer(buf-width, buf-height, buf-components);
buf_tmp = alloc_image_buffer(buf-width, buf-height, 1);
if ( !buf_out || !buf_tmp) {
free_image_buffer(buf_out);
free_image_buffer(buf_tmp);
free_image_buffer(buf);
return -1;
}
start = timeGetTime();
med_filter(buf, buf_out, buf_tmp);
end = timeGetTime();
/* write */
fprintf(stderr, "writing jpg...\n");
if (write_jpg(buf_out, argv[2]) 0) {
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
return -1;
}
/* finish */
fprintf(stderr, "finishing...\n");
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
fprintf(stderr, "done.\n");
printf("%s, time: %d ms.\n", argv[0], end - start);
return 0;
}
/*
* -*- coding: utf-16 -*-
* test.c
* test code for paralel programming.
* by hylom, 2009.
*/
#include stdio.h
#include stdlib.h
#include memory.h
#include windows.h
#include mmsystem.h
extern "C" {
#include "jpeg-6b/jpeglib.h"
}
#include "tbb/task_scheduler_init.h"
#include "tbb/parallel_for.h"
#include "tbb/blocked_range2d.h"
/* tbb名前空間をインポート */
using namespace tbb;
/*
* usage: test input_file output_file
*/
typedef struct {
int width;
int height;
int components;
JSAMPARRAY array;
} image_buffer;
int exit_with_usage()
{
char usage[] = "usage: test input_file output_file\n";
fprintf(stderr, usage);
return -1;
}
image_buffer* alloc_image_buffer(int width, int height, int components)
{
int i, j;
image_buffer* buf;
buf = (image_buffer*)malloc(sizeof(image_buffer));
if (buf == NULL) {
return NULL;
}
buf-width = width;
buf-height = height;
buf-components = components;
buf-array = (JSAMPARRAY)malloc(sizeof(JSAMPROW)*height);
if (buf == NULL) {
free(buf);
return NULL;
}
/* 連続したメモリ領域を割り当てる */
buf-array[0] = (JSAMPROW)calloc(sizeof(JSAMPLE), width*components * height);
if (buf-array[0] == NULL) {
free(buf-array[0]);
free(buf-array);
free(buf);
return NULL;
}
for (i = 1; i height; i++) {
buf-array[i] = buf-array[0] + width*components * i;
}
return buf;
}
void free_image_buffer(image_buffer* buf)
{
int i;
free(buf-array[0]);
free(buf-array);
free(buf);
}
int write_jpg(image_buffer* buf, const char* filename)
{
struct jpeg_compress_struct cinfo;
struct jpeg_error_mgr jerr;
int quality = 90;
FILE* fp;
/* initialize */
cinfo.err = jpeg_std_error(jerr);
jpeg_create_compress(cinfo);
/* open and set output file */
fp = fopen(filename, "wb");
if (fp == NULL) {
fprintf(stderr, "cannot save file: %s.\n", filename);
return -1;
}
jpeg_stdio_dest(cinfo, fp);
/* initialize image attributes */
cinfo.image_width = buf-width;
cinfo.image_height = buf-height;
cinfo.input_components = buf-components;
cinfo.in_color_space = JCS_RGB;
jpeg_set_defaults( cinfo );
/* set quality*/
jpeg_set_quality(cinfo, quality, TRUE);
/* do compression */
jpeg_start_compress( cinfo, TRUE );
jpeg_write_scanlines(cinfo, buf-array, buf-height);
/* finish */
jpeg_finish_compress(cinfo);
jpeg_destroy_compress(cinfo);
fclose(fp);
return 1;
}
image_buffer* read_jpg(const char* filename)
{
/* Allocate and initialize a JPEG decompression object */
struct jpeg_decompress_struct cinfo;
struct jpeg_error_mgr jerr;
FILE* fp;
image_buffer* buf = NULL;
cinfo.err = jpeg_std_error(jerr);
jpeg_create_decompress(cinfo);
/* Specify the source of the compressed data (eg, a file) */
fp = fopen(filename, "rb");
if (fp == NULL) {
fprintf(stderr, "cannot open file: %s.\n", filename);
return NULL;
}
jpeg_stdio_src(cinfo, fp);
/* Call jpeg_read_header() to obtain image info */
jpeg_read_header(cinfo, TRUE);
jpeg_start_decompress(cinfo);
/* allocate buffer */
buf = alloc_image_buffer(cinfo.output_width,
cinfo.output_height,
cinfo.out_color_components);
if (buf == NULL) {
fclose(fp);
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
return NULL;
}
/* copy jpg image's data to memory buffer */
while(cinfo.output_scanline cinfo.output_height) {
jpeg_read_scanlines(cinfo,
buf-array + cinfo.output_scanline,
cinfo.output_height - cinfo.output_scanline);
}
jpeg_finish_decompress(cinfo);
jpeg_destroy_decompress(cinfo);
fclose(fp);
return buf;
}
void monochrome(image_buffer* buf_in, image_buffer* buf_out)
{
int x, y;
for (y = 0; y buf_in-height; y++) {
for (x = 0; x buf_in-width; x++) {
buf_out-array[y][x] =
0.298912*buf_in-array[y][3*x + 0]
+ 0.586611*buf_in-array[y][3*x + 1]
+ 0.114478*buf_in-array[y][3*x + 2];
}
}
}
/*
* フィルタの探索対象9ピクセルをソートする
*
* JSAMPLE* values[]:
* JSAMPARRAY中の対象のJSAMPLEのアドレスが入った配列
*
* int count: 配列の要素数
*
* 戻り値: 中間値となるJSAMPLEのインデックス
*/
int compar(const void* ptr1, const void* ptr2) {
JSAMPLE** p1 = (JSAMPLE**)ptr1;
JSAMPLE** p2 = (JSAMPLE**)ptr2;
return (int)(**p1) - (int)(**p2);
}
class ApplyMedFilter {
public:
const image_buffer* buf_in;
const image_buffer* buf_out;
const image_buffer* buf_tmp;
ApplyMedFilter(const image_buffer* in,
const image_buffer* out,
const image_buffer* tmp) {
buf_in = in;
buf_out = out;
buf_tmp = tmp;
}
private:
/* ソート結果から中央値を持つピクセルのy位置を計算 */
int get_med_position_y(JSAMPLE* tmp_array[], int x, int y) const {
int dif;
const JSAMPLE* ptr_med = tmp_array[4];
const JSAMPLE* ptr_base = (buf_tmp-array[y][x]);
dif = ptr_med - ptr_base;
if (dif -1) {
return -1;
} else if (1 dif) {
return 1;
} else {
return 0;
}
}
/* ソート結果から中央値を持つピクセルのx位置を計算 */
int get_med_position_x(JSAMPLE* tmp_array[], int x, int y, int pos_y) const {
const JSAMPLE* ptr_med = tmp_array[4];
const JSAMPLE* ptr_base = (buf_tmp-array[y][x]);
return ptr_med - (pos_y * buf_in-width) - ptr_base;
}
/* ソートを実行 */
void med_sort(JSAMPLE* values[], int count) const {
qsort(values, count, sizeof(JSAMPLE*), compar);
}
public:
/*
* operator(): 並列実行する処理
*
* const blocked_range2dsize_t r: 処理する範囲
*/
void operator() (const blocked_range2dint r) const {
int x, y, i, j, dx, dy;
JSAMPLE* tmp_array[9];
for (y = r.rows().begin(); y != r.rows().end(); y++) {
for (x = r.cols().begin(); x != r.cols().end(); x++) {
/* 対象ピクセルとその近傍8ピクセルのアドレスをソートバッファにコピー*/
for (j = 0; j 3; j++) {
for (i = 0; i 3; i++) {
tmp_array[3*j + i] = (buf_tmp-array[y-1+j][x-1+i]);
}
}
/* 9ピクセルをソート */
med_sort(tmp_array, 9);
/* 中央値となるピクセルの相対位置を取得 */
dy = get_med_position_y(tmp_array, x, y);
dx = get_med_position_x(tmp_array, x, y, dy);
/* 出力バッファの対象ピクセルの値を中央値となるピクセルの値に設定 */
for (i = 0; i 3; i++) {
buf_out-array[y][3*x + i] = buf_in-array[y+dy][3*(x+dx) + i];
}
}
}
}
};
void med_filter(image_buffer* buf_in, image_buffer* buf_out, image_buffer* buf_tmp)
{
int x, y;
int dx, dy;
int i, j;
JSAMPLE* tmp_array[9];
JSAMPLE* ptr;
const int width = buf_in-width;
const int height = buf_in-height;
fprintf(stderr, "create monochrome image...\n");
/* 作業用バッファにモノクロ化した画像をセット */
monochrome(buf_in, buf_tmp);
/* 簡略化のため、1行目および最終行、1列目および最終列にはフィルタを適用しない */
y = 0; /* 1行目 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
y = height - 1; /* 最終行 */
memcpy(buf_out-array[y], buf_in-array[y], sizeof(JSAMPLE)*3*width);
/* 1列目および最終列 */
for (y = 0; y height; y++) {
for (i = 0; i 3; i++) {
buf_out-array[y][0 + i] = buf_in-array[y][0 + i];
buf_out-array[y][(width - 1) * 3 + i] = buf_in-array[y][(width - 1) * 3 + i];
}
}
fprintf(stderr, "apply filter...\n");
/* フィルタを適用 */
task_scheduler_init init;
ApplyMedFilter filter(buf_in, buf_out, buf_tmp);
blocked_range2dint range(1, height-1, 10000, 1, width-1, 10000);
parallel_for(range, filter);
/* 終了 */
fprintf(stderr, "done.\n");
}
int main(int argc, char* argv[])
{
image_buffer* buf = NULL;
image_buffer* buf_out = NULL;
image_buffer* buf_tmp = NULL;
DWORD start, end;
/* check arguments */
if (argc 3) {
return exit_with_usage();
}
fprintf(stderr, "reading jpeg...\n");
buf = read_jpg(argv[1]);
if (buf == NULL) {
fprintf(stderr, "file read error occured. exit.\n");
return -1;
}
/*
fprintf(stderr, "writing ppm...\n");
if (write_ppm(buf, argv[2]) 0)
return -1;
*/
/* proc */
fprintf(stderr, "apply filter...\n");
buf_out = alloc_image_buffer(buf-width, buf-height, buf-components);
buf_tmp = alloc_image_buffer(buf-width, buf-height, 1);
if ( !buf_out || !buf_tmp) {
free_image_buffer(buf_out);
free_image_buffer(buf_tmp);
free_image_buffer(buf);
return -1;
}
start = timeGetTime();
med_filter(buf, buf_out, buf_tmp);
end = timeGetTime();
/* write */
fprintf(stderr, "writing jpg...\n");
if (write_jpg(buf_out, argv[2]) 0) {
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
return -1;
}
/* finish */
fprintf(stderr, "finishing...\n");
free_image_buffer(buf);
free_image_buffer(buf_tmp);
free_image_buffer(buf_out);
fprintf(stderr, "done.\n");
printf("%s, time: %d ms.\n", argv[0], end - start);
return 0;
}
[PageInfo]
LastUpdate: 2009-11-16 19:38:11, ModifiedBy: hiromichi-m
[Permissions]
view:all, edit:login users, delete/config:members