В этом примере показано, как сгенерировать автономную библиотеку 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-функцию с помощью codegen команда.
codegen histequalizeПрежде, чем сгенерировать код С, необходимо сначала протестировать MEX-функцию в MATLAB, чтобы гарантировать, что это функционально эквивалентно оригинальному коду MATLAB и что никакие ошибки времени выполнения не происходят. По умолчанию, codegen генерирует MEX-функцию под названием histequalize_mex в текущей папке. Это позволяет вам тестировать код MATLAB и MEX-функцию и сравнивать результаты.
Используйте стандартный imread команда, чтобы считать низкоконтрастное изображение.
lcIm = imread('LowContrast.jpg');
image(lcIm);
Передайте низкоконтрастное изображение.
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]
*/