From 0cc860707a52b9d76275f92bda160566e7de4474 Mon Sep 17 00:00:00 2001
From: Jun Saito <dukecyto@gmail.com>
Date: Tue, 19 Sep 2017 09:59:05 -0700
Subject: [PATCH] OpenMP multithreading

---
 setup.py              |  6 +++-
 src/ImageProcessing.h | 78 +++++++++++++++++++------------------------
 src/OpticalFlow.cpp   |  1 +
 3 files changed, 41 insertions(+), 44 deletions(-)

diff --git a/setup.py b/setup.py
index c3ce742..479a59a 100644
--- a/setup.py
+++ b/setup.py
@@ -13,7 +13,11 @@
 
 sourcefiles = ['pyflow.pyx', ]
 sourcefiles.extend(glob("src/*.cpp"))
-extensions = [Extension("pyflow", sourcefiles, include_dirs=[numpy.get_include()])]
+extensions = [Extension("pyflow",
+                        sourcefiles,
+                        include_dirs=[numpy.get_include()],
+                        extra_compile_args=['-fopenmp'],
+                        extra_link_args=['-lgomp'])]
 setup(
     name="pyflow",
     version="1.0",
diff --git a/src/ImageProcessing.h b/src/ImageProcessing.h
index ec9fd3b..7e2fdf5 100644
--- a/src/ImageProcessing.h
+++ b/src/ImageProcessing.h
@@ -249,19 +249,17 @@ template <class T1,class T2>
 void ImageProcessing::hfiltering(const T1* pSrcImage,T2* pDstImage,int width,int height,int nChannels,const double* pfilter1D,int fsize)
 {
 	memset(pDstImage,0,sizeof(T2)*width*height*nChannels);
-	T2* pBuffer;
-	double w;
-	int i,j,l,k,offset,jj;
-	for(i=0;i<height;i++)
-		for(j=0;j<width;j++)
+#pragma omp parallel for
+	for(int i=0;i<height;i++)
+		for(int j=0;j<width;j++)
 		{
-			offset=i*width*nChannels;
-			pBuffer=pDstImage+offset+j*nChannels;
-			for(l=-fsize;l<=fsize;l++)
+			int offset=i*width*nChannels;
+			T2* pBuffer=pDstImage+offset+j*nChannels;
+			for(int l=-fsize;l<=fsize;l++)
 			{
-				w=pfilter1D[l+fsize];
-				jj=EnforceRange(j+l,width);
-				for(k=0;k<nChannels;k++)
+				double w=pfilter1D[l+fsize];
+				int jj=EnforceRange(j+l,width);
+				for(int k=0;k<nChannels;k++)
 					pBuffer[k]+=pSrcImage[offset+jj*nChannels+k]*w;
 			}
 		}
@@ -274,20 +272,18 @@ template <class T1,class T2>
 void ImageProcessing::hfiltering_transpose(const T1* pSrcImage,T2* pDstImage,int width,int height,int nChannels,const double* pfilter1D,int fsize)
 {
 	memset(pDstImage,0,sizeof(T2)*width*height*nChannels);
-	const T1* pBuffer;
-	double w;
-	int i,j,l,k,offset,jj;
-	for(i=0;i<height;i++)
-		for(j=0;j<width;j++)
+#pragma omp parallel for
+	for(int i=0;i<height;i++)
+		for(int j=0;j<width;j++)
 		{
 			int offset0=i*width*nChannels;
-			pBuffer=pSrcImage+(i*width+j)*nChannels;
-			for(l=-fsize;l<=fsize;l++)
+			const T1* pBuffer=pSrcImage+(i*width+j)*nChannels;
+			for(int l=-fsize;l<=fsize;l++)
 			{
-				w=pfilter1D[l+fsize];
-				jj=EnforceRange(j+l,width);
-				offset = offset0 + jj*nChannels;
-				for(k=0;k<nChannels;k++)
+				double w=pfilter1D[l+fsize];
+				int jj=EnforceRange(j+l,width);
+				int offset = offset0 + jj*nChannels;
+				for(int k=0;k<nChannels;k++)
 					pDstImage[offset+k] += pBuffer[k]*w;
 			}
 		}
@@ -339,18 +335,16 @@ template <class T1,class T2>
 void ImageProcessing::vfiltering(const T1* pSrcImage,T2* pDstImage,int width,int height,int nChannels,const double* pfilter1D,int fsize)
 {
 	memset(pDstImage,0,sizeof(T2)*width*height*nChannels);
-	T2* pBuffer;
-	double w;
-	int i,j,l,k,offset,ii;
-	for(i=0;i<height;i++)
-		for(j=0;j<width;j++)
+#pragma omp parallel for
+	for(int i=0;i<height;i++)
+		for(int j=0;j<width;j++)
 		{
-			pBuffer=pDstImage+(i*width+j)*nChannels;
-			for(l=-fsize;l<=fsize;l++)
+			T2* pBuffer=pDstImage+(i*width+j)*nChannels;
+			for(int l=-fsize;l<=fsize;l++)
 			{
-				w=pfilter1D[l+fsize];
-				ii=EnforceRange(i+l,height);
-				for(k=0;k<nChannels;k++)
+				double w=pfilter1D[l+fsize];
+				int ii=EnforceRange(i+l,height);
+				for(int k=0;k<nChannels;k++)
 					pBuffer[k]+=pSrcImage[(ii*width+j)*nChannels+k]*w;
 			}
 		}
@@ -363,19 +357,17 @@ template <class T1,class T2>
 void ImageProcessing::vfiltering_transpose(const T1* pSrcImage,T2* pDstImage,int width,int height,int nChannels,const double* pfilter1D,int fsize)
 {
 	memset(pDstImage,0,sizeof(T2)*width*height*nChannels);
-	const T1* pBuffer;
-	double w;
-	int i,j,l,k,offset,ii;
-	for(i=0;i<height;i++)
-		for(j=0;j<width;j++)
+#pragma omp parallel for
+	for(int i=0;i<height;i++)
+		for(int j=0;j<width;j++)
 		{
-			pBuffer=pSrcImage+(i*width+j)*nChannels;
-			for(l=-fsize;l<=fsize;l++)
+			const T1* pBuffer=pSrcImage+(i*width+j)*nChannels;
+			for(int l=-fsize;l<=fsize;l++)
 			{
-				w=pfilter1D[l+fsize];
-				ii=EnforceRange(i+l,height);
-				offset = (ii*width+j)*nChannels;
-				for(k=0;k<nChannels;k++)
+				double w=pfilter1D[l+fsize];
+				int ii=EnforceRange(i+l,height);
+				int offset = (ii*width+j)*nChannels;
+				for(int k=0;k<nChannels;k++)
 					//pBuffer[k]+=pSrcImage[(ii*width+j)*nChannels+k]*w;
 					pDstImage[offset+k] += pBuffer[k]*w;
 			}
diff --git a/src/OpticalFlow.cpp b/src/OpticalFlow.cpp
index 058d0b3..074b82f 100644
--- a/src/OpticalFlow.cpp
+++ b/src/OpticalFlow.cpp
@@ -381,6 +381,7 @@ void OpticalFlow::SmoothFlowSOR(const DImage &Im1, const DImage &Im2, DImage &wa
 			dv.reset();
 
 			for(int k = 0; k<nSORIterations; k++)
+#pragma omp parallel for
 				for(int i = 0; i<imHeight; i++)
 					for(int j = 0; j<imWidth; j++)
 					{