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);

Figure contains an axes object. The axes object contains an object of type image.

MEX 関数の実行 (ヒストグラム均等化アルゴリズム)

コントラストの低いイメージを渡します。

hcIm = histequalize_mex(lcIm);

結果の表示

image(hcIm);

Figure contains an axes object. The axes object contains an object of type image.

スタンドアロン C コードの生成

codegen -config:lib histequalize
Code generation successful.

-config:lib オプションを指定して codegen を使用すると、スタンドアロン C ライブラリが生成されます。既定では、ライブラリ用に生成されたコードは codegen/lib/histequalize/ フォルダーにあります。

生成された関数の確認

生成されたコードには、複数のスレッドを使用してコードの並列処理を制御する OpenMP プラグマが含まれていることに注意してください。

type codegen/lib/histequalize/histequalize.c
/*
 * File: histequalize.c
 *
 * MATLAB Coder version            : 25.1
 * C/C++ source code generated on  : 13-Jul-2025 16:34:32
 */

/* Include Files */
#include "histequalize.h"
#include "histequalize_data.h"
#include "histequalize_emxutil.h"
#include "histequalize_initialize.h"
#include "histequalize_types.h"
#include "minOrMax.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 tmp_size[3];
  int M;
  int N;
  int loop_ub;
  int plane;
  unsigned int u;
  unsigned int u1;
  int x;
  int y;
  unsigned char tmp_data[24576];
  unsigned char uv[3];
  const unsigned char *originalImage_data;
  unsigned char L_tmp;
  unsigned char r;
  originalImage_data = originalImage->data;
  maximum(originalImage, tmp_data, tmp_size);
  b_maximum(tmp_data, tmp_size, uv);
  L_tmp = c_maximum(uv);
  L = (double)L_tmp + 1.0;
  originalHist_size[0] = 3;
  originalHist_size[1] = L_tmp + 1;
  N = 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, u, u1)

  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 < N; y++) {
      for (x = 0; x < M; x++) {
        r = originalImage_data[(y + originalImage->size[0] * x) +
                               originalImage->size[0] * originalImage->size[1] *
                                   plane];
        u = r + 1U;
        if (r + 1U > 255U) {
          u = 255U;
        }
        u1 = r + 1U;
        if (r + 1U > 255U) {
          u1 = 255U;
        }
        planeHist_data[(int)u - 1] = planeHist_data[(int)u1 - 1] + 1.0;
      }
    }
    loop_ub = originalHist_size[1];
    for (y = 0; y < loop_ub; y++) {
      originalHist_data[plane + 3 * y] = planeHist_data[y];
    }
  }
  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]
 */