本发明属于图像融合
技术领域：
，具体涉及一种基于双色空间的假彩色图像融合方法。
背景技术：
：假彩色图像融合实现不同图像间的信息互补，以增强图像中信息的透明度，对目标清晰、完整、准确的描述，为后续的决策提供可靠的信息支持。基于双色空间的假彩色图像融合技术目前专利未见报道，传统的假彩色图像融合主要运用于红外与可见光、红外与微光的双通道融合中，其使用的主要方法有，基于空域、变换域和传统色彩传递的方法。空域的方法需要对原始图像划分成若干个子图像，选择适当的模糊测度及融合规定来得到融合图像。对于频域的方法，典型的有基于小波变换的融合方法。基于空域的方法尽管有计算简单、复杂度低的优点，但其最大的缺点就是融合后图像容易出现马赛克效应，同时存在一定的模糊失真。基于变换域的方法虽然不存在马赛克效应，但融合图像会出现颜色失真和模糊效应，色彩还原度不高。传统的色彩传递过程主要按照下面公式进行转换：式中：s.t分别表示源图像和参考图像，a，b，c分别表示各颜色空间的分量，和分别表示源图像θ分量和参考图像θ分量的均值，和分别为源图像和参考图像的标准差，根据原始图像的颜色空间分量的期望与标准差值，计算出融合图像的颜色空间分量，最终进行显示，得到假彩色融合图像信息。基于传统的色彩传递融合方法，虽然能够根据不同的参考图像，得到不同的融合效果，但是融合图像颜色与实际场景差距较大，容易造成图像信息被颜色覆盖的情况，造成图像信息丢失。因此如何克服现有技术的不足是目前图像融合
技术领域：
亟需解决的问题。技术实现要素：本发明的目的是为了解决现有技术的不足，提供一种基于双色空间的假彩色图像融合方法，该方法融合图像不存在马赛克现象，色彩还原度高，图像失真小，融合图像信息不会被颜色掩盖，且易于工程实现，能够实时通过fpga进行算法移植。为实现上述目的，本发明采用的技术方案如下：一种基于双色空间的假彩色图像融合方法，其特征在于，包括如下步骤：步骤(1)，选择一幅彩色图像作为参考图像，之后根据参考图像的rgb空间信息t(r0)、t(g0)、t(b0)值，分别计算出参考图像对应的ycbcr空间信息p(y0)、p(cb0)、p(cr0)值；p(y0)＝0.257*t(r0)+0.564*t(g0)+0.098*t(b0)+16；p(cb0)＝-0.148*t(r0)-0.291*t(g0)+0.439*t(b0)+128；p(cr0)＝0.439*t(r0)-0.368*t(g0)-0.071*t(b0)+128；步骤(2)，根据参考图像t(r0)、t(g0)和t(b0)值，分别计算出参考图像r空间信息、g空间信息和b空间信息的均值，分别记为u(r0)、u(g0)、u(b0)；再分别计算出参考图像r空间信息、g空间信息和b空间信息的标准差值，分别记为σ(r0)、σ(g0)、σ(b0)，计算公式如下：同时，根据参考图像p(y0)、p(cb0)和p(cr0)值，分别计算出参考图像y空间信息、cb空间信息和cr空间信息的均值，分别记为u(y0)、u(cb0)、u(cr0)；再分别计算出参考图像y空间信息、cb空间信息和cr空间信息的标准差值，分别记为σ(y0)、σ(cb0)、σ(cr0)；步骤(3)，根据待处理视频图像的灰度值y，计算新的r空间信息、g空间信息、b空间信息、y空间信息、cb空间信息和cr空间信息，分别记为t(rs)、t(gs)、t(bs)、p(ys)、p(cbs)、p(crs)；所述的待处理视频图像有两幅，为同一画面的两个通道的视频图像，记为图像ir和图像llp；计算公式如下：t(rs)＝(y(ir)+y(llp))/2；t(gs)＝max(y(ir)，y(llp))；t(bs)＝y(llp)；p(ys)＝(y(ir)+y(llp))/2；p(cbs)＝max(y(ir)，y(llp))；p(crs)＝y(llp)；步骤(4)，根据步骤(3)计算出来的r空间信息、g空间信息和b空间信息t(rs)、t(gs)和t(bs)，分别计算出对应均值和标准差值，均值分别记为u(rs)、u(gs)、u(bs)，标准差值分别记为σ(rs)、σ(gs)、σ(bs)；同时，根据步骤(3)计算出来的y空间信息、cb空间信息和cr空间信息p(ys)、p(cbs)和p(crs)，分别计算出对应均值和标准差值，均值分别记为u(ys)、u(cbs)、u(crs)，标准差值分别记为σ(ys)、σ(cbs)、σ(crs)；步骤(5)，根据步骤(2)计算所得值u(y0)、u(cb0)、u(cr0)、σ(y0)、σ(cb0)、σ(cr0)和步骤(4)计算值u(ys)、u(cbs)、u(crs)、σ(ys)、σ(cbs)、σ(crs)，分别计算出第二次ycbcr值，记为p(yt)、p(cbt)、p(crt)；计算公式如下：步骤(6)，将步骤(5)计算得到的p(yt)、p(cbt)、p(crt)转成成rgb空间信息，记为t(rt)、t(gt)、t(bt)；转换公式如下：t(rt)＝1.164*(p(yt)-16)+1.596*(p(crt)-128)；t(gt)＝1.164*(p(yt)-16)-0.392*(p(cbt)-128)-0.813*(p(crt)-128)；t(bt)＝1.164*(p(yt)-16)+2.017*(p(cbt)-128)；步骤(7)，根据步骤(2)和步骤(4)计算出来的u(r0)、u(g0)、u(b0)、σ(r0)、σ(g0)、σ(b0)和u(rs)、u(gs)、u(bs)、σ(rs)、σ(gs)、σ(bs)，计算第二次rgb值，记为t(rt1)、t(gt1)、t(bt1)；计算公式如下：步骤(8)，根据步骤(6)计算得到的t(rt)、t(gt)、t(bt)和步骤(7)计算得到的t(rt1)、t(gt1)、t(bt1)，计算第三次rgb值，记为t(rt2)、t(gt2)、t(bt2)；计算公式如下：t(rt2)＝α*t(rt)+(1-α)*t(rt1)t(gt2)＝α*t(gt)+(1-α)*t(gt1)；t(bt2)＝α*t(bt)+(1-α)*t(bt1)其中，0＜α＜1步骤(9)，将步骤(8)得到的t(rt2)、t(gt2)、t(bt2)转为ycbcr空间下，得到图像ir和图像llp的融合假彩色图像；p(yl)＝0.257*t(rt2)+0.564*t(gt2)+0.098*t(bt2)+16；p(cbl)＝-0.148*t(rt2)-0.291*t(gt)+0.439*t(bt2)+128；p(crl)＝0.439*t(rt2)-0.368*t(gt2)-0.071*t(bt2)+128。进一步，优选的是，能够实时处理双通道输入50hz/s的视频图像，参考图像与将要融合输出的彩色图像色彩还原度达到80％。进一步，优选的是，所述的α为0.5。本发明中，参考图像的选择可以根据不同的场景进行选取。本发明能够实时通过fpga实现，融合效果更佳有利于人眼观察，不会产生马赛克现象，融合图像信息损失小，能够在保证色彩度的同事保留图像原始信息。本发明中进行了多次ycbcr空间和rgb空间的转换，主要目的在于对参考图像和待融合图像进行均值和标准差的不同空间统计，有利于在rgb空间实现融合。本发明转换次数为实现中提出的最优次数，通过这些转换得到的融合输出图像色彩还原度更高，如果转换次数再增加，将增加缓存，造成延时，图像输出帧频会小于25hz/s。本发明与现有技术相比，其有益效果为：1、采用本发明方法，视频图像ir和视频图像llp融合结果不会产生明显的滞后情况，实时进行融合；2、通过本发明方法，融合假彩色更加接近真实色彩，根据参考图像的选择，能够实现接近与森林、沙漠、雪地、海洋等色彩模式；3、本发明方法实现的假彩色模式不会出现马赛克现象，图像失真度小，能够最大化保留原始图像信息。4、与现有技术相比，本发明方法能够实时处理双通道输入50hz/s的视频图像，参考图像与将要融合输出的彩色图像色彩还原度可以达到80％。附图说明图1为参考图像c；图2为第一通道的待处理视频图像ir；图3为第二通道的待处理视频图像llp；图4为融合得到的假彩色图像；图5为本发明方法流程图。具体实施方式下面结合实施例对本发明作进一步的详细描述。本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限定本发明的范围。实施例中未注明具体技术或条件者，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用材料或设备未注明生产厂商者，均为可以通过购买获得的常规产品。本发明实施例中均值和标准差的计算通式如下：为了说明具体的实施过程，假设参考图像c的rgb信息是3×3的矩阵大小，其中t(r0)，t(g0)，t(b0)分别表示参考图像的r空间信息、g空间信息、b空间信息，如表1-3所示。表1参考图像c的r空间信息t11(r0)t12(r0)t13(r0)t21(r0)t22(r0)t23(r0)t31(r0)t32(r0)t33(r0)表2参考图像c的g空间信息t11(g0)t12(g0)t13(g0)t21(g0)t22(g0)t23(g0)t31(g0)t32(g0)t33(g0)表3参考图像c的b空间信息t11(b0)t12(b0)t13(b0)t21(b0)t22(b0)t23(b0)t31(b0)t32(b0)t33(b0)为了说明具体的实施过程，假设待处理视频图像有两幅，为同一画面的两个通道的视频图像，记为图像ir和图像llp，均为3×3的矩阵大小，如表4、表5所示，其中y(ir)和y(llp)分别代表输入源图像ir和源图像llp信息的灰度值；表4ir图3×3矩阵y11(ir)y12(ir)y13(ir)y21(ir)y22(ir)y23(ir)y31(ir)y32(ir)y33(ir)表5llp图3×3矩阵y11(llp)y12(llp)y13(llp)y21(llp)y22(llp)y23(llp)y31(llp)y32(llp)y33(llp)下面针对参考图像和源图像矩阵模式，将给出详细的实施步骤，流程图如图5所示：步骤(1)，根据参考图像c(如图1)的rgb矩阵(如表1-3所示)分别计算出参考图像对应的ycbcr空间信息y0，cb0，cr0值，分别记为p(y0)、p(cb0)、p(cr0)详细的计算公式如下：p11(y0)＝0.257*t11(r0)+0.564*t11(g0)+0.098*t11(b0)+16；p12(y0)＝0.257*t12(r0)+0.564*t12(g0)+0.098*t12(b0)+16；p13(y0)＝0.257*t13(r0)+0.564*t13(g0)+0.098*t13(b0)+16；……p33(y0)＝0.257*t33(r0)+0.564*t33(g0)+0.098*t33(b0)+16。p11(cb0)＝-0.148*t11(r0)-0.291*t11(g0)+0.439*t11(b0)+128；p12(cb0)＝-0.148*t12(r0)-0.291*t12(g0)+0.439*t12(b0)+128；p13(cb0)＝-0.148*t13(r0)-0.291*t13(g0)+0.439*t13(b0)+128；……p33(cb0)＝-0.148*t33(r0)-0.291*t33(g0)+0.439*t33(b0)+128。p11(cr0)＝0.439*t11(r0)-0.368*t11(g0)-0.071*t11(b0)+128；p12(cr0)＝0.439*t12(r0)-0.368*t12(g0)-0.071*t12(b0)+128；p13(cr0)＝0.439*t13(r0)-0.368*t13(g0)-0.071*t13(b0)+128；……p33(cr0)＝0.439*t33(r0)-0.368*t33(g0)-0.071*t33(b0)+128。那么，参考图像对应的ycbcr空间信息p(y0)、p(cb0)、p(cr0)矩阵值如表6-8所示。表6ycbcr空间信息y0矩阵p11(y0)p12(y0)p13(y0)p21(y0)p22(y0)p23(y0)p31(y0)p32(y0)p33(y0)表7ycbcr空间信息cb0矩阵p11(cb0)p12(cb0)p13(cb0)p21(cb0)p22(cb0)p23(cb0)p31(cb0)p32(cb0)p33(cb0)表8ycbcr空间信息cr0矩阵步骤(2)，根据参考图像t(r0)、t(g0)和t(b0)值，分别计算出参考图像r空间信息、g空间信息和b空间信息的均值，分别记为u(r0)、u(g0)、u(b0)；再分别计算出参考图像r空间信息、g空间信息和b空间信息的标准差值，分别记为σ(r0)、σ(g0)、σ(b0)，计算公式如下：同时，根据参考图像p(y0)、p(cb0)和p(cr0)值，分别计算出参考图像y空间信息、cb空间信息和cr空间信息的均值，分别记为u(y0)、u(cb0)、u(cr0)；再分别计算出参考图像y空间信息、cb空间信息和cr空间信息的标准差值，分别记为σ(y0)、σ(cb0)、σ(cr0)，计算公式如下：步骤(3)，根据待处理视频图像的灰度值y，计算新的r空间信息、g空间信息、b空间信息、y空间信息、cb空间信息和cr空间信息，分别记为t(rs)、t(gs)、t(bs)、p(ys)、p(cbs)、p(crs)；所述的待处理视频图像有两幅，为同一画面的两个通道的视频图像，记为图像ir和图像llp；计算公式如下所示；……t11(gs)＝max(y11(ir)，y11(llp))；t12(gs)＝max(y12(ir)，y12(llp))；t13(gs)＝max(y13(ir)，y13(llp))；……t33(gs)＝max(y33(ir)，y33(llp))。p11(bs)＝y11(llp)；p12(bs)＝y12(llp)；p13(bs)＝y13(llp)；……p33(bs)＝y33(llp)。那么，t(rs)、t(gs)、t(bs)矩阵值如表7-9所示。表7rs矩阵值表8gs矩阵值t11(gs)t12(gs)t13(gs)t21(gs)t22(gs)t23(gs)t31(gs)t32(gs)t33(gs)表9bs矩阵值t11(bs)t12(bs)t13(bs)t21(bs)t22(bs)t23(bs)t31(bs)t32(bs)t33(bs)……p11(cbs)＝max(y11(ir)，y11(llp))；p12(cbs)＝max(y12(ir)，y12(llp))；p13(cbs)＝max(y13(ir)，y13(llp))；……p33(cbs)＝max(y33(ir)，y33(llp))。p11(crs)＝y11(llp)；p12(crs)＝y12(llp)；p13(crs)＝y13(llp)；……p33(crs)＝y33(llp)。那么，p(ys)、p(cbs)、p(crs)矩阵值如表12-14所示。表12ys矩阵值p11(ys)p12(ys)p13(ys)p21(ys)p22(ys)p23(ys)p31(ys)p32(ys)p33(ys)表13cbs矩阵值p11(cbs)p12(cbs)p13(cbs)p21(cbs)p22(cbs)p23(cbs)p31(cbs)p32(cbs)p33(cbs)表14crs矩阵值p11(crs)p12(crs)p13(crs)p21(crs)p22(crs)p23(crs)p31(crs)p32(crs)p33(crs)步骤(4)，根据步骤(3)计算出来的r空间信息、g空间信息和b空间信息t(rs)、t(gs)和t(bs)，分别计算出对应均值和标准差值，均值分别记为u(rs)、u(gs)、u(bs)，标准差值分别记为σ(rs)，σ(gs)，σ(bs)；计算公式如下：同时，根据步骤(3)计算出来的y空间信息、cb空间信息和cr空间信息p(ys)、p(cbs)和p(crs)，分别计算出对应均值和标准差值，均值分别记为u(ys)、u(cbs)、u(crs)，标准差值分别记为σ(ys)、σ(cbs)、σ(crs)；计算公式如下：步骤(5)，根据步骤(2)计算所得值u(y0)、u(cb0)、u(cr0)、σ(y0)、σ(cb0)、σ(cr0)和步骤(4)计算值u(ys)、u(cbs)、u(crs)、σ(ys)、σ(cbs)、σ(crs)，分别计算第二次ycbcr值，记为p(yt)、p(cbt)、p(crt)；计算公式如下：………………那么，p(yt)、p(cbt)、p(crt)矩阵值如表15-17所示。表15yt矩阵值p11(yt)p12(yt)p13(yt)p21(yt)p22(yt)p23(yt)p31(yt)p32(yt)p33(yt)表16cbt矩阵值p11(cbt)p12(cbt)p13(cbt)p21(cbt)p22(cbt)p23(cbt)p31(cbt)p32(cbt)p33(cbt)表17crt矩阵值步骤(6)，将步骤(5)计算得到的p(yt)、p(cbt)、p(crt)转成成rgb空间信息，记为t(rt)、t(gt)、t(bt)；转换公式如下：t11(rt)＝1.164*(p11(yt)-16)+1.596*(p11(crt)-128)；t12(rt)＝1.164*(p12(yt)-16)+1.596*(p12(crt)-128)；t13(rt)＝1.164*(p13(yt)-16)+1.596*(p13(crt)-128)；……t33(rt)＝1.164*(p33(yt)-16)+1.596*(p33(crt)-128)。t11(gt)＝1.164*(p11(yt)-16)-0.392*(p11(cbt)-128)-0.813*(p11(crt)-128)；t12(gt)＝1.164*(p12(yt)-16)-0.392*(p12(cbt)-128)-0.813*(p12(crt)-128)；t13(gt)＝1.164*(p13(yt)-16)-0.392*(p13(cbt)-128)-0.813*(p13(crt)-128)；……t33(gt)＝1.164*(p33(yt)-16)-0.392*(p33(cbt)-128)-0.813*(p33(crt)-128)。t11(bt)＝1.164*(p11(yt)-16)+2.017*(p11(cbt)-128)；t12(bt)＝1.164*(p12(yt)-16)+2.017*(p12(cbt)-128)；t13(bt)＝1.164*(p13(yt)-16)+2.017*(p13(cbt)-128)；……t33(bt)＝1.164*(p33(yt)-16)+2.017*(p33(cbt)-128)。那么，t(rt)、t(gt)、t(bt)矩阵值如表18-20所示。表18rt矩阵值t11(rt)t12(rt)t13(rt)t21(rt)t22(rt)t23(rt)t31(rt)t32(rt)t33(rt)表19gt矩阵值t11(gt)t12(gt)t13(gt)t21(gt)t22(gt)t23(gt)t31(gt)t32(gt)t33(gt)表20bt矩阵值t11(bt)t12(bt)t13(bt)t21(bt)t22(bt)t23(bt)t31(bt)t32(bt)t33(bt)步骤(7)，根据步骤(2)和步骤(4)计算出来的u(r0)，u(g0)，u(b0)，σ(r0)，σ(g0)，σ(b0)和u(rs)，u(gs)，u(bs)，σ(rs)，σ(gs)，σ(bs)，计算第二次rgb值，记为t(rt1)、t(gt1)、t(bt1)；计算公式如下：………………那么，t(rt1)、t(gt1)、t(bt1)矩阵值如表21-23所示。表21rt1矩阵值t11(rt1)t12(rt1)t13(rt1)t21(rt1)t22(rt1)t23(rt1)t31(rt1)t32(rt1)t33(rt1)表22gt1矩阵值t11(gt1)t12(gt1)t13(gt1)t21(gt1)t22(gt1)t23(gt1)t31(gt1)t32(gt1)t33(gt1)表23bt1矩阵值t11(bt1)t12(bt1)t13(bt1)t21(bt1)t22(bt1)t23(bt1)t31(bt1)t32(bt1)t33(bt1)步骤(8)，根据步骤(6)计算得到的t(rt)、t(gt)、t(bt)和步骤(7)计算得到的t(rt1)、t(gt1)、t(bt1)，计算第二次rgb值，记为t(rt2)、t(gt2)、t(bt2)，其中α＝0.5；计算公式如下：t11(rt2)＝0.5*t11(rt)+0.5*t11(rt1)；t12(rt2)＝0.5*t12(rt)+0.5*t12(rt1)；t13(rt2)＝0.5*t13(rt)+0.5*t13(rt1)；……t33(rt2)＝0.5*t33(rt)+0.5*t133(rt1)。t11(gt2)＝0.5*t11(gt)+0.5*t11(gt1)；t12(gt2)＝0.5*t12(gt)+0.5*t12(gt1)；t13(gt2)＝0.5*t13(gt)+0.5*t13(gt1)；……t33(gt2)＝0.5*t33(gt)+0.5*t133(gt1)。t11(bt2)＝0.5*t11(bt)+0.5*t11(bt1)；t12(bt2)＝0.5*t12(bt)+0.5*t12(bt1)；t13(bt2)＝0.5*t13(bt)+0.5*t13(bt1)；……t33(bt2)＝0.5*t33(bt)+0.5*t133(bt1)。表21rt2矩阵值t11(rt2)t12(rt2)t13(rt2)t21(rt2)t22(rt2)t23(rt2)t31(rt2)t32(rt2)t33(rt2)表22gt2矩阵值t11(gt2)t12(gt2)t13(gt2)t21(gt2)t22(gt2)t23(gt2)t31(gt2)t32(gt2)t33(gt2)表23bt2矩阵值t11(bt2)t12(bt2)t13(bt2)t21(bt2)t22(bt2)t23(bt2)t31(bt2)t32(bt2)t33(bt2)步骤(9)，将步骤(8)得到的t(rt2)、t(gt2)、t(bt2)转为ycbcr空间下进行fpga假彩色显示，ycbcr空间信息记为p(yl)、p(cbl)、p(crl)，得到图像ir和图像llp的融合假彩色图像。p11(yl)＝0.257*t11(t(rt2))+0.564*t11(t(gt2))+0.098*t11(t(bt2))+16；p12(yl)＝0.257*t12(t(rt2))+0.564*t12(t(gt2))+0.098*t12(t(bt2))+16；p13(yl)＝0.257*t13(t(rt2))+0.564*t13(t(gt2))+0.098*t13(t(bt2))+16；……p33(yl)＝0.257*t33(t(rt2))+0.564*t33(t(gt2))+0.098*t33(t(bt2))+16。p11(cbl)＝-0.148*t11(rt2)-0.291*t11(gt2)+0.439*t11(bt2)+128；p12(cbl)＝-0.148*t12(rt2)-0.291*t12(gt2)+0.439*t12(bt2)+128；p13(cbl)＝-0.148*t13(rt2)-0.291*t13(gt2)+0.439*t13(bt2)+128；……p33(cbl)＝-0.148*t33(rt2)-0.291*t33(gt2)+0.439*t33(bt2)+128。p11(crl)＝0.439*t11(rt2)-0.368*t11(gt2)-0.071*t11(bt2)+128；p12(crl)＝0.439*t12(rt2)-0.368*t12(gt2)-0.071*t12(bt2)+128；p13(crl)＝0.439*t13(rt2)-0.368*t13(gt2)-0.071*t13(bt2)+128；……p33(crl)＝0.439*t33(rt2)-0.368*t33(gt2)-0.071*t33(bt2)+128。那么最终fpga进行显示的图像矩阵则为表24-26所示，得到的融合图像如图4所示。表24yl矩阵值p11(yl)p12(yl)p13(yl)p21(yl)p22(yl)p23(yl)p31(yl)p32(yl)p33(yl)表25cbl矩阵值p11(cbl)p12(cbl)p13(cbl)p21(cbl)p22(cbl)p23(cbl)p31(cbl)p32(cbl)p33(cbl)表26crl矩阵值p11(crl)p12(crl)p13(crl)p21(crl)p22(crl)p23(crl)p31(crl)p32(crl)p33(crl)以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。当前第1页12