CN101621701A - 独立于几何校正的任意光滑曲面屏幕多投影仪显示墙色彩校正方法 - Google Patents

Abstract

独立于几何校正的任意光滑曲面屏幕多投影仪显示墙色彩校正方法属计算机***的图像显示处理。步骤是：先利用图像对准算法得到投影仪图像到相机图像间的对应关系；再拍摄四组投影仪图像在RGB三个通道中不同输入下的投影图像；然后利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息；最后联合所有计算结果得到投影仪公共亮度响应区间，通过GPU实时调整输入达到色彩一致性校正效果。本发明脱离几何校正结果，利用自适应的图像对准算法建立投影仪图像到CCD相机图像间像素级的对准关系，可测量投影仪每个像素点的亮度信息，适应任意光滑投影屏幕。对准算法无需任何投影仪和相机以及屏幕曲面参数标定，具有简洁、稳定、自适应和速度快的特点。

Description

独立于几何校正的任意光滑曲面屏幕多投影仪显示墙色彩校正方法

技术领域

本发明涉及计算机***的图像显示处理，特别涉及多通道投影显示墙***的色彩校正方法。

背景技术

随着商业、娱乐和科学研究等领域对超大屏幕高分辨率显示需求的不断增加，多投影显示墙技术成为一个备受关注的研究热点。由于不同投影仪之间的性能差异以及同一台投影仪内部不同区域性能差异，多投影仪显示墙即便在精确的几何对准的情况下也无法投影出视觉上一致性的完整画面。故“色彩校正技术”在构建“多投影仪显示墙”过程中扮演着至关重要的角色。

过去十几年里，众多学者对多投影仪拼接显示的问题做了大量的研究。早期的手工校正工作既繁琐又不精确，已逐渐退出了多投影仪校正领域。光学过渡手段通过在两个投影图像边界引入过渡函数的方式达到图像的光滑过渡；色彩域的匹配手段通过匹配所有投影仪的色彩域来达到全局色彩校正，但理论和实践都显示一些投影仪的色域是无法匹配的；由于投影仪投影图像内部也存在色彩和亮度的变化，一些算法通过“迁就”最差像素点的性质的方法，对每个像素进行变化；这种方法会造成投影图像的对比度严重下降。为此，利用人眼视觉特性或图像统计特性的色彩校正方法相继被提出，这些算法允许图像色彩有一定差异，从而可提高图像的对比度和色彩饱和度。

然而，目前的色彩校正方法不能较好的适应任意光滑的曲面投影屏幕，并且过分依赖几何校正的结果。同时，由于基于线阵CCD相机测量的色彩校正方法需要拍摄大量的图像，其校正速度相对缓慢。

发明内容

本发明的目的是提供一种独立于几何校正的快速投影仪图像色彩校正方法，并使该方法适应于带有任意光滑曲面屏幕的多投影仪显示墙***。

本发明的基本思路是：设计一种细分网格变形技术来对准相机图像和任意光滑曲面屏幕上的投影仪图像，在此基础上提出一种快速、实用的独立于几何校正结果的多投影仪显示墙***色彩校正方法，同时利用加权最小二乘曲线拟合思想和GPU(Graphics Processing Unit图形处理器)像素着色器(Fragment Shader)加速色彩测量和校正计算的过程。

本发明的目的是这样实现的：一种独立于几何校正的任意光滑曲面屏幕多投影仪显示墙色彩校正方法，其特征在于：校正按照以下步骤进行：先利用图像对准算法得到投影仪图像到相机图像间的对应关系，再拍摄四组投影仪图像在RGB(红、绿、蓝)三个通道中不同输入下的投影图像，然后利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息，最后联合所有计算结果得到投影仪公共亮度响应区间，并通过GPU实时调整输入达到色彩一致性校正效果。

所述先利用图像对准算法得到投影仪图像到相机图像间的对应关系的建立方法是：

1)将CCD相机对准投影仪，使投影图像能完整的被相机拍摄到，并固定其位置和朝向；

2)投影特征图像，并使用CCD相机拍摄得到相机图像；

3)利用特征识别的算法建立投影仪图像到CCD相机图像间的稀疏特4点对应关系；

4)利用细分网格加密算法将该稀疏对应点关系加密到稠密的像素级对应关系。

所述拍摄四组投影仪图像在红、绿、蓝RGB三个通道中不同输入下的投影图像的步骤是：

1)令投影仪投影R通道最大亮度图像，并进行多次不同曝光时间的拍摄；

2)令投影仪投影G通道最大亮度图像，并进行多次不同曝光时间的拍摄；

3)令投影仪投影B通道最大亮度图像，并进行多次不同曝光时间的拍摄；

4)令投影仪投影RGB通道最小亮度图像，并进行多次不同曝光时间的拍摄；

5)令投影仪投影带有R通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

6)令投影仪投影带有G通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

7)令投影仪投影带有B通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

8)令投影仪投影带有色块位置特征的图像。

利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息的步骤是：

1)利用图像对准结果从拍摄图像中提取出投影仪有效区域的色彩信息，并使用高动态范围技术HDR计算投影仪每个像素点最大和最小亮度信息；

2)利用特征点定位方式得到投影仪RGB三个通道在不同输入下亮度离散点列，并应用加权最小二乘曲线拟合方法得到投影仪RGB三个通道的亮度响应曲线；

