使用 PARFOR 加速图像对比度增强算法

打开实时脚本

此示例说明如何从 MATLAB® 代码生成独立的 C 库，该库对图像应用简单的直方图均衡化函数来提高图像对比度。该示例使用 parfor 在各个单独线程上分别处理三个标准 RGB 图像平面之一。示例还说明如何在生成 C 代码之前在 MATLAB 中生成和运行 MEX 函数，以验证该 MATLAB 代码适合进行代码生成。

MATLAB Coder™ 使用 OpenMP 可移植共享内存并行编程标准来实现它对 parfor 的支持。请参阅 The OpenMP API Specification for Parallel Programming。MATLAB 通过创建多个工作进程会话来支持 parfor，而 MATLAB Coder 使用 OpenMP 来创建在同一计算机上运行的多个线程。

前提条件

为了支持并行化，编译器必须支持 OpenMP 共享内存并行编程标准。如果您的编译器不支持此功能，您仍可以运行此示例，但生成的代码将以串行方式运行。

关于 `histequalize` 函数

histequalize.m 函数接受一个图像（表示为 N×M×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.

将 codegen 与 -config:lib 选项结合使用生成独立 C 库。默认情况下，为库生成的代码位于文件夹 codegen/lib/histequalize/ 中。

检查生成的函数

请注意，生成的代码包含 OpenMP pragma，这些 pragma 使用多个线程控制代码的并行化。

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.2
 * C/C++ source code generated on  : 20-Jul-2024 12:22:12
 */

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