В этом примере показано, как сгенерировать автономную библиотеку C из кода MATLAB, который применяет простую функцию эквализации гистограммы к изображениям, чтобы улучшить контрастность изображений. Пример использует parfor
обработать каждую из стандартных трех плоскостей изображения RGB на отдельных потоках. Пример также показывает, как сгенерировать и запустить MEX-функцию в MATLAB до генерации кода С, чтобы проверить, что код MATLAB подходит для генерации кода.
MATLAB Coder использует OpenMP портативный стандарт параллельного программирования общей памяти, чтобы реализовать его поддержку parfor
. См. спецификацию API OpenMP для Параллельного программирования. Принимая во внимание, что MATLAB поддерживает parfor
путем создания нескольких сеансов рабочего MATLAB Coder использует OpenMP, чтобы создать несколько потоков, работающих на той же машине.
Для того, чтобы поддержать распараллеливание, компилятор должен поддержать стандарт параллельного программирования общей памяти OpenMP. Если ваш компилятор не имеет этой поддержки, то можно все еще запустить этот пример, но сгенерированный код запустится последовательно.
histequalize
Функцияhistequalize.m
функционируйте берет изображение (представленный как матрица NxMx3) и возвращает изображение с расширенным контрастом.
type histequalize
function 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.1 * C/C++ source code generated on : 17-Aug-2020 20:30:47 */ /* 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] * 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,planeHist_size,loop_ub,y,x,b_i) for (plane = 0; plane < 3; plane++) { planeHist_size[1] = (short)*L; 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); 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] = originalImage->size[2]; 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] */