3)联合所有投影仪的计算结果，得到公共亮度响应曲线，并由此计算每个投影仪的亮度衰减模板、平移模板和线性化查表函数。

所述利用GPU的像素着色器达到对投影仪画面逐像素的实时校正计算包括高动态范围成像HDR(High Dynamic Range Imaging)计算、亮度离散点列的读取、加权最小二乘曲线拟合和GPU实时校正方式。

所述利用细分网格加密算法将该稀疏对应点关系加密到稠密的像素级对应关系，是将在矩形域上的任何一个网格都看成是由两个正交方向上的曲线族构成，给出曲线构造形式，在每个网格个点加入一个细分参数，将二维网格拓广为三维网格，沿网格的两个方向分别构造自适应4点插值细分曲线来加密网格，同时加入面点来形成新的网格，保持网格方向的位移平衡和细分后的网格疏密分布特征；对于边界上的点采用镜像的办法将网格内点映射到边界以外，成为虚拟点，消除边界。

所述HDR计算步骤是：

1)读取图像分辨：HDR技术数次曝光；读取相机响应曲线CRF；读取文件名向量；

2)申请权系数矩阵空间：申请亮度矩阵空间；建立加权函数；

3)针对所有的矩阵进行运算，得到运算结果。

所述亮度离散点列读取的步骤是：

1)确定采样数量；确定领域读取半径；

2)根据采样数量生成色块特征图像和RGB三个通道色块亮度图像；

3)投影色块特征图像，并使用自动曝光时间拍摄投影后的图像；

4)分别投影RGB三个通道的色块特征图像，并使用手动控制曝光时间拍摄四种不同曝光时间的投影后图像；

5)利用变形网格恢复拍摄的变形图像；

6)读取投影区域最大、最小亮度；

7)利用HDR技术计算色块亮度矩阵；

8)合成各通道与空间位置无关亮度比率矩阵：

9)利用连通集合的手段定位色块中心：

10)读取色块亮度点列。

所述加权最小二乘曲线拟合是：事先排除噪声和测量误差的影响，利用曲线的光滑性做二阶导数的平滑处理：在最小二乘意义下给出投影仪的连读响应曲线函数，然后构造一新的光滑函数逼近目标；

所述GPU实时校正实现方式是：设待投影图像存储在纹理图形处理器GL_TEXTURE0里，RGB各通道的衰减模板存储在纹理图形处理器GL_TEXTURE1的RGB通道里，平移模板存储在纹理图形处理器GL_TEXTURE1的A通道里，RGB各通道的逆线性化查表存储在纹理图形处理器GL_TEXTURE2的RGB通道里，存储在纹理GL_TEXTURE3的RGB通道里，利用GPU的Fragment Shader实时的对图像进行校正。

本发明具有明显的优点和有益效果：

首先，可以脱离几何校正的结果，利用自适应的图像对准算法建立投影仪图像到CCD相机图像间像素级的对准关系，从而可以测量投影仪每个像素点的亮度信息，并适应任意光滑投影屏幕。该对准方法无需任何投影仪和相机以及屏幕曲面参数标定，并具有简洁、稳定和自适应的特点。

其次，由于本发明采用加权最小二乘曲线拟合思想计算投影仪亮度响应曲线，使对投影仪亮度响应曲线的测量数据量从传统的256张HDR图像降低到了仅需单张HDR图像。在程序实现方面，本发明采用GPU实时校正像素着色器加速校正计算的过程，可达到实时校正的计算速度。

附图说明

图1是本发明主流程图；

图2是基于细分网格加密的图像对准过程示意图。

图3是原始测量点列及其平滑后的曲线图；

图4是本发明在类柱面上的实验结果；

图5是本发明在类柱面和椭球面混合交界处的实验结果。

具体实施方式

按照本发明提出的方法，先利用图像对准算法得到投影仪图像到相机图像间的对应关系，再拍摄四组投影仪图像在RGB三个通道中不同输入下的投影图像，然后利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息，最后联合所有计算结果得到投影仪公共亮度响应区间，并通过GPU实时调整输入达到色彩一致性校正效果。

在图像对准关系的建立中，采用了如前所述的步骤，基于细分网格加密的图像对准过程示意图如图2所示，其细分加密算法由如下描述：

在矩形域上的任何一个网格都可以看成是由两个正交方向上的曲线族构成，

记为 $U_i^k=\left\{\begin{array}{cccc}{p}_{i0}^k& p_{i1}^k& \·\·\·& p_{{\mathrm{in}}^k}^k\end{array}\right\}$ (i＝0，1…m^k) $V_j^k=\left\{\begin{array}{cccc}{p}_{0j}^k& p_{1j}^k& \·\·\·& p_{m^{k}j}^k\end{array}\right\}$ (j＝0，1…n^k)；

其中U_i ^k为第i条u线，V_j ^k为第j条v线， $p_{\mathrm{ij}}^k={(x_{\mathrm{ij}}^k,y_{\mathrm{ij}}^k)}^T\∈R^2$ 为网格格点。整个网格记为M^k，k表示细分次数。在每个网格个点加入一个细分参数ω_ij ^k，将二维网格拓广为三维网格 ${\stackrel{\‾}{M}}^k={\{P_{\mathrm{ij}}^k={(x_{\mathrm{ij}}^k,y_{\mathrm{ij}}^k,{\mathrm{\ω}}_{\mathrm{ij}}^k)}^T\∈R^3\}}_{i,j=1}^{m^k,n^k},$ 沿u/v两个方向分别构造自适应4点插值细分曲线来加密网格，边点构造模式如下：

