В этом примере показано, как создать автономную библиотеку C из кода MATLAB ®, которая применяет простую функцию выравнивания гистограммы к изображениям для улучшения контрастности изображения. В примере используетсяparfor для обработки каждой из трех стандартных плоскостей изображения RGB на отдельных потоках. В примере также показано, как генерировать и запускать функцию MEX в MATLAB перед генерацией кода C для проверки того, что код 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 histequalizeCode generation successful.
Перед созданием кода C необходимо сначала протестировать функцию 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
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]
*/