Этот пример показывает, как сгенерировать автономную библиотеку 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 : 4.2
* C/C++ source code generated on : 21-Feb-2019 18:00:45
*/
/* Include Files */
#include <string.h>
#include <math.h>
#include "histequalize.h"
#include "histequalize_emxutil.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 int
originalHist_size[2], 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;
unsigned char maxval_data[24576];
int xPageOffset;
unsigned char maxval;
unsigned char b_maxval[3];
int i;
int xOffset;
unsigned int sz_idx_0;
unsigned int sz_idx_1;
int plane;
unsigned char r;
double planeHist_data[256];
int y;
int x;
int i0;
int i1;
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];
npages *= 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;
sz_idx_0 = (unsigned int)originalImage->size[0];
sz_idx_1 = (unsigned int)originalImage->size[1];
vlen = (int)sz_idx_0;
npages = (int)sz_idx_1;
#pragma omp parallel for \
num_threads(omp_get_max_threads()) \
private(r,planeHist_data,y,x,i0,i1)
for (plane = 0; plane < 3; plane++) {
memset(&planeHist_data[0], 0, (unsigned int)((int)*L * (int)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];
i0 = (int)(r + 1U);
i1 = i0;
if ((unsigned int)i0 > 255U) {
i1 = 255;
i0 = 255;
}
planeHist_data[i1 - 1] = planeHist_data[i0 - 1] + 1.0;
}
}
y = (int)*L;
for (i0 = 0; i0 < y; i0++) {
originalHist_data[plane + 3 * i0] = planeHist_data[i0];
}
}
}
/*
* Arguments : double L
* const double originalHist_data[]
* const int originalHist_size[2]
* const emxArray_uint8_T *originalImage
* emxArray_uint8_T *equalizedImage
* Return Type : void
*/
static void equalize(double L, const double originalHist_data[], const int
originalHist_size[2], 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;
double planeHist_data[256];
int loop_ub;
int i2;
int x;
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,planeHist_data,loop_ub,i2,x,j)
for (plane = 0; plane < 3; plane++) {
loop_ub = originalHist_size[1];
for (i2 = 0; i2 < loop_ub; i2++) {
planeHist_data[i2] = originalHist_data[plane + 3 * i2];
}
for (loop_ub = 0; loop_ub < N; loop_ub++) {
for (x = 0; x < M; x++) {
r = originalImage->data[(loop_ub + originalImage->size[0] * x) +
originalImage->size[0] * originalImage->size[1] * plane];
s = 0.0;
i2 = r;
for (j = 0; j <= i2; j++) {
s += planeHist_data[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[(loop_ub + 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];
computeHistogram(originalImage, &L, originalHist_data, originalHist_size);
equalize(L, originalHist_data, originalHist_size, originalImage,
equalizedImage);
}
/*
* File trailer for histequalize.c
*
* [EOF]
*/