MATLAB Examples

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.

Contents

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.

Create a New Folder and Copy Relevant Files

The following code will create a folder in your current working folder (pwd). The new folder will only contain the files that are relevant for this example. If you do not want to affect the current folder (or if you cannot generate files in this folder), you should change your working folder.

Run Command: Create a New Folder and Copy Relevant Files

coderdemo_setup('coderdemo_contrast_enhancer');

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
Warning: The selected compiler does not support the OpenMP library; this loop
will not be parallelized.

Warning in ==> histequalize Line: 22 Column: 5
Code generation successful (with warnings): To view the report, open('codegen/mex/histequalize/html/index.html').

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
Warning: The selected compiler does not support the OpenMP library; this loop
will not be parallelized.

Warning in ==> histequalize Line: 22 Column: 5
Code generation successful (with warnings): To view the report, open('codegen/lib/histequalize/html/index.html').

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
/*
 * File: histequalize.c
 *
 * MATLAB Coder version            : 3.4
 * C/C++ source code generated on  : 25-Jul-2017 16:03:46
 */

/* Include Files */
#include "rt_nonfinite.h"
#include "histequalize.h"
#include "histequalize_emxutil.h"

/* Function Declarations */
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    : 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 y;
  double normalizer;
  unsigned int sz[3];
  int loop_ub;
  int plane;
  double planeHist_data[256];
  int x;
  double s;
  int j;
  unsigned char u0;
  y = 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, y);
  for (y = 0; y < 3; y++) {
    sz[y] = (unsigned int)originalImage->size[y];
  }

  normalizer = (L - 1.0) / ((double)sz[0] * (double)sz[1]);
  loop_ub = originalHist_size[1];
  for (plane = 0; plane < 3; plane++) {
    for (y = 0; y < loop_ub; y++) {
      planeHist_data[y] = originalHist_data[plane + originalHist_size[0] * y];
    }

    for (y = 0; y < (int)sz[0]; y++) {
      for (x = 0; x < (int)sz[1]; x++) {
        s = 0.0;
        for (j = 0; j <= originalImage->data[(y + originalImage->size[0] * x) +
             originalImage->size[0] * originalImage->size[1] * plane]; j++) {
          s += planeHist_data[j];
        }

        s *= normalizer;
        s = rt_roundd_snf(s);
        if (s < 256.0) {
          if (s >= 0.0) {
            u0 = (unsigned char)s;
          } else {
            u0 = 0;
          }
        } else if (s >= 256.0) {
          u0 = MAX_uint8_T;
        } else {
          u0 = 0;
        }

        equalizedImage->data[(y + equalizedImage->size[0] * x) +
          equalizedImage->size[0] * equalizedImage->size[1] * plane] = u0;
      }
    }
  }
}

/*
 * 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)
{
  int maxval_size_idx_1;
  int npages;
  int ix;
  int i;
  int ixstop;
  unsigned char maxval_data[24576];
  unsigned char mtmp;
  unsigned char maxval[3];
  int originalHist_size[2];
  unsigned int sz[3];
  double originalHist_data[768];
  double planeHist_data[256];
  int x;
  int i0;
  maxval_size_idx_1 = originalImage->size[1];
  npages = 1;
  for (ix = 0; ix < 2; ix++) {
    npages *= originalImage->size[ix + 1];
  }

  for (ix = 0; ix + 1 <= npages; ix++) {
    ixstop = originalImage->size[0];
    maxval_data[ix] = originalImage->data[ixstop * ix];
    for (i = 1; i + 1 <= originalImage->size[0]; i++) {
      ixstop = originalImage->size[0];
      if (originalImage->data[i + ixstop * ix] > maxval_data[ix]) {
        ixstop = originalImage->size[0];
        maxval_data[ix] = originalImage->data[i + ixstop * ix];
      }
    }
  }

  for (i = 0; i < 3; i++) {
    ix = i * maxval_size_idx_1;
    ixstop = ix + maxval_size_idx_1;
    mtmp = maxval_data[ix];
    if ((maxval_size_idx_1 > 1) && (ix + 1 < ixstop)) {
      for (ix = i * maxval_size_idx_1 + 1; ix + 1 <= ixstop; ix++) {
        if (maxval_data[ix] > mtmp) {
          mtmp = maxval_data[ix];
        }
      }
    }

    maxval[i] = mtmp;
  }

  mtmp = maxval[0];
  for (ix = 0; ix < 2; ix++) {
    if (maxval[ix + 1] > mtmp) {
      mtmp = maxval[ix + 1];
    }
  }

  originalHist_size[0] = 3;
  originalHist_size[1] = mtmp + 1;
  for (ix = 0; ix < 3; ix++) {
    sz[ix] = (unsigned int)originalImage->size[ix];
  }

  ixstop = mtmp + 1;
  npages = mtmp + 1;
  for (i = 0; i < 3; i++) {
    memset(&planeHist_data[0], 0, (unsigned int)(ixstop * (int)sizeof(double)));
    for (maxval_size_idx_1 = 0; maxval_size_idx_1 < (int)sz[0];
         maxval_size_idx_1++) {
      for (x = 0; x < (int)sz[1]; x++) {
        ix = (int)(originalImage->data[(maxval_size_idx_1 + originalImage->size
                    [0] * x) + originalImage->size[0] * originalImage->size[1] *
                   i] + 1U);
        if ((unsigned int)ix > 255U) {
          ix = 255;
        }

        i0 = (int)(originalImage->data[(maxval_size_idx_1 + originalImage->size
                    [0] * x) + originalImage->size[0] * originalImage->size[1] *
                   i] + 1U);
        if ((unsigned int)i0 > 255U) {
          i0 = 255;
        }

        planeHist_data[ix - 1] = planeHist_data[i0 - 1] + 1.0;
      }
    }

    for (ix = 0; ix < npages; ix++) {
      originalHist_data[i + 3 * ix] = planeHist_data[ix];
    }
  }

  equalize((double)mtmp + 1.0, originalHist_data, originalHist_size,
           originalImage, equalizedImage);
}

/*
 * File trailer for histequalize.c
 *
 * [EOF]
 */

Cleanup

Remove files and return to original folder

Run Command: Cleanup

cleanup