对于网格中第i条u线，细分曲线递归地定： ${\stackrel{\‾}{U}}_{2i}^{k+1}=\left\{\begin{array}{cccc}{P}_{2i,0}^{k+1}& P_{2i,1}^{k+1}& \·\·\·& P_{2i,{2n}^k}^{k+1}\end{array}\right\}$ (i＝0，1…m^k)：

$\left\{\begin{array}{c}{P}_{2i,2j}^{k+1}=P_{i,j}^k\\ P_{2i,2j+1}^{k+1}=(1+\mathrm{\α})\\ \mathrm{\α}=({\mathrm{\ω}}_{i,j}^k+{\mathrm{\ω}}_{i,j+1}^k)/2\end{array}\left(\frac{P_{i,j}^k+P_{i,j+1}^k}{2}\right)-\mathrm{\α}\left(\frac{P_{i,j-1}^k+P_{i,j+2}^k}{2}\right)\right.---\left(1\right)$

对于网格中第j条v线，细分曲线递归地定义 ${\stackrel{\‾}{V}}_{2j}^{k+1}=\left\{\begin{array}{cccc}{P}_{\mathrm{0,2}j}^{k+1}& P_{\mathrm{1,2}j}^{k+1}& \·\·\·& P_{{2m}^k,2j}^{k+1}\end{array}\right\}$ (j＝0，1…n^k)：

$\left\{\begin{array}{c}{P}_{2i,2j}^{k+1}=P_{i,j}^k\\ P_{2i+1,2j}^{k+1}=(1+\mathrm{\α})\\ \mathrm{\α}=({\mathrm{\ω}}_{i,j}^k+{\mathrm{\ω}}_{i+1,j}^k)/2\end{array}\left(\frac{P_{i,j}^k+P_{i+1,j}^k}{2}\right)-\mathrm{\α}\left(\frac{P_{i-1,j}^k+P_{i+2,j}^k}{2}\right)\right.---\left(2\right)$

为了平衡u/v两个方向的位移，并保持细分后的网格疏密分布特征，加入面点与来形成新的网格：

$P_{2i+\mathrm{1,2}j+1}^{k+1}=\frac{1}{2}(P_{2i+\mathrm{1,2}j}^{k+1}+P_{2i,2j+1}^{k+1}+P_{2i+\mathrm{1,2}j+2}^{k+1}+P_{2i+\mathrm{2,2}j+1}^{k+1})-\frac{1}{4}(P_{2i,2j}^{k+1}+P_{2i,2j+2}^{k+1}+P_{2i+\mathrm{2,2}j}^{k+1}+P_{2i+\mathrm{2,2}j+2}^{k+1})$

由于本发明使用的是开网格结构，故存在边界上的点不能满足细分模式，需要特殊处理。为了保持一致的细分模式，并在边界处同样保持原网格的几何形状和疏密分布特征，本发明采用镜像的办法将网格内点映射到边界以外，成为虚拟点，消除边界。

设网格中某条u曲线的左边界三个顶点为P_i1 ^k，P_i2 ^k，P_i3 ^k，构造虚拟顶点P_i0 ^k，使其为P_i3 ^k关于P_i1 ^k，P_i2 ^k中垂线的镜像，则 $P_{i0}^k=P_{i3}^k+(P_{i1}^k-P_{i2}^k)+\frac{2\left|(P_{i2}^k-P_{i3}^k,P_{i1}^k-P_{i2}^k)\right|}{{\left({\left|\right|P_{i1}^k-P_{i2}^k\left|\right|}_2\right)}^2}(P_{i1}^k-P_{i2}^k)$

其中(*，*)表示向量内积，‖*‖₂表示向量的2范数，v线边界的构造式类似。

细分网格的自适应性由初始细分参数的确定来体现，设投影仪图像到相机图像间的稀疏特征点对应关系中某初始网格为 $M^0={\left\{p_{\mathrm{ij}}^0\right\}}_{i,j=1}^{m^0,n^0},$ 可根据网格自身的性质来确定自适应细分参数矩阵 ${\mathrm{\ω}}^0={\left({\mathrm{\ω}}_{\mathrm{ij}}^0\right)}_{m^0\×n^0},$ 以便得到初始拓广后的网格 ${\stackrel{\‾}{M}}^0={\{P_{\mathrm{ij}}^0={(x_{\mathrm{ij}}^0,y_{\mathrm{ij}}^0,{\mathrm{\ω}}_{\mathrm{ij}}^0)}^T\}}_{i,j=1}^{m^0,n^0}.$

设M^-1为对初始网格M⁰进行间隔采样后的粗网格，期望通过优化的方式寻找合适的细分参数矩阵ω^-1，使得由M^-1按照(1)(2)式加密一次的结果尽量的逼近M⁰。

设 $S\left(\mathrm{\ω}\right):M^{-1}\stackrel{\mathrm{\ω}}{\→}{\{p_{\mathrm{ij}}^{\′}={(x_{\mathrm{ij}}^{\′},y_{\mathrm{ij}}^{\′})}^T\}}_{i,j=1}^{m^0,n^0}$ 即为M^-1在细分参数ω下按照(1)(2)式加密一次的结果，则可寻找最优的ω作为细分参数矩阵ω^-1：

${\mathrm{\ω}}^{-1}=\arg \underset{\mathrm{\ω}}{\min }\underset{i=1}{\overset{m^0}{\mathrm{\Σ}}}\underset{j=1}{\overset{n^0}{\mathrm{\Σ}}}{\left|\right|p_{\mathrm{ij}}^0-p_{\mathrm{ij}}^{\′}\left|\right|}^2---\left(3\right)$

由于初始网格格点数本身并不多，间隔采样后格点点数更少，故(3)式是一个低维最佳平方逼近问题，可采用一般的数值优化算法求解。

将(3)式求得的ω^-1与M^-1结合得到拓广后的粗网格 ${\stackrel{\‾}{M}}^{-1}={\{P_{\mathrm{ij}}^{-1}={(x_{\mathrm{ij}}^{-1},y_{\mathrm{ij}}^{-1},{\mathrm{\ω}}_{\mathrm{ij}}^{-1})}^T\}}_{i,j=1}^{m^{-1},n^{-1}}$ 并按照(1)(2)式加密一次即可得到拓广后的初始网格M⁰。

在利用GPU的像素着色器达到对投影仪画面逐像素的实时校正计算主要包括了HDR计算、亮度离散点列的读取、加权最小二乘曲线拟合和GPU实时校正方式。

HDR计算算法的详细描述如下

1、读取图像分辨：(M N)，HDR技术曝光次数：K；

读取相机响应曲线(CRF)：f_camera[i]i＝0，1…255；

读取文件名向量：fileName[i]i＝1，2…K。

2、申请权系数矩阵空间：

$\mathrm{new}{W}_r={\left(w_{\mathrm{ij}}^r\right)}_{M\×N},$ $\mathrm{new}{W}_g={\left(w_{\mathrm{ij}}^g\right)}_{M\×N}{\mathrm{newW}}_b={\left(w_{\mathrm{ij}}^b\right)}_{M\×N};$ 申请亮度矩阵空间：

$\mathrm{new}{L}_r={\left(l_{\mathrm{ij}}^r\right)}_{M\×N},$ $\mathrm{new}{L}_g={\left(l_{\mathrm{ij}}^g\right)}_{M\×N},$ $\mathrm{new}{L}_b={\left(l_{\mathrm{ij}}^b\right)}_{M\×N};$ new L＝(l_ij)_M×N；

建立加权函数： $f_{\mathrm{weight}}\left(x\right)=\left\{\begin{array}{cc}x/128& x<128\\ (255-x)/128& 128\&le;x<255\end{array}\right..$

3、for i＝1 to K

Image＝readImage(filename[i])；

exposeTime＝readExposeTime(filename[i])；

Image＝subdivisionMesh(Image)；

exposeTime＝ln(exposeTime)；

$W_r=W_r\&CirclePlus;f_{\mathrm{weight}}[\mathrm{Image}.R];$

$W_g=W_g\&CirclePlus;f_{\mathrm{weight}}[\mathrm{Image}.G];$

$W_b=W_b\&CirclePlus;f_{\mathrm{weight}}[\mathrm{Image}.B];$

end for

4、 $L_r=L_r\&CircleTimes;{\left(W_r\right)}^{-1};$ $L_g=L_g\&CircleTimes;{\left(W_g\right)}^{-1};$ $L_b=L_b\&CircleTimes;{\left(W_b\right)}^{-1}.$

5、 $L=0.299\&CenterDot;\exp \left(L_r\right)\&CirclePlus;0.587\&CenterDot;\exp \left(L_g\right)\&CirclePlus;0.114\&CenterDot;\exp \left(L_b\right),$ 输出L。

上述计算过程中的所有针对矩阵的运算含义为：对任意矩A＝(a_ij)_M×N和B＝(b_ij)_M×N，

A^-1＝(1/a_ij)_M×N，f(A)＝[f(a_ij)]_M×N，k·A＝(k·a_ij)_M×N， $A\&CircleTimes;B={(a_{\mathrm{ij}}\&CenterDot;b_{\mathrm{ij}})}_{M\&times;N},$ $A\&CirclePlus;B={(a_{\mathrm{ij}}+b_{\mathrm{ij}})}_{M\&times;N}.$

亮度离散点列的读取：

1、确定采样数量：K² 2≤K≤16；确定领域读取半径：δ；

2、根据K生成色块特征图像和RGB三个通道色块亮度图像；

3、投影色块特征图像，并使用自动曝光时间拍摄投影后的图像；

4、分别投影RGB三个通道的色块特征图像，并使用手动控制曝光时间拍摄10中不同曝光时间的投影后图像；

5、利用变形网格恢复拍摄的变形图像；

6、读取投影区域最大、最小亮度：L_max ^r，L_max ^g，L_max ^b，L_min ^r，L_min ^g，L_min ^b；

7、利用HDR技术计算色块亮度矩阵：L_block ^r，L_block ^g，L_block ^b；

8、合成各通道与空间位置无关亮度比率矩阵： ${\mathrm{\&alpha;}}_r=(L_{\mathrm{block}}^r-L_{\min }^r)\&CircleTimes;{(L_{\max }^r-L_{\min }^r)}^{-1};$ ${\mathrm{\&alpha;}}_g=(L_{\mathrm{block}}^g-L_{\min }^g)\&CircleTimes;{(L_{\max }^g-L_{\min }^g)}^{-1};$ ${\mathrm{\&alpha;}}_b=(L_{\mathrm{block}}^b-L_{\min }^b)\&CircleTimes;{(L_{\max }^b-L_{\min }^b)}^{-1};$

9、利用连通集合的手段定位色块中心：

10、读取色块亮度点列： $r={\{\underset{I=x_i-\mathrm{\&delta;}}{\overset{x_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}\underset{m=y_i-\mathrm{\&delta;}}{\overset{y_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_r^{\mathrm{Im}}/{\mathrm{\&delta;}}^2\}}_{i=1}^{K^2}$ $g={\{\underset{I=x_i-\mathrm{\&delta;}}{\overset{x_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}\underset{m=y_i-\mathrm{\&delta;}}{\overset{y_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_g^{\mathrm{Im}}/{\mathrm{\&delta;}}^2\}}_{i=1}^{K^2}$ $b={\{\underset{I=x_i-\mathrm{\&delta;}}{\overset{x_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}\underset{m=y_i-\mathrm{\&delta;}}{\overset{y_i+\mathrm{\&delta;}}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_b^{\mathrm{Im}}/{\mathrm{\&delta;}}^2\}}_{i=1}^{K^2}.$

加权最小二乘曲线拟合其原始测量点列及其平滑后的曲线如图3所示。具体做法是：

设 $P={\{p_i=(x_i,y_i)\}}_{i=0}^{255}$ 为实际测量的离散点列，由于各种测量误差和噪声的影响，P并不是一个光滑的点列。投影仪的亮度响应曲线是由投影仪内部的电子线路和DLP光敏特性所决定的，一般情况下是一光滑的连续函数y＝f(x)。如果想通过测量点列P来推测出f(x)，必须事先排除噪声和测量误差的影响，即需要对P做平滑处理。本发明利用曲线的光滑性(二阶导数能量最小)，在最小二乘意义下给出一种离散点列拟合的方法。设y＝f(x)为投影仪的连读响应曲线，我们构造一新的光滑函数y＝g(x)逼近f(x)，逼近的标准为满足如下目标函数：

$\min \mathrm{\&Theta;}={\left|\right|f\left(x\right)-g\left(x\right)\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|g^{\&prime;\&prime;}\left(x\right)\left|\right|}_2^2.$

GPU实时校正实现方式是：设待投影图像存储在纹理GL_TEXTURE0里，RGB各通道的衰减模板存储在纹理图形处理器GL_TEXTURE1的RGB通道里，平移模板存储在纹理图形处理器GL_TEXTURE1的A通道里，RGB各通道的逆线性化查表存储在纹理图形处理器GL_TEXTURE2的RGB通道里，存储在纹理图形处理器GL_TEXTURE3的RGB通道里。

像素着色器伪代码为：

    uniform sampler2DRectT0；//GL_TEXTURE0

    uniform sampler2DRectT1；//GL_TEXTURE1

    uniform sampler1DRectT2；//GL_TEXTURE2

    uniform sampler1DRectT3；//GL_TEXTURE3

    uniform float scale；

    void main()

    {

       vec2 TexCoord＝gl_TexCoord[0].st；

       vec4 PInput＝texture2D(T0，TexCoord)；

       vec4 mask＝texture2D(T1，TexCoord)；

       vec4 lookup；

       vec4 lookupInverse；

       vec3 Poutput；

       PInput＝RGB2HSV(PInput)；

       PInput[1]＝PInput[1]*scale；//提升色彩饱和度

       PInput＝HSV2RGB(PInput)；

       POutput.r＝texture1D(T2，PInput.r).r；//查表

       POutput.g＝texture1D(T2，PInput.g).g；

       POutput.b＝texture1D(T2，PInput.b).b；

    POutput＝vec3(mask.rgb)*POutput+vec3(mask.a，mask.a，mask.a)；//校

正

       POutput.r＝texture1D(T3，POutput.r).r；//查表

       POutput.g＝texture1D(T3，POutput.g).g；

       POutput.b＝texture1D(T3，POutput.b).b；

       gl_FragColor＝vec4(POutput，1.0)；//返回校正后颜色

    }

在进行实际校正时，详细操作步骤如下：

一、测试图像生成

参见图1。生成屏幕分辨率的特征图像，该图像背景为黑色，四周特征为白色方块，中间特征为白色圆块。生成屏幕分辨的红色(R＝255；G＝0；B＝0)、绿色(R＝0；G＝255；B＝0)和蓝色(R＝0；G＝0；B＝255)三幅最大亮度图像。生成屏幕分辨率的红色色块图像，该图像背景为黑色，前景为B、G值为0，R值从0到255均匀递增的100个红色方块，并均分布在图像中。生成屏幕分辨率的绿色色块图像，该图像背景为黑色，前景为R、B值为0，G值从0到255均匀递增的100个绿色方块，并均分布在图像中。生成屏幕分辨率的蓝色色块图像，该图像背景为黑色，前景为R、G值为0，B值从0到255均匀递增的100个蓝色方块，并均分布在图像中。生成屏幕分辨率的色块特征图像，该图像背景为黑色，前景为任意与背景具有较大区分度的均匀分布的100个色块。生成屏幕分辨率的全黑色图像。

二、图像对准步骤

通过远程网络控制软件VNC控制投影仪投影特征图像，并通过CCD相机SDK包控制相机在自动曝光时间下拍摄投影图像。利用图像分割、角点识别和重心提取方法定位相机图像中特征点对应位置，从而建立相机图像到投影仪图像间的稀疏特征点对应关系。通过本发明给出的细分加密算法将该稀疏对应关系加密到像素级，并保存计算结果。在后续的拍摄过程中，保证CCD相机的物理位置、朝向和内部参数中除曝光时间外的其他值不再变化。

三、图像拍摄

远程控制投影仪投影红色背景图像，并控制相机在10组递增的曝光时间下拍摄投影画面，并利用像素级对应关系从相机图像中提取出投影图像所在区域。以同样的方式拍摄绿色、蓝色的背景图像以及红色、蓝色、绿色的色块图像和黑色背景图像。控制投影仪投射色块特征图像，并在CCD相机自动曝光时间下拍摄其投影图像。

四、投影仪亮度信息计算

利用在技术方案中提到的HDR计算方法计算拍摄的背景亮度和色块亮度信息。提取色块特征图像中色块中心位置，并提取该位置上最大和最小亮度以及色块亮度信息，从而得到投影仪三个通道亮度响应函数的离散点列。利用加权最小二乘曲线拟合出投影仪亮度响应曲线。在处理完所有投影仪亮度信息计算后，寻找全部投影仪像素点最大亮度的最小值作为公共亮度响应区间的上限，寻找全部投影仪像素点最小亮度的最大值作为公共亮度响应区间的下限，并通过逐像素的计算得到每个投影仪色彩校正的衰减、平移模板。

五、实时校正

将每个投影仪色彩校正的衰减、平移模板和线性化查表函数以纹理的形式载入各投影仪控制计算机的显存，启动像素着色器将每一帧待投影的图像进行逐像素的色彩校正计算。

本发明已在到四川大学西区展厅柱面360度全景显示墙***和四川川大智胜股份有限公司的带穹顶的机场虚拟搭台仿真***中进行了实施测试。由于加工技术的限制和自身重力的影响，其投影屏幕只是近似柱面和椭球面，并存在一定程度的变形。所有计算都使用单台PC机(AMD Athlon 64 X2 Dual Core Processor4400+；2048M memory；NVIDIA GeForce 7800 GTX GPU)上进行计算，使用OLYMPUSSP350 Digital Camera以及相应的SDK开发包作为色彩测量设备，使用VNC软件实现远程控制投影仪显示画面的功能。本例中，拍摄一幅最高精度的图像需要15秒左右，对每台投影仪共拍摄72张，对每个投影仪的色彩校正全过程，包括拍摄、亮度计算和校正计算，仅需要20分钟左右的时间。

图4为本发明在四川大学西区展厅柱面360度全景显示墙***中的应用实例，其投影仪为EIKI系列投影仪。图5为本发明在四川川大智胜股份有限公司的带穹顶的机场虚拟搭台仿真***中的应用实例，其投影仪为Panasonic系列投影仪。

Claims (8)

1、一种独立于几何校正的任意光滑曲面屏幕多投影仪显示墙色彩校正方法，其特征在于：校正步骤：先利用图像对准算法得到投影仪图像到相机图像间的对应关系，再拍摄四组投影仪图像在RGB三个通道中不同输入下的投影图像，然后利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息，最后联合所有计算结果得到投影仪公共亮度响应区间，并通过GPU实时调整输入达到色彩一致性校正效果。

2、如权利要求1所述的校正方法，其特征在于：所述先利用图像对准算法得到投影仪图像到相机图像间的对应关系的建立方法是：

1)将CCD相机对准投影仪，使投影图像能完整的被相机拍摄到，并固定其位置和朝向；

2)投影特征图像，并使用CCD相机拍摄得到相机图像；

3)利用特征识别的算法建立投影仪图像到CCD相机图像间的稀疏特征4点对应关系；

4)利用细分网格加密算法将该稀疏对应点关系加密到稠密的像素级对应关系。

3、如权利要求1所述的校正方法，其特征在于：所述拍摄四组投影仪图像在RGB三个通道中不同输入下的投影图像的步骤是：

1)令投影仪投影R通道最大亮度图像，并进行多次不同曝光时间的拍摄；

2)令投影仪投影G通道最大亮度图像，并进行多次不同曝光时间的拍摄；

3)令投影仪投影B通道最大亮度图像，并进行多次不同曝光时间的拍摄；

4)令投影仪投影RGB通道最小亮度图像，并进行多次不同曝光时间的拍摄；

5)令投影仪投影带有R通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

6)令投影仪投影带有G通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

7)令投影仪投影带有B通道不同亮度值的色块图像，并进行多次不同曝光时间的拍摄；

