Использование 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 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.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]
*/