利用OpenMP加速拉伸图像操作

前面的博客《OpenCV拉伸图像》介绍了如何利用OpenCV的现成函数实现图像的透视变换。本文受到了《http://blog.csdn.net/xiaowei_cqu/article/details/26471527》的启发，利用OpenMP实现对透视变换的加速。

1）启用OpenMP。openmp的使用很简单，只要是vs2010或以上的环境都可以使用。在需要使用的文件开头包含 <omp.h>，然后在“项目属性>c/c++>语言" 里选择openMP支持：

2）根据http://blog.csdn.net/xiaowei_cqu/article/details/26471527的代码，我做了一些修改，利用宏定义 "#prgma omp parallel for"并行加速。代码的具体含义请参见http://blog.csdn.net/xiaowei_cqu/article/details/26471527。

#include <vector>
#include <omp.h>

class PerspectiveTransform
{
public:
	PerspectiveTransform::PerspectiveTransform(float inA11, float inA21, 
                                           float inA31, float inA12, 
                                           float inA22, float inA32, 
                                           float inA13, float inA23, 
                                           float inA33);

	float a11, a21, a31, a12, a22, a32, a13, a23, a33;

	static PerspectiveTransform quadrilateralToQuadrilateral(float x0, float y0, float x1, float y1,
														float x2, float y2, float x3, float y3, float x0p, float y0p, 
														float x1p, float y1p, float x2p, float y2p,
														float x3p, float y3p);
	static PerspectiveTransform squareToQuadrilateral(float x0, float y0, float x1, float y1, float x2,
												float y2, float x3, float y3);
	static PerspectiveTransform quadrilateralToSquare(float x0, float y0, float x1, float y1, float x2,
												float y2, float x3, float y3);
	PerspectiveTransform buildAdjoint();
	PerspectiveTransform times(PerspectiveTransform other);
	void transformPoints(std::vector<float> &points);
	void pntProjection(unsigned char *pSrc, int iSrcWidth, int iSrcHeight,
										 unsigned char * pDst, int iDstWidth, int iDstHeight);
};

#include "Perspective.h"


PerspectiveTransform::PerspectiveTransform(float inA11, float inA21, 
                                           float inA31, float inA12, 
                                           float inA22, float inA32, 
                                           float inA13, float inA23, 
                                           float inA33) : 
  a11(inA11), a12(inA12), a13(inA13), a21(inA21), a22(inA22), a23(inA23),
  a31(inA31), a32(inA32), a33(inA33) {}

PerspectiveTransform PerspectiveTransform::quadrilateralToQuadrilateral(float x0, float y0, float x1, float y1,
    float x2, float y2, float x3, float y3, float x0p, float y0p, float x1p, float y1p, float x2p, float y2p,
    float x3p, float y3p) {
  PerspectiveTransform qToS = PerspectiveTransform::quadrilateralToSquare(x0, y0, x1, y1, x2, y2, x3, y3);
  PerspectiveTransform sToQ =
    PerspectiveTransform::squareToQuadrilateral(x0p, y0p, x1p, y1p, x2p, y2p, x3p, y3p);
  return sToQ.times(qToS);
}

PerspectiveTransform PerspectiveTransform::squareToQuadrilateral(float x0, float y0, float x1, float y1, float x2,
    float y2, float x3, float y3) {
  float dx3 = x0 - x1 + x2 - x3;
  float dy3 = y0 - y1 + y2 - y3;
  if (dx3 == 0.0f && dy3 == 0.0f) {
    PerspectiveTransform result(PerspectiveTransform(x1 - x0, x2 - x1, x0, y1 - y0, y2 - y1, y0, 0.0f,
                                     0.0f, 1.0f));
    return result;
  } else {
    float dx1 = x1 - x2;
    float dx2 = x3 - x2;
    float dy1 = y1 - y2;
    float dy2 = y3 - y2;
    float denominator = dx1 * dy2 - dx2 * dy1;
    float a13 = (dx3 * dy2 - dx2 * dy3) / denominator;
    float a23 = (dx1 * dy3 - dx3 * dy1) / denominator;
    PerspectiveTransform result(PerspectiveTransform(x1 - x0 + a13 * x1, x3 - x0 + a23 * x3, x0, y1 - y0
                                     + a13 * y1, y3 - y0 + a23 * y3, y0, a13, a23, 1.0f));
    return result;
  }
}

PerspectiveTransform PerspectiveTransform::quadrilateralToSquare(float x0, float y0, float x1, float y1, float x2,
    float y2, float x3, float y3) {
  // Here, the adjoint serves as the inverse:
  return squareToQuadrilateral(x0, y0, x1, y1, x2, y2, x3, y3).buildAdjoint();
}