8)令投影仪投影带有色块位置特征的图像，并进行一次拍摄。

4、如权利要求1所述的校正方法，其特征在于：利用图像对准关系和拍摄结果计算投影仪每个像素点的色彩信息的步骤是：

1)利用图像对准结果从拍摄图像中提取出投影仪有效区域的色彩信息，并使用高动态范围技术HDR计算投影仪每个像素点最大和最小亮度信息；

2)利用特征点定位方式得到投影仪RGB三个通道在不同输入下亮度离散点列，并应用加权最小二乘曲线拟合方法得到投影仪RGB三个通道的亮度响应曲线；

3)联合所有投影仪的计算结果，得到公共亮度响应曲线，并由此计算每个投影仪的亮度衰减模板、平移模板和线性化查表函数。

5、如权利要求1所述的校正方法，其特征在于：所述利用GPU的像素着色器达到对投影仪画面逐像素的实时校正计算包括HDR计算、亮度离散点列的读取、加权最小二乘曲线拟合和GPU实时校正方式。

6、如权利要求2所述的校正方法，其特征在于：所述利用细分网格加密算法将该稀疏对应点关系加密到稠密的像素级对应关系，是将在矩形域上的任何一个网格都看成是由两个正交方向上的曲线族构成，给出曲线构造形式，在每个网格个点加入一个细分参数，将二维网格拓广为三维网格，沿网格的两个方向分别构造自适应4点插值细分曲线来加密网格，同时加入面点来形成新的网格，保持网格两个方向的位移平衡和细分后的网格疏密分布特征；对于边界上的点采用镜像的办法将网格内点映射到边界以外，成为虚拟点，消除边界。

