本发明涉及一种合成高动态范围图像的方法，属于数字图像处理技术领域。
背景技术：
随着人们生活水平的提高，人们对于图像的质量的要求也越来越高，普通的图像即低动态范围图像由于其深度信息表达范围的局限性，已经不足以满足人们的需求。比如现实场景中的亮度动态范围一般为105:1，通常人眼在对环境进行观察时，可以察觉到亮度为10-3cd/m2的星光，也能感受到亮度为105cd/m2的阳光，且人眼在很暗或很亮的情况下都能察觉到现实场景中的各种细节。而一般的图像没有这么大的表现范围。于是人们开始寻找一种可以满足人类更多需求，可以展现更多的场景细节，可以更好地满足人类视觉需求的一种图像表达方式。高动态范围图像(high dynamic range image)应运而生，它是一种可以表示实际场景中亮度大范围变化的图像类型，因此，可以更好地同时表现场景中亮区和暗区的光学特性，高动态范围图像还具备像素值正比于场景中对应点的实际亮度值的特性。在图像处理的很多应用领域，因为普通图像仅仅记录了色调相对关系，不记录场景的具体亮度，缺乏理想的对比度和亮度范围，通常要使用高动态范围图像。目前已经有很多方法合成高动态范围图像，如Debevec与Malik算法、Nayar算法，以及它们的各种改算。
合成高动态范围图像的关键步骤是相机响应曲线的恢复，目前存在的计算相机响应曲线方法，多数情况是通过使用一组曝光度不同的图像并自动选取若干个点的方法，仅仅利用图像的RGB值和曝光时间，来生成相机响应曲线。虽然这种方法可以快速、高效的合成高动态范围图像，但是由于所选取点的真实亮度值未知，而是通过一些理论公式计算得来的，不可能适用于全部的像素点，因此会产生很多的误差，这大大降低了相机响应曲线的准确性。
技术实现要素：
本发明的目的在于提供一种合成高动态范围图像的方法，为了快速、高精度地生成高动态范围图像，其通过拍摄一组曝光度不同的干净天空的图片，使用Lumincept软件模拟该图片所对应的亮度信息模型，该模型记录了所拍摄的这一组曝光度不同的干净天空所对应的真实光照信息，图片中的RGB值和模型对应的真实亮度值，构成了一组非线性关系，即为要恢复的相机响应曲线，通过使用CIE模型中的三个可以反应人类眼睛敏感度的参数来计算合成高动态范围图像的三个通道的值。
本发明的技术方案是这样实现的：一种合成高动态范围图像的方法，其特征在于：在计算机中装有Lumincept9.31.1软件，按以下步骤完成：
步骤1、选择合适的场景，场景中只包括一个塔尖形状的物体，塔尖形状的物体远高于场景中的其他物体,场景中还包括没有云彩的天空，这样的场景记为Scene；
步骤2、获取塔尖形状的物体的底部到顶部之间的距离，用变量Bulid_H表示，距离的单位为米；获取塔尖形状的物体在WGS84坐标系下的经度和纬度，分别用变量Location_longitude1、Location_latitude1表示。
步骤3、调节相机的位置和姿态，使Scene中塔尖形状的物体的最高处位于相机的取景器的中心位置。始终保持相机的位置和姿态不变，获取相机在WGS84坐标系下的经度和纬度，分别用变量Location_longitude2、Location_latitude2表示；由小到大的改变相机的曝光时间A次，相机的其他参数不变，通过相机获取每个曝光时间对应的一张Scene照片，第I张照片的曝光时间设为S(I)，曝光时间的单位为秒，第I张照片用Image(I)表示，I=1,2,...,A；记录拍摄第一张照片时的当地时间，用变量Time表示。
步骤4、查看谷歌地图，分别找到相机和Scene中塔尖形状的物体在谷歌地图中的位置，获得他们之间的距离，用变量Distance表示，利用Distance和Build_H的值，在Lumicept软件中模拟相机的位置和方向。
步骤5、Lumicept软件中使用CIE(Commission Internationale de L’Eclairage )标准天空模型，根据Time、Bulid_H、Location_longitude1、Location_latitude1和Location_longitude2、Location_latitude2，在Lumicept软件中复现Image(1)中的场景，将复现的场景记为Simulated_model。
步骤6、设，寻找Image(k)中所有的天空区域，用Sky表示，将Sky中所有像素点按照从左到右，从上到下的顺序构成矩阵Z；将矩阵Z的行数设为m,列数设为n，Z(c,d)=(R(c,d),G(c,d),B(c,d))为矩阵Z的第c行第d列元素，其中R(c,d)、G(c,d)、B(c,d)分别为Z(c,d)的RGB值，c=1,2,…,m；d=1,2,…,n。
步骤7、构造矩阵Gy，设元素Gy（c,d）是矩阵Gy的第c行第d列元素，Gy(c,d)=B(c,d)，c=1,2,...,m；d=1,2,...,n。
设矩阵Gy中元素的最大值用maxGy表示，矩阵Gy中元素的最小值用minGy表示。
步骤8、根据矩阵Gy中的元素的最大值和最小值来确定采样点的个数N，只需满足方程：
步骤9、构造向量Lrange，Lrange(t)为向量Lrange中的第t个元素，
。
步骤10、构造向量RandNumber，RandNumber(e)为RandNumber的第e个元素，且RandNumber(e)为整数，;
步骤11、设在矩阵Gy中与RandNumber(e)相同元素的个数为size(e)，e=1,2….,N；构造矩阵Hp，设元素Hp(i,j)是矩阵Hp的第i行第j列元素，i=1,2,…,size(p)；j=1,2；p=1,2,…,N；构造向量w，w(q)为向量w的第q个元素，w(q)=1，q=1,2,…,size(p)，p=1,2,…,N。
步骤12、按照从上到下，从左到右的顺序，如果Gy(c,d)等于 RandNumber(e)，那么Hp(w(p),1)=c；Hp(w(p),2)=d，w(p)=w(p)+1；构造向量random (i)为向量random的第i个值，random(i)的值是大于1并且小于的任意一个整数，c=1,2,…,m；d=1,2,…,n，i=1,2,…,size(p)，e=p=1,2,…,N。
步骤13、构造矩阵RandTest，其中RandTest(i,j)是矩阵RandTest中的第i行第j列元素，i=1,2,…,N；j=1,2；RandTest(i,1)=Hp(random(i),1)，RandTest(i,2)=Hp(random(i),2)，i=1,2,…,N。
步骤14、将Image(L)上所有像素点按照从左到右，从上到下的顺序构成矩阵ZL，将矩阵ZL的行数设为r,列数设为s,设元素ZL(i,j)=(R_L(i,j),G_L(i,j),B_L(i,j))，ZL(i,j)是矩阵ZL的第i行第j列元素，其中R_L(i,j)、G_L(i,j)、B_L(i,j)分别为ZL(i,j)的RGB值i=1,2,…,r；j=1,2,…,s；L=1,2,…,A。
步骤15、构造一个向量U_R，其中U_R(X)为向量U_R中第X个元素；U_R(X)=R_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤16、构造一个向量U_G，其中U_G(X)为向量U_G中第X个元素;当U_G(X)=G_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤17、构造一个向量U_B，其中U_B (X)为向量U_B中第X个元素;当U_B(X)=B_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤18、构造向量Luminance_R，Luminance_R(Y)为向量Luminance_R中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1),RandTest(i,2))的像素点R通道Luminance值，记为LuminanceR(RandTest(i,1),RandTest(i,2))；则Luminance_R(Y)=S(I)*LuminanceR(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤19、构造向量Luminance_G，Luminance_G(Y)为向量Luminance_G中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1), RandTest(i,2))的像素点的G通道Luminance值，记为LuminanceG(RandTest(i,1),RandTest(i,2))；则Luminance_G(Y)=S(I)*LuminanceG(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤20、构造向量Luminance_B，Luminance_B(Y)为向量Luminance_B中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1), RandTest(i,2))的像素点的B通道Luminance值，记为LuminanceB(RandTest(i,1),RandTest(i,2))；则Luminance_B(Y)=S(I)*LuminanceB(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤21、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4；jj=0,1,2,3,4。构造向量F1，F1(ii)表示向量F1的第ii个元素，
，ii=0,1,2,3,4；
则可根据下式
计算的值，R通道的相机响应曲线可用函数表示，
步骤22、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4； jj=0,1,2,3,4。构造向量F2，F2(ii)为向量F2的第ii个元素，
，ii=0,1,2,3,4;
则可根据下式
计算的值，G通道的相机响应曲线可用函数表示，
步骤23、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4；jj=0,1,2,3,4。构造向量F3，F3(ii)为向量F3的第ii个元素，
，ii=0,1,2,3,4;
则可根据下式
计算值，B通道的相机响应曲线可用函数表示，
步骤24、用相机1拍摄任意场景，获得一张照片，用变量Picture_cloud表示；拍摄Picture_cloud时相机1的曝光时间用变量ExposureTime_cloud表示。
步骤25、将Picture_cloud中所有像素点按照从左到右，从上到下的顺序构成矩阵；设矩阵P_C的行数为r，列数为s，P_C(i,j)=(P_CR(i,j), P_CG(i,j), P_CB(i,j))为矩阵P_C第i行第j列元素，P_CR(i,j)、P_CG(i,j)、P_CB(i,j))分别为P_C(i,j)的RGB值，i=1,2,…,r；j=1,2,…,s。
步骤26、令CIER=0.25100031242007470，CIEG=0.6646797942289187，CIEB=0.08431989，构造矩阵L_total，L_total(i，j)为矩阵L_total的第i行第j列元素，则L_total(i,j)可表示为：
i=1,2,…,r；j=1,2,…,s。
步骤27、Picture_cloud对应的高动态范围图像的RGB值分别用矩阵表示，分别为的第i行第j列元素，
。
通过上述步骤可以快速、高精度地生成高动态范围图像。
本发明的积极效果是利用Lumicept软件去模拟一组干净天空的图片，构造一个记录当时天空情况的真实照度模型。对这些一系列曝光时间不同的干净天空的图片，通过上述步骤自动选取少量有限采样点的方法生成相机响应曲线，克服了传统方法中将图片的所有像素点都用于进行相机响应曲线的计算所带来的巨大开销，然而本方法根据模板图片等分灰度区间，并在等分的灰度区间内的随机选点，对每幅图片进行采样点的选取，从而大大减少了计算相机响应曲线时的像素点个数，进而降低了生成相机响应曲线所用的开销。由于所选出的点都有真实照度值与之对应，所以可以得到更加准确的相机响应曲线。最后合成高动态范围图像时，使用了CIE天空标准中的各个通道所占比例，得到的图像的真实照度值更加贴近现实。
具体实施方式
下面结合实施例对本发明做进一步的描述：一种合成高动态范围图像的方法，其特征在于：采用相机1、计算机2，并在计算机2中装有Lumincept9.31.1软件3，按以下的步骤完成：
步骤1、选择合适的场景，场景中只包括一个塔尖形状的物体，塔尖形状的物体远高于场景中的其他物体,场景中还包括没有云彩的天空，这样的场景记为Scene；
步骤2、获取塔尖形状的物体的底部到顶部之间的距离，用变量Bulid_H表示，距离的单位为米，Bulid_H=123.2m；获取塔尖形状的物体在WGS84坐标系下的经度和纬度，分别用变量Location_longitude1、Location_latitude1表示，Location_longitude1=59。55＇48＂，Location_latitude1=30。19＇12＂。
步骤3、调节相机1的位置和姿态，使Scene中塔尖形状的物体的最高处位于相机1的取景器的中心位置。始终保持相机1的位置和姿态不变，获取相机1在WGS84坐标系下的经度和纬度，分别用变量Location_longitude2、Location_latitude2表示，Location_longitude2=59。57＇，Location_latitude2=30。19＇12＂；由小到大的改变相机1的曝光时间A次，A=6，相机1的其他参数不变，通过相机1获取每个曝光时间对应的一张Scene照片，第I张照片的曝光时间设为S(I)，曝光时间的单位为秒，第I张照片用Image(I)表示，I=1,2,...,A，S(1)=1/6400，S(2)=1/2500，S(3)=1/1600，S(4)=1/1000，S(5)=1/800，S(6)=1/640；记录拍摄第一张照片时的当地时间，用变量Time表示，Time为2016年1月1号 12:36。
步骤4、查看谷歌地图，分别找到相机1和Scene中塔尖形状的物体在谷歌地图中的位置，获得他们之间的距离，用变量Distance表示，Distance=208.69，利用Distance和Build_H的值，在Lumicept软件中模拟相机1的位置和方向。
步骤5、Lumicept软件中使用CIE(Commission Internationale de L’Eclairage )标准天空模型，根据Time、Bulid_H、Location_longitude1、Location_latitude1和Location_longitude2、Location_latitude2，在Lumicept软件中复现Image(1)中的场景，将复现的场景记为Simulated_model。
步骤6、设，寻找Image(k)中所有的天空区域，用Sky表示，将Sky中所有像素点按照从左到右，从上到下的顺序构成矩阵Z；将矩阵Z的行数设为m,列数设为n，
Z(c,d)=(R(c,d),G(c,d),B(c,d))为矩阵Z的第c行第d列元素，其中R(c,d)、G(c,d)、B(c,d)分别为Z(c,d)的RGB值，c=1,2,…,m；d=1,2,…,n，本文k=3。
步骤7、构造矩阵Gy，设元素Gy（c,d）是矩阵Gy的第c行第d列元素，Gy(c,d)=B(c,d)，c=1,2,...,m；d=1,2,...,n。
设矩阵Gy中元素的最大值用maxGy表示，矩阵Gy中元素的最小值用minGy表示。
步骤8、根据矩阵Gy中的元素的最大值和最小值来确定采样点的个数N，只需满足方程：
，本文取N=16。
步骤9、构造向量Lrange，Lrange(t)为向量Lrange中的第t个元素，
。
步骤10、构造向量RandNumber，RandNumber(e)为RandNumber的第e个元素，且RandNumber(e)为整数，，e=1,2….,N-1;
步骤11、设在矩阵Gy中与RandNumber(e)相同元素的个数为size(e)，e=1,2….,N；构造矩阵Hp，设元素Hp(i,j)是矩阵Hp的第i行第j列元素，i=1,2,…,size(p)；j=1,2；p=1,2,…,N；构造向量w，w(q)为向量w的第q个元素，w(q)=1，q=1,2,…,size(p)，p=1,2,…,N。
步骤12、按照从上到下，从左到右的顺序，如果Gy(c,d)等于 RandNumber(e)，那么Hp(w(p),1)=c；Hp(w(p),2)=d，w(p)=w(p)+1；构造向量random (i)为向量random的第i个值，random(i)的值是大于1并且小于的任意一个整数，c=1,2,…,m；d=1,2,…,n，i=1,2,…,size(p)，e=p=1,2,…,N。
步骤13、构造矩阵RandTest，其中RandTest(i,j)是矩阵RandTest中的第i行第j列元素，i=1,2,…,N；j=1,2；RandTest(i,1)=Hp(random(i),1)，RandTest(i,2)=Hp(random(i),2)，i=1,2,…,N。
步骤14、将Image(L)上所有像素点按照从左到右，从上到下的顺序构成矩阵ZL，将矩阵ZL的行数设为r,列数设为s,设元素ZL(i,j)=(R_L(i,j),G_L(i,j),B_L(i,j))，ZL(i,j)是矩阵ZL的第i行第j列元素，其中R_L(i,j)、G_L(i,j)、B_L(i,j)分别为ZL(i,j)的RGB值i=1,2,…,r；j=1,2,…,s；L=1,2,…,A。
步骤15、构造一个向量U_R，其中U_R(X)为向量U_R中第X个元素；U_R(X)=R_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤16、构造一个向量U_G，其中U_G(X)为向量U_G中第X个元素;当U_G(X)=G_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤17、构造一个向量U_B，其中U_B (X)为向量U_B中第X个元素;当U_B(X)=B_L(RandTest(i,1),RandTest(i,2))，其中X=(L-1)*N+i；L=1,2,…,A，i=1,2,…,N。
步骤18、构造向量Luminance_R，Luminance_R(Y)为向量Luminance_R中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1),RandTest(i,2))的像素点R通道Luminance值，记为LuminanceR(RandTest(i,1),RandTest(i,2))；则Luminance_R(Y)=S(I)*LuminanceR(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤19、构造向量Luminance_G，Luminance_G(Y)为向量Luminance_G中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1), RandTest(i,2))的像素点的G通道Luminance值，记为LuminanceG(RandTest(i,1),RandTest(i,2))；则Luminance_G(Y)=S(I)*LuminanceG(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤20、构造向量Luminance_B，Luminance_B(Y)为向量Luminance_B中的第Y个元素。使用Lumicept，读取Simulated_model中坐标为(RandTest(i,1), RandTest(i,2))的像素点的B通道Luminance值，记为LuminanceB(RandTest(i,1),RandTest(i,2))；则Luminance_B(Y)=S(I)*LuminanceB(RandTest(i,1),RandTest(i,2))，其中Y=(I-1)*N+i，I=1,2,…,A，i=1,2,…,N。
步骤21、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4；jj=0,1,2,3,4。构造向量F1，F1(ii)表示向量F1的第ii个元素，
，ii=0,1,2,3,4；
则可根据下式
计算的值，R通道的相机响应曲线可用函数表示，
步骤22、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4； jj=0,1,2,3,4。构造向量F2，F2(ii)为向量F2的第ii个元素，
，ii=0,1,2,3,4;
则可根据下式
计算的值，G通道的相机响应曲线可用函数表示，
步骤23、构造一个矩阵表示矩阵的第ii行第jj列元素，
，
ii=0,1,2,3,4；jj=0,1,2,3,4。构造向量F3，F3(ii)为向量F3的第ii个元素，
，ii=0,1,2,3,4;
则可根据下式
计算值，B通道的相机响应曲线可用函数表示，
步骤24、用相机1拍摄任意场景，获得一张照片，用变量Picture_cloud表示；拍摄Picture_cloud时相机1的曝光时间用变量ExposureTime_cloud表示，ExposureTime_cloud=1/2000。
步骤25、将Picture_cloud中所有像素点按照从左到右，从上到下的顺序构成矩阵P_C；设矩阵P_C的行数为r，列数为s，P_C(i,j)=(P_CR(i,j), P_CG(i,j), P_CB(i,j))为矩阵P_C第i行第j列元素，P_CR(i,j)、P_CG(i,j)、P_CB(i,j))分别为P_C(i,j)的RGB值，i=1,2,…,r；j=1,2,…,s。
步骤26、令CIER=0.25100031242007470，CIEG=0.6646797942289187，CIEB=0.08431989，构造矩阵L_total，L_total(i，j)为矩阵L_total的第i行第j列元素，则L_total(i,j)可表示为：
i=1,2,…,r；j=1,2,…,s。
步骤27、Picture_cloud对应的高动态范围图像的RGB值分别用矩阵表示，分别为的第i行第j列元素，
。
通过上述步骤可以快速、高精度地生成高动态范围图像。