柱面投影的C++实现（一）

本文介绍柱面投影具体作法：若要实现柱面投影变换，除了须要获取源图像数据之外，还须要已知相机镜头的一些信息——
（信息A：相机焦距+目标距离）
（信息B：图像变换尺度s，与相机焦距f成正比）
（信息C：相机水平方向的视角，即相机到源图像两边延伸点的射线之间的夹角）

以上三种信息获取其中一种，便可求出平面图像的对应柱面投影。其中，信息A除特殊状况基本不考虑，信息B的正向变换及原理在【注1】中老哥的博客里有详细介绍，本文主要介绍已知相机视场角状况(α = pi / 6)下的实现方法。

首先，已知柱面半径 r = 相机焦距 f，经过l=α·r和源图像水平像素个数就能够获得 s = W / (tan(α))，因此目标柱面的水平方向弧长 L = 2sα, 任一在目标图像的像素点根据这些值来进行坐标变换。

讨论如下四种状况，以图像中心（相机的光心）将水平竖直方向分割图片:
x < W / 2, y < H / 2.
x < W / 2, y > H / 2.
x > W / 2, y < H / 2.
x > W / 2, y > H / 2.

因此根据公式，设源图像像素(x0, y0)，柱面变换后目标图像像素(x1, y1)，则有：
x1 = (L / 2) - s * arctan((W / 2 - x0) / s) , x0 <= (W / 2)
x1 = (L / 2) + s * arctan((x0 - W / 2) / s) , (W / 2) < x0 <= W
一样不难推出：
y1 = (H / 2) - s * (H / 2 - y0) /√ ((W / 2 - x0）^ 2 + s ^ 2) , y0 <= (H / 2)
y1 = (H / 2) + s * (y0 - H / 2) /√ ((W / 2 - x0）^ 2 + s ^ 2) , (H / 2) < y0 <= H

#include "pch.h"
#include <iostream>
#include <opencv2/opencv.hpp>
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

#define PI 3.141592653

using namespace std;
using namespace cv;

int main()
{
	Mat img0, img1;
	img0 = imread("cathedrale.jpg");
	//img1 = img0.clone();
	//img1.setTo(Scalar::all(0));
	int W = img0.cols;							//W为图像的平面宽度
	int H = img0.rows;							//H为图像的柱面高度
	float L;									//目标柱面的宽度
	float a = PI / 6;							//已知相机能看到的半角大小
	float s = W / 2 / (tan(a));					//s为相机的焦距
	int x0, y0;									//图像平面上某一点的横纵坐标（x0, y0）
	float x1, y1;								//目标曲面上对应点的横纵坐标（x1, y1）
	int zsx, zsy;								//整数x1和y1

	L = 2 * s * a;								//求解L的总长度
	
	int springX = (int)L * tan(a) / a +1;
	img1.create(cvSize((int)L,H), img0.type());
	img1.setTo(Scalar::all(0));
	//cout << img0.size() << endl;
	//cout << springX << endl;


	double t0 = (double)getTickCount();
	//分类讨论x0<=(W/2)和(W/2)<x0<=(W)以及y0<=(H/2)和(H/2)<y0<=(H)的状况，以后指令向量化优化
	for (y0 = 1; y0 <= H; y0++)
	{
		for (x0 = 1; x0 <= W; x0++)
		{
			if (x0 <= (W / 2))
				x1 = (L / 2) - s * cvFastArctan((float)(W / 2 - x0), s) * PI / 180;
			else
				if (x0 > (W / 2))
					x1 = (L / 2) + s *cvFastArctan((float)(x0 - (W / 2)), s) * PI / 180;
			if (y0 <= (H / 2))
				y1 = (H / 2) - s * (H / 2 - y0) / sqrt(pow((W / 2) - x0, 2) + pow(s, 2));
			else
				if (y0 > (H / 2))
					y1 = (H / 2) + s * (y0 - H / 2) / sqrt(pow((W / 2) - x0, 2) + pow(s, 2));
			
			//化为整数坐标
			zsx = (int)(x1);
			zsy = (int)(y1);

			//cout << L << endl;
			//cout << zsx << endl;

			if (zsx < W && zsy < H && x0<W && y0<H)
			{
				//cout << zsy << endl;
				//cout << zsx << endl;
				//cout << y0 << endl;
				//cout << x0 << endl;
				img1.at<Vec3b>(zsy, zsx)[0] = img0.at<Vec3b>(y0, x0)[0];
				img1.at<Vec3b>(zsy, zsx)[1] = img0.at<Vec3b>(y0, x0)[1];
				img1.at<Vec3b>(zsy, zsx)[2] = img0.at<Vec3b>(y0, x0)[2];
			}
		}
	}
	double t_proc = ((double)getTickCount() - t0) / getTickFrequency();
	cout << "柱面投影：" << t_proc << endl;
	cout << "变换尺度：" << s << endl;
	
	cout << "柱面变换结束..." << endl;
	imshow("cylindralCathedrale.jpg", img1);
	imwrite("cylindralCathedrale.jpg", img1);
	waitKey();

	/*		向img1赋值
			B0 = *(img0.data + img0.step[0] * y0 + img0.step[1] * x0 + img0.elemSize1() * 0);
			*(img1.data + img1.step[0] * zsy + img1.step[1] * zsx + img1.elemSize1() * 0) = B0;
			G0 = *(img0.data + img0.step[0] * y0 + img0.step[1] * x0 + img0.elemSize1() * 1);
			*(img1.data + img1.step[0] * zsy + img1.step[1] * zsx + img1.elemSize1() * 1) = G0;
			R0 = *(img0.data + img0.step[0] * y0 + img0.step[1] * x0 + img0.elemSize1() * 2);
			*(img1.data + img1.step[0] * zsy + img1.step[1] * zsx + img1.elemSize1() * 2) = R0;
	*/

	/*
	//测试cvFastArctan(分子, 分母),可用，获得的是角度制须要乘以PI/180转化为弧度制算弧长;
	float testsum;
	float testitem = 1.7320508075689 * W;
	testsum = cvFastArctan(W, testitem);
	cout << "结果在这里： " << endl;
	cout << testsum << endl;
	*/
}

注1：这篇文章介绍了已知尺度(s for scale) s=1000状况下的前向映射作法
https://blog.csdn.net/czl389/article/details/54599253
注2：这篇文章介绍了柱面投影变换的计算过程
https://www.cnblogs.com/cheermyang/p/5431170.html
注3：在作完投影变换后会有空洞出现，可是本文因为做者很懒 篇幅问题不作赘述
注4：柱面投影的后向变换下期再见