7、如权利要求6所述的校正方法，其特征在于：所述细分加密算法的具体描述：

在矩形域上的任何一个网格都可以看成是由两个正交方向上的曲线族构成，

记为 $U_i^k=\left\{\begin{array}{cccc}{p}_{i0}^k& p_{i1}^k& ...& p_{{\mathrm{in}}^k}^k\end{array}\right\}(i=\mathrm{0,1}...m^k)$ $V_j^k=\left\{\begin{array}{cccc}{p}_{0j}^k& p_{1j}^k& ...& p_{m^{k}j}^k\end{array}\right\}(j=\mathrm{0,1}...n^k)$

其中U_i ^k为第i条u线，V_j ^k为第j条v线， $p_{\mathrm{ij}}^k={(x_{\mathrm{ij}}^k,y_{\mathrm{ij}}^k)}^T\&Element;R^2$ 为网格格点，整个网格记为M^k，k表示细分次数。在每个网格个点加入一个细分参数ω_ij ^k，将二维网格拓广为三维网格 ${\stackrel{\&OverBar;}{M}}^k={\{P_{\mathrm{ij}}^k={(x_{\mathrm{ij}}^k,y_{\mathrm{ij}}^k,{\mathrm{\&omega;}}_{\mathrm{ij}}^k)}^T\&Element;R^3\}}_{i,j=1}^{m^k,n^k},$ 沿u/v两个方向分别构造自适应4点插值细分曲线来加密网格，边点构造模式如下：

