Using PARFOR to Speed Up an Image Contrast Enhancement Algorithm
This example shows how to generate a standalone C library from MATLAB® code that applies a simple histogram equalization function to images to improve image contrast. The example uses parfor
to process each of the standard three RGB image planes on separate threads. The example also shows how to generate and run a MEX function in MATLAB prior to generating C code to verify that the MATLAB code is suitable for code generation.
MATLAB Coder™ uses the OpenMP portable shared memory parallel programming standard to implement its support for parfor
. See The OpenMP API Specification for Parallel Programming. Whereas MATLAB supports parfor
by creating multiple worker sessions, MATLAB Coder uses OpenMP to create multiple threads running on the same machine.
Prerequisites
In order to support parallelization, the compiler must support the OpenMP shared memory parallel programming standard. If your compiler does not have this support, then you can still run this example, but the generated code will run serially.
About the histequalize
Function
The histequalize.m
function takes an image (represented as an NxMx3 matrix) and returns an image with enhanced contrast.
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
Generate the MEX Function
Generate a MEX function using the codegen
command.
codegen histequalize
Code generation successful.
Before generating C code, you should first test the MEX function in MATLAB to ensure that it is functionally equivalent to the original MATLAB code and that no run-time errors occur. By default, codegen
generates a MEX function named histequalize_mex
in the current folder. This allows you to test the MATLAB code and MEX function and compare the results.
Read in the Original Image
Use the standard imread
command to read the low-contrast image.
lcIm = imread('LowContrast.jpg');
image(lcIm);
Run the MEX Function (The Histogram Equalization Algorithm)
Pass the low-contrast image.
hcIm = histequalize_mex(lcIm);
Display the Result
image(hcIm);
Generate Standalone C Code
codegen -config:lib histequalize
Code generation successful.
Using codegen
with the -config:lib
option produces a standalone C library. By default, the code generated for the library is in the folder codegen/lib/histequalize/
.
Inspect the Generated Function
Notice that the generated code contains OpenMP pragmas that control parallelization of the code using multiple threads.
type codegen/lib/histequalize/histequalize.c
/* * Prerelease License - for engineering feedback and testing purposes * only. Not for sale. * File: histequalize.c * * MATLAB Coder version : 24.2 * C/C++ source code generated on : 20-Jul-2024 12:22:12 */ /* Include Files */ #include "histequalize.h" #include "histequalize_data.h" #include "histequalize_emxutil.h" #include "histequalize_initialize.h" #include "histequalize_types.h" #include "minOrMax.h" #include "omp.h" #include <math.h> #include <string.h> /* Function Declarations */ static double computeHistogram(const emxArray_uint8_T *originalImage, 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 originalHist_data[] * int originalHist_size[2] * Return Type : double */ static double computeHistogram(const emxArray_uint8_T *originalImage, double originalHist_data[], int originalHist_size[2]) { double planeHist_data[256]; double L; int tmp_size[3]; int M; int N; int i; int loop_ub; int plane; int x; int y; unsigned char tmp_data[24576]; unsigned char uv[3]; const unsigned char *originalImage_data; unsigned char L_tmp; unsigned char r; originalImage_data = originalImage->data; maximum(originalImage, tmp_data, tmp_size); b_maximum(tmp_data, tmp_size, uv); L_tmp = c_maximum(uv); L = (double)L_tmp + 1.0; originalHist_size[0] = 3; originalHist_size[1] = L_tmp + 1; N = originalImage->size[0]; M = originalImage->size[1]; #pragma omp parallel for num_threads(omp_get_max_threads()) private( \ r, planeHist_data, loop_ub, y, x, i) for (plane = 0; plane < 3; plane++) { loop_ub = (int)L; memset(&planeHist_data[0], 0, (unsigned int)loop_ub * sizeof(double)); 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]; i = (int)(r + 1U); if (r + 1U > 255U) { i = 255; } loop_ub = (int)(r + 1U); if (r + 1U > 255U) { loop_ub = 255; } planeHist_data[i - 1] = planeHist_data[loop_ub - 1] + 1.0; } } loop_ub = originalHist_size[1]; for (i = 0; i < loop_ub; i++) { originalHist_data[plane + 3 * i] = planeHist_data[i]; } } return L; } /* * 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; const unsigned char *originalImage_data; unsigned char r; unsigned char *equalizedImage_data; originalImage_data = originalImage->data; 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); equalizedImage_data = equalizedImage->data; N = originalImage->size[0]; M = originalImage->size[1]; normalizer = (L - 1.0) / ((double)originalImage->size[0] * (double)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(); } L = computeHistogram(originalImage, originalHist_data, originalHist_size); equalize(L, originalHist_data, originalImage, equalizedImage); } /* * File trailer for histequalize.c * * [EOF] */