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

Скрипт Open Live Script

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

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

Предпосылки

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

О `histequalize` Функция

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

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

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

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

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

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

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

hcIm = histequalize_mex(lcIm);

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

image(hcIm);

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

codegen -config:lib histequalize

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

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

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

type codegen/lib/histequalize/histequalize.c

/*
 * File: histequalize.c
 *
 * MATLAB Coder version            : 5.0
 * C/C++ source code generated on  : 29-Jan-2020 14:23:48
 */

/* Include Files */
#include "histequalize.h"
#include "histequalize_data.h"
#include "histequalize_emxutil.h"
#include "histequalize_initialize.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])
{
  int maxval_size_idx_1;
  int vlen;
  int npages;
  int p;
  int xPageOffset;
  unsigned char maxval_data[24576];
  unsigned char maxval;
  unsigned char b_maxval[3];
  int i;
  int xOffset;
  int plane;
  unsigned char r;
  double planeHist_data[256];
  int loop_ub;
  int y;
  int x;
  int b_i;
  maxval_size_idx_1 = originalImage->size[1];
  if (originalImage->size[1] == 0) {
    for (vlen = 0; vlen < 3; vlen++) {
      for (npages = 0; npages < maxval_size_idx_1; npages++) {
        maxval_data[npages + maxval_size_idx_1 * vlen] = 0U;
      }
    }
  } else {
    vlen = originalImage->size[0];
    npages = originalImage->size[1] * originalImage->size[2];
    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,loop_ub,y,x,b_i)

  for (plane = 0; plane < 3; plane++) {
    loop_ub = (int)*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);
        loop_ub = b_i;
        if ((unsigned int)b_i > 255U) {
          loop_ub = 255;
          b_i = 255;
        }

        planeHist_data[(unsigned char)loop_ub - 1] = planeHist_data[(unsigned
          char)b_i - 1] + 1.0;
      }
    }

    loop_ub = (short)(int)*L;
    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)
{
  int N;
  unsigned int sz_idx_0;
  unsigned int sz_idx_1;
  int M;
  double normalizer;
  int plane;
  double s;
  unsigned char r;
  int y;
  int x;
  int i;
  int j;
  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] = originalImage->size[2];
  emxEnsureCapacity_uint8_T(equalizedImage, N);
  sz_idx_0 = (unsigned int)originalImage->size[0];
  sz_idx_1 = (unsigned int)originalImage->size[1];
  N = (int)sz_idx_0;
  M = (int)sz_idx_1;
  normalizer = (L - 1.0) / ((double)sz_idx_0 * (double)sz_idx_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 L;
  double originalHist_data[768];
  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]
 */

Документация

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

Предпосылки

О `histequalize` Функция

Сгенерируйте MEX-функцию

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

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

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

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

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

Документация MATLAB Coder

Поддержка

Документация

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

Предпосылки

О histequalize Функция

Сгенерируйте MEX-функцию

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

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

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

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

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

Документация MATLAB Coder

Поддержка

О `histequalize` Функция