对于网格中第i条u线，细分曲线递归地定义 ${\stackrel{\&OverBar;}{U}}_{2i}^{k+1}=\left\{\begin{array}{cccc}{P}_{2i,0}^{k+1}& P_{2i,1}^{k+1}& ...& P_{2i,{2n}^k}^{k+1}\end{array}\right\}(i=\mathrm{0,1}...m^k):$

$\left\{\begin{array}{c}{P}_{2i,2j}^{k+1}=P_{i,j}^k\\ P_{2i,2j+1}^{k+1}=(1+\mathrm{\&alpha;})\left(\frac{P_{i,j}^k+P_{i,j+1}^k}{2}\right)-\mathrm{\&alpha;}\left(\frac{P_{i,j-1}^k+P_{i,j+2}^k}{2}\right)\\ \mathrm{\&alpha;}=({\mathrm{\&omega;}}_{i,j}^k+{\mathrm{\&omega;}}_{i,j+1}^k)/2\end{array}\right.---\left(1\right)$

对于网格中第j条v线，细分曲线递归地定义 ${\stackrel{\&OverBar;}{V}}_{2j}^{k+1}=\left\{\begin{array}{cccc}{P}_{\mathrm{0,2}j}^{k+1}& P_{\mathrm{1,2}j}^{k+1}& ...& P_{2m^k,2j}^{k+1}\end{array}\right\}(j=\mathrm{0,1}...n^k):$

$\left\{\begin{array}{c}{P}_{2i,2j}^{k+1}=P_{i,j}^k\\ P_{2i+1,2j}^{k+1}=(1+\mathrm{\&alpha;})\left(\frac{P_{i,j}^k+P_{i+1,j}^k}{2}\right)-\mathrm{\&alpha;}\left(\frac{P_{i-1,j}^k+P_{i+2,j}^k}{2}\right)\\ \mathrm{\&alpha;}=({\mathrm{\&omega;}}_{i,j}^k+{\mathrm{\&omega;}}_{i+1,j}^k)/2\end{array}\right.---\left(2\right)$

为了平衡u/v两个方向的位移，并保持细分后的网格疏密分布特征，加入面点与来形成新的网格：

$P_{2i+\mathrm{1,2}j+1}^{k+1}=\frac{1}{2}(P_{2i+\mathrm{1,2}j}^{k+1}+P_{2i,2j+1}^{k+1}+P_{2i+\mathrm{1,2}j+2}^{k+1}+P_{2i+\mathrm{2,2}j+1}^{k+1})-\frac{1}{4}(P_{2i,2j}^{k+1}+P_{2i,2j+2}^{k+1}+P_{2i+\mathrm{2,2}j}^{k+1}+P_{2i+\mathrm{2,2}j+2}^{k+1});$

所述对于边界上的点采用镜像的办法将网格内点映射到边界以外，成为虚拟点，消除边界是：设网格中某条u曲线的左边界三个顶点为P_i1 ^k，P_i2 ^k，P_3i ^k，构造虚拟顶点P_i0 ^k，使其为P_i3 ^k关于P_i1 ^k，P_i2 ^k中垂线的镜像，则 $P_{i0}^k=P_{i3}^k+(P_{i1}^k-P_{i2}^k)+\frac{2\left|(P_{i2}^k-P_{i3}^k,P_{i1}^k-P_{i2}^k)\right|}{{\left({\left|\right|P_{i1}^k-P_{i2}^k\left|\right|}_2\right)}^2}(P_{i1}^k-P_{i2}^k)$

其中(*，*)表示向量内积，‖*‖₂表示向量的2范数，v线边界的构造式类似；

细分网格的自适应性由初始细分参数的确定来体现，设投影仪图像到相机图像间的稀疏特征点对应关系中某初始网格为 $M^0={\left\{p_{\mathrm{ij}}^0\right\}}_{i,j=1}^{m^0,n^0},$ 根据网格自身的性质来确定自适应细分参数矩阵 ${\mathrm{\&omega;}}^0={\left({\mathrm{\&omega;}}_{\mathrm{ij}}^0\right)}_{m^0\&times;n^0},$ 得到初始拓广后的网格 ${\stackrel{\&OverBar;}{M}}^0={\{P_{\mathrm{ij}}^0={(x_{\mathrm{ij}}^0,y_{\mathrm{ij}}^0,{\mathrm{\&omega;}}_{\mathrm{ij}}^0)}^T\}}_{i,j=1}^{m^0,n^0}$

设M^-1为对初始网格M⁰进行间隔采样后的粗网格，期望通过优化的方式寻找合适的细分参数矩阵ω^-1，使得由M^-1按照(1)(2)式加密一次的结果尽量的逼近M⁰；

设 $S\left(\mathrm{\&omega;}\right):M^{-1}\stackrel{\mathrm{\&omega;}}{\&RightArrow;}{\{p_{\mathrm{ij}}^{\&prime;}={(x_{\mathrm{ij}}^{\&prime;},y_{\mathrm{ij}}^{\&prime;})}^T\}}_{i,j=1}^{m^0,n^0}$ 即为M^-1在细分参数ω下按照(1)(2)式加密一次的结果，则可寻找最优的ω作为细分参数矩阵ω^-1：

${\mathrm{\&omega;}}^{-1}=\arg \underset{\mathrm{\&omega;}}{\min }\underset{i=1}{\overset{m^0}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n^0}{\mathrm{\&Sigma;}}}{\left|\right|p_{\mathrm{ij}}^0-p_{\mathrm{ij}}^{\&prime;}\left|\right|}^2---\left(3\right)$

对(3)式采用数值优化算法求解：

将(3)式求得的ω^-1与M^-1结合得到拓广后的粗网格 ${\stackrel{\&OverBar;}{M}}^{-1}={\{P_{\mathrm{ij}}^{-1}={(x_{\mathrm{ij}}^{-1},y_{\mathrm{ij}}^{-1},{\mathrm{\&omega;}}_{\mathrm{ij}}^{-1})}^T\}}_{i,j=1}^{m^{-1},n^{-1}}$ 并按照(1)(2)式加密一次即可得到拓广后的初始网格M⁰。

8、如权利要求5所述的校正方法，其特征在于：

所述HDR计算步骤是：

1)读取图像分辨：HDR技术数次曝光；读取相机响应曲线CRF；读取文件名向量；

