Использование PARFOR для ускорения алгоритма улучшения контрастности изображений

В этом примере показано, как сгенерировать автономную библиотеку C из кода MATLAB ®, которая применяет простую функцию эквализации гистограммы к изображениям для улучшения контрастности изображения. В примере используются parfor для обработки каждой из стандартных трех плоскостей изображений RGB на отдельных потоках. Пример также показывает, как сгенерировать и запустить MEX-функцию в MATLAB до генерации кода С, чтобы убедиться, что код MATLAB подходит для генерации кода.

MATLAB Coder™ использует стандарт параллельного программирования общей памяти OpenMP для реализации своей поддержки parfor. Смотрите Спецификацию API OpenMP для параллельного программирования. В то время как MATLAB поддерживает parfor создав несколько рабочие сеансы, MATLAB Coder использует OpenMP для создания нескольких потоков, выполняемых на одной машине.

Необходимые условия

В порядок поддержки параллелизации компилятор должен поддерживать стандарт параллельного программирования общей памяти OpenMP. Если у вашего компилятора нет этой поддержки, то вы все равно можете запустить этот пример, но сгенерированный код будет запускаться последовательно.

О histequalize Функция

The histequalize.m функция принимает изображение (представленное как NxMx3 матрица) и возвращает изображение с усиленной контрастностью.

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-функцию

Сгенерируйте MEX-функцию с помощью codegen команда.

codegen histequalize
Code generation successful.

Перед генерацией Кода С необходимо сначала протестировать MEX-функцию в MATLAB, чтобы убедиться, что она функционально эквивалентна оригинальному коду MATLAB и что никаких ошибок времени выполнения не происходит. По умолчанию codegen генерирует MEX-функцию с именем histequalize_mex в текущей папке. Это позволяет вам протестировать код MATLAB и MEX-функцию и сравнить результаты.

Чтение в оригинальном изображении

Используйте стандартную imread команда для чтения низкоконтрастного изображения.

lcIm = imread('LowContrast.jpg');
image(lcIm);

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

Запустите MEX-функцию (алгоритм Эквализации гистограммы)

Передайте низкоконтрастное изображение.

hcIm = histequalize_mex(lcIm);

Отображение результата

image(hcIm);

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

Сгенерируйте автономный код С

codegen -config:lib histequalize
Code generation successful.

Использование codegen с -config:lib опция создает автономную библиотеку C. По умолчанию код, сгенерированный для библиотеки, находится в папке codegen/lib/histequalize/.

Смотрите сгенерированную функцию

Заметьте, что сгенерированный код содержит прагмы OpenMP, которые управляют параллелизацией кода с помощью нескольких потоков.

type codegen/lib/histequalize/histequalize.c
/*
 * File: histequalize.c
 *
 * MATLAB Coder version            : 5.2
 * C/C++ source code generated on  : 21-Apr-2021 01:26:27
 */

/* Include Files */
#include "histequalize.h"
#include "histequalize_data.h"
#include "histequalize_emxutil.h"
#include "histequalize_initialize.h"
#include "histequalize_types.h"
#include <math.h>
#include <string.h>

/* Function Declarations */
static void computeHistogram(const emxArray_uint8_T *originalImage, double *L,
                             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 *L
 *                double originalHist_data[]
 *                int originalHist_size[2]
 * Return Type  : void
 */
static void computeHistogram(const emxArray_uint8_T *originalImage, double *L,
                             double originalHist_data[],
                             int originalHist_size[2])
{
  double planeHist_data[256];
  int b_i;
  int i;
  int loop_ub;
  int maxval_size_idx_1;
  int npages;
  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 b_maxval[3];
  unsigned char maxval;
  unsigned char r;
  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 (0 <= vlen - 1) {
      memset(&maxval_data[0], 0, vlen * sizeof(unsigned char));
    }
  } else {
    vlen = originalImage->size[0];
    npages = originalImage->size[1] * 3;
    for (p = 0; p < npages; 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;
    b_maxval[p] = maxval_data[xPageOffset];
    for (i = 2; i <= maxval_size_idx_1; i++) {
      xOffset = (xPageOffset + i) - 1;
      if (b_maxval[p] < maxval_data[xOffset]) {
        b_maxval[p] = maxval_data[xOffset];
      }
    }
  }
  maxval = b_maxval[0];
  if (b_maxval[0] < b_maxval[1]) {
    maxval = b_maxval[1];
  }
  if (maxval < b_maxval[2]) {
    maxval = b_maxval[2];
  }
  *L = (double)maxval + 1.0;
  originalHist_size[0] = 3;
  originalHist_size[1] = maxval + 1;
  vlen = originalImage->size[0];
  npages = 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;
    if (0 <= loop_ub - 1) {
      memset(&planeHist_data[0], 0, loop_ub * sizeof(double));
    }
    for (y = 0; y < vlen; y++) {
      for (x = 0; x < npages; 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[(unsigned char)b_i - 1] =
            planeHist_data[(unsigned char)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];
    }
  }
}

/*
 * 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;
  unsigned char r;
  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);
  N = originalImage->size[0];
  M = originalImage->size[1];
  normalizer = (L - 1.0) / ((double)(unsigned int)originalImage->size[0] *
                            (double)(unsigned int)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();
  }
  computeHistogram(originalImage, &L, originalHist_data, originalHist_size);
  equalize(L, originalHist_data, originalImage, equalizedImage);
}

/*
 * File trailer for histequalize.c
 *
 * [EOF]
 */