title: Applying a (Gaussian) Filter on an Image in Parallel - MATLAB
body: <p>I created the following (MATLAB) function to apply Gaussian Filter blur on an image:</p>
  
  <pre><code>function [ mBlurredImage ] = ApplyGaussianBlur( mInputImage, gaussianKernelStd, stdToRadiusFactor )
  
  gaussianBlurRadius  = ceil(stdToRadiusFactor * gaussianKernelStd); % Imitating Photoshop - See Reference
  
  vGaussianKernel = exp(-([-gaussianBlurRadius:gaussianBlurRadius] .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
  vGaussianKernel = vGaussianKernel / sum(vGaussianKernel);
  
  mInputImagePadded   = padarray(mInputImage, [gaussianBlurRadius, gaussianBlurRadius], 'replicate', 'both');
  
  mBlurredImage = conv2(vGaussianKernel, vGaussianKernel.', mInputImagePadded, 'valid');
  
  
  end
  </code></pre>
  
  <p>I'm trying to find the best approach to create a parallel version of it.<br>
  I want to find a method /strategy which applies to OpenMP as well.</p>
  
  <p>I tried padding the image, then divide it to 4 sections and apply the blur on each.<br>
  Then I gathered all pieces.<br>
  Here's the code:</p>
  
  <pre><code>function [ mBlurredImage ] = ApplyGaussianBlurParallel( mInputImage, gaussianKernelStd, stdToRadiusFactor, numThreads )
  
  numRows = size(mInputImage, 1);
  numCols = size(mInputImage, 2);
  
  % mBlurredImage = zeros(numRows, numCols);
  
  gaussianKernelRadius  = ceil(stdToRadiusFactor * gaussianKernelStd); % Imitating Photoshop - See Reference
  
  vGaussianKernel = exp(-([-gaussianKernelRadius:gaussianKernelRadius] .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
  vGaussianKernel = vGaussianKernel / sum(vGaussianKernel);
  
  numRowsPadded = numRows + (2 * gaussianKernelRadius);
  numColsPadded = numCols + (2 * gaussianKernelRadius);
  
  mInputImagePadded   = padarray(mInputImage, [gaussianKernelRadius, gaussianKernelRadius], 'replicate', 'both');
  
  vColIdxImageBlock = round(linspace(1, numCols, (numThreads + 1)));
  
  vFirstColIdxImageBlock = vColIdxImageBlock(1:numThreads);
  % Going form Image Axis to Padded Image Axis
  vFirstColIdxImagePaddedBlock = vFirstColIdxImageBlock + gaussianKernelRadius;
  % Adding Pixels to the left
  vFirstColIdxImagePaddedBlock = vFirstColIdxImagePaddedBlock - gaussianKernelRadius;
  
  vLastColIdxImageBlock = vColIdxImageBlock(2:(numThreads + 1));
  % Going form Image Axis to Padded Image Axis
  vLastColIdxImagePaddedBlock = vLastColIdxImageBlock + gaussianKernelRadius;
  % Adding Pixels to the right
  vLastColIdxImagePaddedBlock = vLastColIdxImagePaddedBlock + gaussianKernelRadius;
  
  vRowsIdxImage = 1:numRows;
  vRowsImagePadded = 1:numRowsPadded;
  
  cImageBlock             = cell(numThreads, 1);
  cImageBlockProcessed    = cell(numThreads, 1);
  
  for iBlockIdx = 1:numThreads
      firstColIdxImagePaddedBlock    = vFirstColIdxImagePaddedBlock(iBlockIdx);
      lastColIdxImagePaddedBlock     = vLastColIdxImagePaddedBlock(iBlockIdx);
      vColsIdxImagePadded = [firstColIdxImagePaddedBlock:lastColIdxImagePaddedBlock];
  
      cImageBlock{iBlockIdx} =  mInputImagePadded(vRowsImagePadded, vColsIdxImagePadded);
  end
  
  parfor iBlockIdx = 1:numThreads
     cImageBlockProcessed{iBlockIdx} = conv2(vGaussianKernel, vGaussianKernel.', cImageBlock{iBlockIdx}, 'valid');
  end
  
  mBlurredImage = zeros(numRows, numCols);
  
  for iBlockIdx = 1:numThreads
      firstColIdxImageBlock    = vFirstColIdxImageBlock(iBlockIdx);
      lastColIdxImageBlock     = vLastColIdxImageBlock(iBlockIdx);
      vColsIdxImage = [firstColIdxImageBlock:lastColIdxImageBlock];
  
      mBlurredImage(vRowsIdxImage, vColsIdxImage) = cImageBlockProcessed{iBlockIdx};
  end
  
  
  end
  </code></pre>
  
  <p>I also created the following script to analyze the performance:</p>
  
  <pre><code>% `ApplyGaussianBlurParallel` Test Case
  clear();
  
  vInputImageSize = [720, 1280, 1920, 2560];
  numIterations   = 20;
  
  vRunTimeParallelGaussianBlur    = zeros(numIterations, length(vInputImageSize));
  vRunTimeSerialGaussianBlur      = zeros(numIterations, length(vInputImageSize));
  
  gaussianKernelStd   = 10;
  stdToRadiusFactor   = 3.5;
  numThreads          = 4;
  
  for iImageSizeIdx = 1:length(vInputImageSize);
      imageSize = vInputImageSize(iImageSizeIdx);
  
      mInputImage = randn(imageSize, 'single');
  
      maxNumCompThreads(1);
  
      for iIter = 1:numIterations
          hTimeStart = tic();
          mBlurredImage1 = ApplyGaussianBlur(mInputImage, gaussianKernelStd, stdToRadiusFactor);
          vRunTimeSerialGaussianBlur(iIter, iImageSizeIdx) = toc(hTimeStart);
      end
  
      maxNumCompThreads(numThreads);
  
      for iIter = 1:numIterations
          hTimeStart = tic();
          mBlurredImage1 = ApplyGaussianBlurParallel(mInputImage, gaussianKernelStd, stdToRadiusFactor, numThreads);
          vRunTimeParallelGaussianBlur(iIter, iImageSizeIdx) = toc(hTimeStart);
      end
  
  end
  
  vRunTimeParallelGaussianBlurMean    = mean(vRunTimeParallelGaussianBlur);
  vRunTimeParallelGaussianBlurStd     = std(vRunTimeParallelGaussianBlur);
  vRunTimeParallelGaussianBlurMedian  = median(vRunTimeParallelGaussianBlur);
  
  vRunTimeSerialGaussianBlurMean    = mean(vRunTimeSerialGaussianBlur);
  vRunTimeSerialGaussianBlurStd     = std(vRunTimeSerialGaussianBlur);
  vRunTimeSerialGaussianBlurMedian  = median(vRunTimeSerialGaussianBlur);
  
  figure();
  plot(vInputImageSize, [vRunTimeParallelGaussianBlurMean(:), vRunTimeSerialGaussianBlurMean(:)], ...
      'LineStyle', 'none', 'Marker', 'o');
  title('Mean Runtime');
  legend({['Parallel'], ['Serial']});
  
  figure();
  plot(vInputImageSize, [vRunTimeParallelGaussianBlurMedian(:), vRunTimeSerialGaussianBlurMedian(:)], ...
      'LineStyle', 'none', 'Marker', 'o');
  title('Median Runtime');
  legend({['Parallel'], ['Serial']});
  </code></pre>
  
  <p>Yet what I get is:</p>
  
  <p><img src="https://i.stack.imgur.com/6AS0J.png" alt="Mean of Run Time">
  <img src="https://i.stack.imgur.com/CNBIH.png" alt="Median of Run Time"></p>
  
  <p>Namely, I can't make it efficient enough.<br>
  Can anyone think on a better more efficient approach or doing it better?  </p>
  
  <p>Thank You.</p>