2)申请权系数矩阵空间：申请亮度矩阵空间；建立加权函数；

3)针对所有的矩阵进行运算，得到运算结果；

所述亮度离散点列读取的步骤是：

1)确定采样数量；确定领域读取半径；

2)根据采样数量生成色块特征图像和RGB三个通道色块亮度图像；

3)投影色块特征图像，并使用自动曝光时间拍摄投影后的图像；

4)分别投影RGB三个通道的色块特征图像，并使用手动控制曝光时间拍摄四种不同曝光时间的投影后图像；

5)利用变形网格恢复拍摄的变形图像；

6)读取投影区域最大、最小亮度；

7)利用HDR技术计算色块亮度矩阵；

8)合成各通道与空间位置无关亮度比率矩阵：

9)利用连通集合的手段定位色块中心：

10)读取色块亮度点列；

所述加权最小二乘曲线拟合是事先排除噪声和测量误差的影响，利用曲线的光滑性做二阶导数的平滑处理：在最小二乘意义下给出投影仪的连读响应曲线函数，然后构造一新的光滑函数逼近目标；

所述GPU实时校正实现方式是：设待投影图像存储在纹理图像处理器GL_TEXTURE0里，RGB各通道的衰减模板存储在纹理图像处理器GL_TEXTURE1的RGB通道里，平移模板存储在纹理图像处理器GL_TEXTURE1的A通道里，RGB各通道的逆线性化查表存储在纹理图像处理器GL_TEXTURE2的RGB通道里，存储在纹理图像处理器GL_TEXTURE3的RGB通道里，利用GPU的Fragment Shader实时的对图像进行校正。

Images