Используя PARFOR, чтобы ускорить алгоритм улучшения контрастности изображений
В этом примере показано, как сгенерировать автономную библиотеку C из кода MATLAB®, который применяет простую функцию эквализации гистограммы к изображениям, чтобы улучшить контрастность изображений. Пример использует parfor обработать каждую из стандартных трех плоскостей изображения RGB на отдельных потоках. Пример также показывает, как сгенерировать и запустить MEX-функцию в MATLAB до генерации кода С, чтобы проверить, что код MATLAB подходит для генерации кода.
MATLAB Coder™ использует OpenMP портативный стандарт параллельного программирования общей памяти, чтобы реализовать его поддержку parfor. См. спецификацию API OpenMP для Параллельного программирования. Принимая во внимание, что MATLAB поддерживает parfor путем создания нескольких сеансов рабочего MATLAB Coder использует OpenMP, чтобы создать несколько потоков, работающих на той же машине.
Необходимые условия
Для того, чтобы поддержать распараллеливание, компилятор должен поддержать стандарт параллельного программирования общей памяти OpenMP. Если ваш компилятор не имеет этой поддержки, то можно все еще запустить этот пример, но сгенерированный код запустится последовательно.
О histequalize Функция
histequalize.m функционируйте берет изображение (представленный как матрица NxMx3) и возвращает изображение с расширенным контрастом.
type histequalizefunction 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 histequalizeCode generation successful.
Прежде, чем сгенерировать код С, необходимо сначала протестировать 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
Code generation successful.
Используя codegen с -config:lib опция производит автономную библиотеку C. По умолчанию код, сгенерированный для библиотеки, находится в папке codegen/lib/histequalize/.
Смотрите сгенерированную функцию
Заметьте, что сгенерированный код содержит прагмы OpenMP, которые управляют распараллеливанием кода с помощью нескольких потоков.
type codegen/lib/histequalize/histequalize.c/*
* File: histequalize.c
*
* MATLAB Coder version : 5.3
* C/C++ source code generated on : 25-Aug-2021 05:34:03
*/
/* 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];
const unsigned char *originalImage_data;
unsigned char 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 (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;
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)(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]
*/