PerspectiveTransform PerspectiveTransform::buildAdjoint() {
  // Adjoint is the transpose of the cofactor matrix:
  PerspectiveTransform result(PerspectiveTransform(a22 * a33 - a23 * a32, a23 * a31 - a21 * a33, a21 * a32
                                   - a22 * a31, a13 * a32 - a12 * a33, a11 * a33 - a13 * a31, a12 * a31 - a11 * a32, a12 * a23 - a13 * a22,
                                   a13 * a21 - a11 * a23, a11 * a22 - a12 * a21));
  return result;
}

PerspectiveTransform PerspectiveTransform::times(PerspectiveTransform other) {
  PerspectiveTransform result(PerspectiveTransform(a11 * other.a11 + a21 * other.a12 + a31 * other.a13,
                                   a11 * other.a21 + a21 * other.a22 + a31 * other.a23, a11 * other.a31 + a21 * other.a32 + a31
                                   * other.a33, a12 * other.a11 + a22 * other.a12 + a32 * other.a13, a12 * other.a21 + a22
                                   * other.a22 + a32 * other.a23, a12 * other.a31 + a22 * other.a32 + a32 * other.a33, a13
                                   * other.a11 + a23 * other.a12 + a33 * other.a13, a13 * other.a21 + a23 * other.a22 + a33
                                   * other.a23, a13 * other.a31 + a23 * other.a32 + a33 * other.a33));
  return result;
}

void PerspectiveTransform::transformPoints(std::vector<float> &points) {
  int max = points.size();
  for (int i = 0; i < max; i += 2) {
    float x = points[i];
    float y = points[i + 1];
    float denominator = a13 * x + a23 * y + a33;
    points[i] = (a11 * x + a21 * y + a31) / denominator;
    points[i + 1] = (a12 * x + a22 * y + a32) / denominator;
  }
}

void PerspectiveTransform::pntProjection(unsigned char *pSrc, int iSrcWidth, int iSrcHeight,
										 unsigned char * pDst, int iDstWidth, int iDstHeight)
{
#pragma omp parallel for
	for(int iX = 0; iX < iSrcWidth; iX++)
	{
		for(int iY = 0; iY < iSrcHeight; iY++)
		{
			float denominator = a13 * iX + a23 * iY + a33;
			int iDstX = (a11 * iX + a21 * iY + a31) / denominator;
			int iDstY = (a12 * iX + a22 * iY + a32) / denominator;

			//仅处理映射在目标区域之内的点
			if(iDstX >= 0 && iDstX < iDstWidth && iDstY >= 0 && iDstY < iDstHeight)
			{
				pDst[iDstX + iDstY * iDstWidth] = pSrc[iX + iSrcWidth * iY];
			}
		}
	}
}

然后用main.cpp引用以上函数，实现透视变换：

#include "Perspective.h"
#include <opencv2/opencv.hpp>
#include <ctime>
#include <iostream>
#include <Windows.h>
#include <WinBase.h>


#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_core249d.lib")  
#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_imgproc249d.lib")  
#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_highgui249d.lib")  
#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_features2d249d.lib")   
#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_contrib249d.lib")  
#pragma comment(lib, "E:\\cv\\opencv\\build\\x86\\vc11\\lib\\opencv_legacy249d.lib")  

using namespace cv;
using namespace std;


int main()
{
	SYSTEMTIME start, end;
	Mat img=imread("E:\\hehe.png");
	
	int img_height = img.rows;
	int img_width = img.cols;
	Mat img_trans = Mat::zeros(img_height,img_width,CV_8UC1);/*CV_8UC3*/
	PerspectiveTransform tansform = PerspectiveTransform::quadrilateralToQuadrilateral
		(
		/*前面的8个坐标是源图的位置*/
		0,0,
		img_width-1,0,
		0,img_height-1,
		img_width-1,img_height-1,

		/*后面8个坐标是目标图的位置*/
		0,0, // top left
		670,200, // top right
		0,300,// bottom left
		350,440
		);

	unsigned char *p1 = img.ptr<unsigned char>(0);
	unsigned char *p2 = new unsigned char[img_width * img_height];

	//由于img是rgb三通道的，而我的变换仅限于1通道灰度图，所以变换之前，先从img中提取一个通道
	for(int i = 0; i < img_width; i++)
	{
		for(int j =0; j < img_height; j++)
		{
			p2[j * img_width + i] = p1[j * img_width * 3 + i * 3];
		}
	}

	GetLocalTime(&start);
	std::cout<<start.wMilliseconds<<"\n";

	tansform.pntProjection(p2, img_width, img_height,
										 img_trans.ptr<unsigned char>(0), img_width, img_height);
	
	GetLocalTime(&end);
	std::cout<<end.wMilliseconds<<"\n";

	delete [] p2;
	std::cin.get();
	imwrite("E:\\trans.png",img_trans);
	return 0;
}

源图：

目标图（结果）：

由于有一些目标图的像素没有被源图赋值（被拉长的区域肯定有对不上源图的点），导致目标图出现一些黑点。从加速的结果看，效果与#pragma omp parallel for的位置关系很大。目前还在总结规律。