parfor を使用したイメージ コントラストを改善するアルゴリズムの高速化
この例では、イメージのコントラストを改善するために、イメージに簡単なヒストグラム均等化関数を適用する MATLAB® コードからスタンドアロンの C ライブラリを生成する方法を示します。例では、3 つの標準 RGB イメージ平面をそれぞれ個別のスレッドで処理するために parfor
を使用しています。また、MATLAB コードがコードの生成に適していることを確認するために、C コードの生成前に MEX 関数を生成して実行する方法も示します。
MATLAB Coder™ は OpenMP (移植可能な共有メモリ並列プログラミング標準) を使用して parfor
のサポートを実装します。The OpenMP API Specification for Parallel Programming を参照してください。MATLAB は複数のワーカー セッションを作成して parfor
をサポートしますが、MATLAB Coder は OpenMP を使用して同一マシンで実行される複数のスレッドを作成します。
必要条件
並列処理をサポートするには、コンパイラが OpenMP の共有メモリ並列プログラミング標準をサポートしていなければなりません。コンパイラでのサポートがない場合でもこの例は実行できますが、生成されたコードは逐次実行となります。
関数 histequalize
について
関数 histequalize.m
はイメージ (N x M x 3 の行列として表現) を受け取り、コントラストが強調されたイメージを返します。
type histequalize
function equalizedImage = histequalize(originalImage) %#codegen % equalizedImage = histequalize(originalImage) % Histogram equalization (or linearization) for improving image contrast. % Given an NxMx3 image, equalizes the histogram of each of the three image % planes in order to improve image contrast. assert(size(originalImage,1) <= 8192); assert(size(originalImage,2) <= 8192); assert(size(originalImage,3) == 3); assert(isa(originalImage, 'uint8')); [L, originalHist] = computeHistogram(originalImage); equalizedImage = equalize(L, originalHist, originalImage); end function [L, originalHist] = computeHistogram(originalImage) L = double(max(max(max(originalImage)))) + 1; originalHist = coder.nullcopy(zeros(3,L)); sz = size(originalImage); N = sz(1); M = sz(2); parfor plane = 1:sz(3) planeHist = zeros(1,L); for y = 1:N for x = 1:M r = originalImage(y,x,plane); planeHist(r+1) = planeHist(r+1) + 1; end end originalHist(plane,:) = planeHist; end end function equalizedImage = equalize(L, originalHist, originalImage) equalizedImage = coder.nullcopy(originalImage); sz = size(originalImage); N = sz(1); M = sz(2); normalizer = (L - 1)/(N*M); parfor plane = 1:sz(3) planeHist = originalHist(plane,:); for y = 1:N for x = 1:M r = originalImage(y,x,plane); s = 0; for j = 0:int32(r) s = s + planeHist(j+1); end s = normalizer * s; equalizedImage(y,x,plane) = s; end end end end
MEX 関数の生成
codegen
コマンドを使用して MEX 関数を生成します。
codegen histequalize
Code generation successful.
C コードを生成する前に、MATLAB で MEX 関数をテストして、その関数が元の MATLAB コードと機能的に等価であることと実行時のエラーが発生しないことを確認しなければなりません。既定で、codegen
は、現在のフォルダーに histequalize_mex
という名前の MEX 関数を生成します。これにより、MATLAB コードと MEX 関数をテストして結果を比較することができます。
元のイメージの読み取り
標準の imread
コマンドを使用してコントラストの低いイメージを読み取ります。
lcIm = imread('LowContrast.jpg');
image(lcIm);
MEX 関数の実行 (ヒストグラム均等化アルゴリズム)
コントラストの低いイメージを渡します。
hcIm = histequalize_mex(lcIm);
結果の表示
image(hcIm);
スタンドアロン C コードの生成
codegen -config:lib histequalize
Code generation successful.
-config:lib
オプションを指定して codegen
を使用すると、スタンドアロン C ライブラリが生成されます。既定では、ライブラリ用に生成されたコードは codegen/lib/histequalize/
フォルダーにあります。
生成された関数の確認
生成されたコードには、複数のスレッドを使用してコードの並列処理を制御する OpenMP プラグマが含まれていることに注意してください。
type codegen/lib/histequalize/histequalize.c
/* * Prerelease License - for engineering feedback and testing purposes * only. Not for sale. * File: histequalize.c * * MATLAB Coder version : 24.1 * C/C++ source code generated on : 25-Jan-2024 15:07:00 */ /* Include Files */ #include "histequalize.h" #include "histequalize_data.h" #include "histequalize_emxutil.h" #include "histequalize_initialize.h" #include "histequalize_types.h" #include "omp.h" #include <math.h> #include <string.h> /* Function Declarations */ static double computeHistogram(const emxArray_uint8_T *originalImage, double originalHist_data[], int originalHist_size[2]); static void equalize(double L, const double originalHist_data[], const emxArray_uint8_T *originalImage, emxArray_uint8_T *equalizedImage); static double rt_roundd_snf(double u); /* Function Definitions */ /* * Arguments : const emxArray_uint8_T *originalImage * double originalHist_data[] * int originalHist_size[2] * Return Type : double */ static double computeHistogram(const emxArray_uint8_T *originalImage, double originalHist_data[], int originalHist_size[2]) { double planeHist_data[256]; double L; int M; int b_i; int i; int loop_ub; int maxval_idx_1; int p; int plane; int vlen; int x; int xOffset; int xPageOffset; int y; unsigned char maxval_data[24576]; unsigned char maxval[3]; const unsigned char *originalImage_data; unsigned char b_maxval; unsigned char r; originalImage_data = originalImage->data; maxval_idx_1 = originalImage->size[1]; if (originalImage->size[1] == 0) { vlen = originalImage->size[1] * 3; if (vlen - 1 >= 0) { memset(&maxval_data[0], 0, (unsigned int)vlen * sizeof(unsigned char)); } } else { vlen = originalImage->size[0]; M = (unsigned short)(originalImage->size[1] * 3); for (p = 0; p < M; p++) { xPageOffset = p * vlen; maxval_data[p] = originalImage_data[xPageOffset]; for (i = 2; i <= vlen; i++) { xOffset = (xPageOffset + i) - 1; if (maxval_data[p] < originalImage_data[xOffset]) { maxval_data[p] = originalImage_data[xOffset]; } } } } for (p = 0; p < 3; p++) { xPageOffset = p * maxval_idx_1; maxval[p] = maxval_data[xPageOffset]; for (i = 2; i <= maxval_idx_1; i++) { xOffset = (xPageOffset + i) - 1; if (maxval[p] < maxval_data[xOffset]) { maxval[p] = maxval_data[xOffset]; } } } b_maxval = maxval[0]; if (maxval[0] < maxval[1]) { b_maxval = maxval[1]; } if (b_maxval < maxval[2]) { b_maxval = maxval[2]; } L = (double)b_maxval + 1.0; originalHist_size[0] = 3; originalHist_size[1] = b_maxval + 1; vlen = originalImage->size[0]; M = originalImage->size[1]; #pragma omp parallel for num_threads(omp_get_max_threads()) private( \ r, planeHist_data, loop_ub, y, x, b_i) for (plane = 0; plane < 3; plane++) { loop_ub = (int)L; memset(&planeHist_data[0], 0, (unsigned int)loop_ub * sizeof(double)); for (y = 0; y < vlen; y++) { for (x = 0; x < M; x++) { r = originalImage_data[(y + originalImage->size[0] * x) + originalImage->size[0] * originalImage->size[1] * plane]; b_i = (int)(r + 1U); if (r + 1U > 255U) { b_i = 255; } loop_ub = (int)(r + 1U); if (r + 1U > 255U) { loop_ub = 255; } planeHist_data[b_i - 1] = planeHist_data[loop_ub - 1] + 1.0; } } loop_ub = originalHist_size[1]; for (b_i = 0; b_i < loop_ub; b_i++) { originalHist_data[plane + 3 * b_i] = planeHist_data[b_i]; } } return L; } /* * Arguments : double L * const double originalHist_data[] * const emxArray_uint8_T *originalImage * emxArray_uint8_T *equalizedImage * Return Type : void */ static void equalize(double L, const double originalHist_data[], const emxArray_uint8_T *originalImage, emxArray_uint8_T *equalizedImage) { double normalizer; double s; int M; int N; int i; int j; int plane; int x; int y; const unsigned char *originalImage_data; unsigned char r; unsigned char *equalizedImage_data; originalImage_data = originalImage->data; N = equalizedImage->size[0] * equalizedImage->size[1] * equalizedImage->size[2]; equalizedImage->size[0] = originalImage->size[0]; equalizedImage->size[1] = originalImage->size[1]; equalizedImage->size[2] = 3; emxEnsureCapacity_uint8_T(equalizedImage, N); equalizedImage_data = equalizedImage->data; N = originalImage->size[0]; M = originalImage->size[1]; normalizer = (L - 1.0) / ((double)originalImage->size[0] * (double)originalImage->size[1]); #pragma omp parallel for num_threads(omp_get_max_threads()) private( \ s, r, y, x, i, j) for (plane = 0; plane < 3; plane++) { for (y = 0; y < N; y++) { for (x = 0; x < M; x++) { r = originalImage_data[(y + originalImage->size[0] * x) + originalImage->size[0] * originalImage->size[1] * plane]; s = 0.0; i = r; for (j = 0; j <= i; j++) { s += originalHist_data[plane + 3 * j]; } s *= normalizer; s = rt_roundd_snf(s); if (s < 256.0) { if (s >= 0.0) { r = (unsigned char)s; } else { r = 0U; } } else if (s >= 256.0) { r = MAX_uint8_T; } else { r = 0U; } equalizedImage_data[(y + equalizedImage->size[0] * x) + equalizedImage->size[0] * equalizedImage->size[1] * plane] = r; } } } } /* * Arguments : double u * Return Type : double */ static double rt_roundd_snf(double u) { double y; if (fabs(u) < 4.503599627370496E+15) { if (u >= 0.5) { y = floor(u + 0.5); } else if (u > -0.5) { y = u * 0.0; } else { y = ceil(u - 0.5); } } else { y = u; } return y; } /* * equalizedImage = histequalize(originalImage) * Histogram equalization (or linearization) for improving image contrast. * Given an NxMx3 image, equalizes the histogram of each of the three image * planes in order to improve image contrast. * * Arguments : const emxArray_uint8_T *originalImage * emxArray_uint8_T *equalizedImage * Return Type : void */ void histequalize(const emxArray_uint8_T *originalImage, emxArray_uint8_T *equalizedImage) { double originalHist_data[768]; double L; int originalHist_size[2]; if (!isInitialized_histequalize) { histequalize_initialize(); } L = computeHistogram(originalImage, originalHist_data, originalHist_size); equalize(L, originalHist_data, originalImage, equalizedImage); } /* * File trailer for histequalize.c * * [EOF] */