Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

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
/*
 * Prerelease License - for engineering feedback and testing purposes
 * only. Not for sale.
 * File: histequalize.c
 *
 * MATLAB Coder version            : 23.2
 * C/C++ source code generated on  : 27-Jul-2023 15:00:26
 */

/* 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_size_idx_1;
  int p;
  int plane;
  int vlen;
  int x;
  int xOffset;
  int xPageOffset;
  int y;
  short planeHist_size[2];
  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_size_idx_1 = originalImage->size[1];
  if (originalImage->size[1] == 0) {
    maxval_size_idx_1 = originalImage->size[1];
    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_size_idx_1;
    maxval[p] = maxval_data[xPageOffset];
    for (i = 2; i <= maxval_size_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, planeHist_size, loop_ub, y, x, b_i)

  for (plane = 0; plane < 3; plane++) {
    loop_ub = (int)L;
    planeHist_size[1] = (short)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 = planeHist_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]
 */