HPC/並列プログラミングポータルでは、HPC(High Performance Computing)プログラミングや並列プログラミングに関する情報を集積・発信しています。

新着トピックス

リスト1 並列化されていないメディアンフィルタプログラム
/*
 * -*- 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;
}

リスト2 スレッドを利用して並列化したメディアンフィルタプログラム
/*
 * -*- 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;
}





リスト3 OpenMPを利用して並列化したメディアンフィルタプログラム
/*
 * -*- 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;
}



リスト4 TBBを利用して並列化したメディアンフィルタプログラム
/*
 * -*- 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;
}