一种视频通信伽玛特性的校正方法和装置的制作方法

文档序号:7963569阅读:104来源:国知局
专利名称:一种视频通信伽玛特性的校正方法和装置的制作方法
技术领域
本发明涉及视频通讯技术领域,具体涉及一种视频通信伽玛特性的校正方法和装置。
背景技术
视频通信,尤其是多方视频通信,目前正在随着宽带网络的迅速发展,得到日益广泛的应用。在国内和国际上,视频会议和可视电话业务正在成为NGN(Next Generation Network,下一代网络)上的基本业务。各国的电信运营商也非常重视这个市场机会。预期在未来几年中,视频通信业务将成为运营商重要的业务增长点。
发展视频通讯业务的一个关键问题是提高端到端(End-to-end)的用户体验(User Experience或者Quality of Experience)。用户体验中除了网络的QoS如丢包、延迟、抖动、R因子等参数外,对于视频,由于各个环节引起的Gamma(伽玛)非线性问题造成的对亮度信号的畸变(Distortion),也是影响最终用户体验的重要因素。
目前,对于提高端到端用户体验的方法和技术主要集中在保证网络QoS和视频压缩编码相关的前后处理(Pre-processing,post-processing)方面。对于Gamma特性引起的亮度畸变问题,缺乏关注和系统的解决方法。
下面对Gamma特性进行简要介绍。
视频通信的过程为需要被传送的场景如人物、背景、文件等的光信号进入到视频通信终端(下文简称终端)如摄像机/摄像头等,经过A/D转换成数字图像信号,再经过压缩编码,传送出去,到达对方终端,然后,经过去压缩(decompression)解码还原为数字图像信号,再在显示设备上显示出来,最终又变成光信号被人眼感知。在上述过程中,图像亮度信号经过了多个环节。根据定义,Gamma特性就是一个环节的图像亮度信号的输入-输出关系不是线性的,而是一种非线性的关系。这里的图像亮度信号(Luminance)是一种广义的亮度信号,即一开始的光信号,到电信号,再到数字化的图像亮度/灰度信号,每个阶段的信号都含有图像亮度信号的信息,因此,广义地说,图像亮度信号经过了多个环节。
单个环节Gamma特性的一般模型如附图1所示。
图1中,输入亮度信号和输出亮度信号的非线性的关系可以表示为Lout=G(Lin),其中,Lout为输出亮度信号,Lin为输入亮度信号,函数G(.)为一个非线性函数。
典型的Gamma特性示例如附图2所示。
图2中的每一个方块中标注的数字为亮度值,方块的灰度表示亮度信号的明亮程度。图2(a)中,上面的一行灰度方块的亮度是线性递增的,即从0.1递增到1.0,下面一行灰度方块的亮度是按照幂函数规律递增的,也就是说,下面一行灰度方块的亮度经过了Gamma非线性的失真影响。图2(b)中给出的是以曲线表示的Gamma特性。
当多个环节级联(cascading)起来或者说串联起来时,则总的Gamma特性等于各个环节Gamma函数的复合(composition)。
多个环节级联的Gamma特性的一般模型如附图3所示。
图3中,每个环节的输入亮度信号和输出亮度信号的非线性的关系分别为Lout=G(1)(Lin)、Lout=G(2)(Lin)、Lout=G(3)(Lin)。
由此可以得知各个环节Gamma函数的复合为公式(1)所示GCT(.)=G(1)(.)οG(2)(.)οG(3)(.)......G(n-1)(.)οG(n)(.)(1)lout=GCT(lin)=G(n)(G(n-1)(G(n-2)(.......G(2)(G(1)(lin)))))其中,“ο”表示函数的复合运算。CT表示cascaded total,即级联总Gamma的意思。
在实际中,Gamma非线性是由不同原因引起的。显示设备如CRT显示器的Gamma特性在理想状况下是Lout=Lin2.2(2)而对应的摄像机/摄像头的理想的Gamma特性是Lout=Lin0.45(3)从Gamma问题的起源来看,Gamma问题起源于CRT显示器,因为其Gamma值是2.2,为了补偿掉这个非线性,在摄像机中引入了Gamma值0.45。在这个例子中,Gamma特性的形式是一个幂函数(Power Function)。需要说明的是,这里的输入和输出亮度信号都是在各自的坐标空间中进行了规一化(Normalized)的,即0≤Lout≤1,0≤Lin≤1。而其它类型的显示器比如液晶等,其Gamma函数的形式或者会有所不同,或者虽然在形式上也是幂函数,但是参数不同。
理想的情况是输入亮度信号和输出亮度信号之间存在线性关系,即Lout=Lin。要获得线性关系,必须对于具有非线性Gamma特性的环节进行Gamma校正(Gamma Correction)。
对一个Gamma环节的Gamma校正原理图如附图4所示。
图4中,对于一个环节来说,其Gamma特性给定即Lout=Gg(Lin),这样,可以用另外一个校正环节Lout=Gc(Lin)和它进行级联,来使得总的Gamma特性成为真正的线性关系,从而达到校正掉给定环节的非线性的目的。
显然,Gg(.)和Gc(.)互为反函数。在一般情况下,对于一个函数,要获得其反函数不一定有解,而且,即使存在反函数,也无法用计算的方法获得。
实际应用中,更多的情况是存在多个Gamma环节的情况,对多个Gamma环节的Gamma校正原理图如附图5所示。
图5中,校正环节需要插入到前后两个给定Gamma环节之间。前给定环节的Gamma特性即Lout=Ga(Lin),后给定环节的Gamma特性即Lout=Gp(Lin)。此时,校正环节中的Gc(.)非常复杂,Gc(.)和Ga(.)或者Gp(.)之间不再是简单的反函数关系了。
实现上述Gamma校正方法,其前提是能够对于一个给定的Gamma环节或者多个给定的Gamma环节的级联确定Gamma特性参数。这里的Gamma特性参数就是Gamma特性函数曲线的参数。
在通信的一般情况下,校正需要涉及到两个以上的通信终端。比如在一个两方视频通信中,终端A的视频传送到终端B,那么这路视频的校正就同时涉及到终端A上的Gamma环节和终端B上的Gamma环节。
目前,确定Gamma特性参数的方法主要有两种方法一仪器测量方法。即通过专用仪器测量出Gamma特性函数曲线上的一些点,然后,采用数据拟合的方法来进行曲线拟合,以确定Gamma特性参数。
方法二采用输入亮度信号和输出亮度信号的方法。即对于单个给定的Gamma环节,只要Gg(.)满足一定条件,就可以找到对Gg(.)进行Gamma校正的Gc(.);对于多个给定的Gamma环节,只要Ga(.)、Gp(.)满足一定条件,就可以找到对Ga(.)和Gp(.)进行Gamma校正的Gc(.)。
从上述两种方法的描述可以看出,目前确定Gamma特性参数的实现方法都有一个前提条件,即明确知道需要被确定Gamma特性参数的Gamma环节的输入亮度信号和输出亮度信号的具体数值,也就是明确知道Gamma环节的输入亮度信号和输出亮度信号的全部知识,因此,上述两种确定Gamma特性参数的实现方法均属于非盲测量方法。
非盲测量方法的实现原理如附图6所示。
图6中,Gamma特性参数测量系统利用Gamma环节的输入亮度信号全部知识、输出亮度信号全部知识测量出Gamma特性参数,这里的Gamma环节可以是单个Gamma环节,也可以是级联的Gamma环节。
但是,非盲测量方法在实际应用中的适用范围是非常有限的,也就是说,上述前提条件在很多情况下是不能满足的,下面例举目前常见的三种不能满足上述前提条件的应用。
应用情景一对于IPTV(Internet Protocol Television)等流媒体业务和应用,由于节目制作过程中,已经受到了视频输入设备的Gamma特性的影响,在节目播出的时候,尤其是点播等情况,已经无法获得原来节目制作时候用于采集视频信号的视频输入设备的Gamma特性,而且,也不可能对视频输入设备的Gamma特性参数进行测量了。
应用情景二对于数据会议等应用也存在上述问题。目前,视频会议的发展和数据会议的发展同步,两者完善的结合,对于协作应用(collaborativeapplications)有很大的意义。在企业等环境中,上述协作应用业务有强烈的市场需求。但是,在数据会议应用中,很多多媒体资料比如图片等的来源是不可考的,很难获得当时生成这些数据的视频输入设备的Gamma特性,而且,也不可能对视频输入设备的Gamma特性参数进行测量了。
应用情景三对于面向千万家庭用户的公众视频通信业务来说,为了降低成本和视频通信业务使用门槛,往往大量的采用廉价摄像头,尤其是那些非常便宜的USB接口摄像头。这些廉价的视频输入设备的Gamma特性曲线和标准的Lout=Lin0.45相差很远,甚至根本不是幂函数的形式。而且从这些廉价摄像头的出厂技术资料中一般无法获取其Gamma特性参数。甚至有些廉价的摄像头根本就没有出厂技术资料。用户在家里使用这些摄像头时,也不可能通过上述确定Gamma特性参数的方法来获得Gamma特性参数。
以上三种应用是非常重要的,而且,上述三种应用都有很大的市场潜力,尤其是IPTV和协作数据会议市场的发展非常快。视频通信要真正用于巨大的市场,必须依靠走公众运营的道路来吸引千家万户,这样就要求入门条件一定要非常低,视频输入设备的价格要非常低廉。目前确定Gamma特性参数的非盲测量方法使Gamma校正难于应用,而且,通过专用仪器测量的方法提高了视频通信业务的实现成本。

发明内容
本发明的目的在于,提供一种伽玛特性的校正方法和装置,以实现降低视频通信成本,提高Gamma校正易用性的目的。
为达到上述目的,本发明提供的一种伽玛特性的校正方法,包括a、获取输入、输出图像的亮度直方图;b、建立输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系;c、根据所述输入、输出图像的亮度直方图对所述数学关系进行求解,以确定伽玛特性参数;d、根据所述伽玛特性参数对伽玛环节进行伽玛校正。
以下技术方案为可选的技术方案。
所述伽玛环节为单个给定伽玛环节、或者多个给定伽玛环节的级联组合。
所述步骤a中获取输入图像的亮度直方图的步骤包括根据已知统计特性的测试图像获取输入图像的亮度直方图;或者根据图像统计特性相似性从预先存储的直方图信息中选取输入图像的亮度直方图。
所述步骤b包括输入、输出亮度信号各自的亮度分布概率密度函数fe(x,t)、fr(x,t)、以及Gamma特性函数之间的数学关系为d(e;p)fr(r)=fe(e);其中r=g(e,p),e∈
;亮度分布概率密度函数和亮度直方图之间的数学关系为1Nfs(2k+12N,t)=Prob{kN≤s(t)≤k+1N},]]>其中k=0,1,2,....,N-1;将输入、输出信号各自的亮度分布概率密度函数和Gamma特性函数之间的数学关系转化为输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))=he(k);]]>
其中k=0,1,2,...,N-1,N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,he(k)为输入亮度信号的亮度直方图,he(k)为输出亮度信号的亮度直方图,函数[.]表示取整函数。
所述步骤c包括利用直接求解方程法、或者非线性函数优化方法进行求解,以确定伽玛特性参数。
非线性函数优化方法包括根据输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系构造代价函数J(p)=Σk=0N-1(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))-he(k))2(d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+;]]>其中N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,M为伽玛特性参数的个数;在M个参数中选取M-1个参数,将参数空间的维度降低为M-1维;确定M-1个参数的方法为确定参数向量ptrue,使对于任意p∈RM-1,关系J(p)>=J(pture)成立,即寻找到代价函数的全局最小点,参数向量ptrue即为M-1个参数;利用关系g(1;p)=1结合ptrue,确定第M个伽玛特性参数。
确定ptrue的方法包括组合优化方法、神经网络方法、蛮力搜索方法。
所述蛮力搜索方法包括步骤将维度降低一维的伽玛特性参数空间划分为多个超立方体;选取初始搜索点,并根据预定顺序从该初始搜索点所在的超立方体开始进行遍历搜索;计算搜索过程中进入的每个超立方体的几何中心坐标Q,并根据所述代价函数的表达式以及搜索过程中进入的超立方体计算J(Q);如果J(Q)小于等于预定门限,或者搜索过的超立方体满足预定条件,则将本次搜索进入的超立方体的几何中心坐标Q作为ptrue,搜索过程结束;否则,继续搜索过程。
所述初始搜索点根据实际应用中伽玛特性参数的经验数值来设置。
所述遍历搜索包括以初始搜索点所在的超立方体作为初始超立方体,将包围初始超立方体的超立方体按照距离初始超立方体的远近划分为依次包裹前一层的多层超立方体阵列,并逐层进行搜索。
在蛮力搜索方法中,维度降低一维的伽玛特性参数空间被划分为多个粗粒度的超立方体,所述遍历搜索包括以初始搜索点所在的超立方体作为初始超立方体,将包围初始超立方体的超立方体按照距离初始超立方体的远近划分为依次包裹前一层的多层超立方体阵列,并逐层进行搜索;将搜索到的满足条件的超立方体作为新的伽玛特性参数空间,并将其划分为更细粒度的超立方体,依次类推,进行从粗到细的逐层搜索。
本发明还提供一种伽玛特性的校正装置,所述装置包括获取直方图模块、存储模块、计算模块和伽玛校正模块;获取直方图模块用于获取输入、输出图像的亮度直方图;存储模块用于存储输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系;计算模块用于根据输入、输出图像的亮度直方图对存储模块中的所述数学关系进行求解,以确定伽玛特性参数;伽玛校正模块用于根据所述伽玛特性参数对伽玛环节进行伽玛校正。
所述装置位于视频数据源设备中、和/或位于视频通信网络的中间设备中、和/或位于视频数据目的设备。
通过上述技术方案的描述可知,本发明提供的Gamma校正方法仅需要输入亮度信号的直方图和输出亮度信号的直方图即可,这样,本发明的Gamma参数确定方法可以称为半盲测量方法;由于本发明不需要输入亮度信号的全部知识,而且,本发明可以通过已知统计特性的测试图像、或者预先存储的直方图信息中确定输入图像的亮度直方图,因此,本发明的Gamma校正方法具有很高的应用可行性;本发明的伽玛校正方法特别适用于IPTV、协作数据会议、广泛使用低端视频输入设备的公众视频通信;本发明在采用蛮力搜索方法来确定伽玛特性函数的参数时,采用了逐层搜索、从粗到细搜索、根据实际情况选取初始搜索点等方法,提高了搜索效率;从而通过本发明提供的技术方案实现了提高Gamma校正易用性,拓宽Gamma校正的应用范围,提高用户体验和服务质量的目的。


图1是环节Gamma特性的模型示意图;图2(a)是Gamma特性示意图一;图2(b)是Gamma特性示意图二;图3是多个环节级联的Gamma特性的模型示意图;图4是对一个Gamma环节的Gamma校正原理示意图;图5是对多个Gamma环节的Gamma校正原理示意图;图6是现有技术中的非盲测量方法的实现原理示意图;图7是本发明实施例的半盲Gamma特性参数确定方法的实现原理示意图;图8是视频信号的亮度直方图示例图;图9是本发明实施例的Gamma校正示意图;图10是本发明实施例的蛮力搜索方法的原理示意图;图11是本发明实施例的初始超立方体和其外围多层超立方体在二维情况下的示意图;图12是本发明实施例的蛮力搜索方法示意图。
具体实施例方式
在很多应用场景中,输出亮度信号的具体数值即全部知识是可知的,而输入亮度信号的具体数值即全部知识是不可知的。虽然输入亮度信号的全部知识不可知,但是,输入亮度信号的部分知识即输入亮度信号的先验知识(a prioriknowledge)是可以获得的,如输入亮度信号的概率分布统计规律等,本发明正是利用了输入亮度信号的先验知识和输出亮度信号的知识来确定Gamma特性参数,以进行Gamma校正的。由于在本发明的技术方案中,利用了输入亮度信号的部分知识,所以,本发明中确定Gamma特性参数的方法可以称为半盲Gamma特性参数确定方法。
本发明的半盲Gamma特性参数确定方法的实现原理如附图7所示。
图7中,对于需要确定Gamma特性的环节,首先需要确定该环节输入亮度信号的一些先验知识,如输入亮度信号的概率统计分布规律等。然后,本发明根据已经获知的输入亮度信号的一些先验知识和已知的输出亮度信号的全部知识来确定该环节的Gamma特性参数。在本发明的半盲Gamma特性参数确定方法中,输出亮度信号的全部知识是已知的,但是,本发明并不一定利用输出亮度信号的全部知识。
在确定了Gamma特性参数后,就能够根据Gamma特性参数对Gamma环节进行Gamma校正了。在本发明中,需要确定Gamma特性的环节可以为单个给定的Gamma环节,也可以为多个给定的Gamma环节的级联组合。
下面结合附图对本发明提供的技术方案进行详细描述。
首先,本发明需要获取输入亮度信号的亮度直方图和输出亮度信号的亮度直方图,亮度直方图如附图8所示。图8中,设定图像的亮度为0到255个等级,则不同亮度等级均对应于一个亮度分布概率。
直方图是图像处理技术领域的技术术语,其实,直方图就是一种离散形式的分布概率密度函数。
视频信号是由一帧一帧的图像组成的,输出亮度信号的直方图可以从某一帧图像中获得。从图像中获得亮度信号的直方图的方法属于常规技术,在此不再详细描述。
获得输入亮度信号的亮度直方图的方法有多种,下面主要介绍两种获得输入亮度信号的亮度直方图的方法。
方法一、利用标准参考信号来获得输入亮度信号的亮度直方图的方法。
在视频通信如IPTV中,在正式节目播放之前,首先播放一些标准的测试图像信号,这些测试图像作为输入信号的统计特性是完全可以知道的,这样,视频通信网络中的视频设备可以根据标准的测试图像确定输入亮度信号的亮度直方图。
方法二、利用具有相似统计特性的图像来获得输入亮度信号的亮度直方图。在视频通信如数据会议中,可能会不存在标准的测试图像信号,每一类图像如建筑物、自然场景、人物等,尽管具体的景物对象不同,但是统计特性可以非常接近。比如两个人物的图像,人可以完全不相似,但是,其亮度直方图可以是非常相似的。基于这个原理,目前很多多媒体检索系统,均采用根据直方图的方法进行图像检索,其检索的效果还是很好的。本发明也利用了这个原理来获得输入亮度信号的亮度直方图。这样,就要求在执行Gamma校正的执行点上,存储按照直方图相似原则分类的各种图像的直方图信息,执行Gamma校正的执行点可以在接收到图像后,根据图像的具体内容从其存储的直方图信息中查找对应的直方图信息,并近似将其作为输入信号的直方图信息。
图9给出了基于以上两种获取直方图方法进行Gamma校正的应用场景。
图9中,Gamma校正模块可以存在于视频数据源设备如IPTV的内容源设备中,也可以存在于网络中间设备如IPTV内容数据在传送过程中经过的媒体网关等中间设备中,还可以存在于视频数据目的设备如IPTV用户终端中。在视频数据源设备中可以采用标准的参考信号来确定输入亮度直方图。视频数据源设备、中间设备和视频数据目的设备也可以采用和统计特性相似场景的分类直方图数据库的方法来确定输入亮度直方图。一般来说,虽然以上三种设备中都可能存在Gamma校正模块,即上述三种设备都可以具有Gamma校正功能,但是在具体应用中,只要使用其中某一个设备的校正功能就可以满足Gamma校正的要求了。具体的Gamma校正功能的分配和安排可以通过事先配置的方式来实现,也可以通过通信过程中的协商等方式来实现。
直方图信息和连续的分布概率密度函数之间存在着密切的关系。一般来说,由连续的分布概率密度函数可以直接得到亮度直方图;反过来,由亮度直方图,也可以通过数据插值或者拟合等方法来得到连续的分布概率密度函数。事实上,直方图信息和连续的分布概率密度函数存在严格的比例数量关系。以上关系说明如下。
对于单个Gamma环节或者多个Gamma环节的级联组合来说,输入亮度信号的全体集合是{s(t)|t∈R,0≤s(t)≤1},其中,R表示全体实数集合。也就是说,输入亮度信号的全体集合是全体信号幅值(amplitude)小于等于1的非负值时间信号的集合。这里的输入亮度信号的取值为非负,输入亮度信号的取值是根据物理意义确定的,因为负亮度没有物理意义。任何信号在经过规一化处理之后,一定会满足信号幅值小于等于1的条件。
由于存在随机干扰,所以,这些输入亮度信号可以看成是随机过程。这些输入亮度信号的统计特性可能各不相同,但是,按照信号的统计特性,特别是按照分布概率特性,可以对输入信号进行分类。任何信号作为一个随机过程都有一个分布概率密度函数与之对应,如果随机过程是平稳的(这里的平稳是严格意义上的平稳),那么这个分布概率密度函数和时间无关;如果随机过程不是平稳的,那么这个分布概率密度函数可能和时间有关。因此,一般来说,对于一个随机过程s(t)(t∈R,0≤s(t)≤1)来说,可以用fs(x,t),t∈R表示其分布概率密度函数。如果是严格意义上的平稳的随机过程,则fs(x,t),t∈R和t无关,即分布概率密度函数不随时间变化而变化,此时,fs(x,t)=fs(x)。
信号的规一化处理方法如下如果一个信号s(t)不满足条件t∈R,0≤s(t)≤1,那么,需要对该信号进行规一化处理,使其满足t∈R,0≤s(t)≤1。例如如果信号实际的取值范围是
,则规一化处理后的信号sn(t)为sn(t)=s(t)/Smax(4)公式(4)中的下标n表示英文normalized,意思为规一化。
相应地,如果将信号从规一化的值还原到实际的取值,即对信号进行逆规一化处理,其计算公式如下s(t)=Smaxsn(t) (5)根据分布概率密度函数的定义,分布概率密度函数有如下属性∫-∞+∞fs(x,t)dx=1,]]>对于任何t并且 (6)fs(x,t)≥0,对于任何t而且,对于信号幅值小于等于1的非负值信号,满足fs(x,t)=0,x<0或者x>1 (7)也就是说,信号值大于1或者小于0是不可能的,概率为零。
作为一个自然推论就是∫01fs(x,t)dx=1,]]>对于任何t (8)按照概率密度函数的定义,对于很小的区间长度δ和区间
上一点x0来说,fs(x0,t)δ≈Prob{x0≤s(t)≤x0+δ} (9)或者等效地fs(x0,t)δ≈Prob{x0-12δ≤s(t)≤x0+12δ}---(10)]]>其中符号Prob表示概率(Probability)。
其直观意义是说,在时刻t,亮度信号落在区间[x0,x0+δ]或者 的概率近似等于fs(x0,t)δ。这其实是一种把连续分布概率密度函数变成离散概率密度的方法。由此可知,由连续概率密度函数,通过这样的离散化可以得到信号的亮度直方图。
对于规一化的亮度信号,可以把
区间等分成N个子区间,每个子区间的长度是1/N。第k(k=0,1,2,....,N-1)个子区间是[k/N,(k+1)/N]。如果N足够大,1/N足够小,那么,可以认为1Nfs(2k+12N,t)=Prob{kN≤s(t)≤k+1N},k=0,1,2,....,N-1---(11)]]>于是,可以形成一个概率序列(sequence){1Nfs(2k+12N,t)|k=0,1,2...,N-1}---(12)]]>如果信号还原到其非规一化的信号空间中,如在视频通信中通常亮度信号取0-255的整数,共256级亮度,当然,也可以将亮度信号一般化为2D级亮度的情况,此时,需要将单位区间即
线性映射成集合{0,1,2,3,...,2D-2,2D-1},每个子区间相应扩大2D倍,成为(1/N)2D。于是相应的概率序列变成连续概率密度函数{h(k)|h(k)=1Nfs(2k+12N,t),k=0,1,2...,N-1}---(13)]]>根据公式(8)和公式(10),显然可以得出Σi=02D-1h(i)=1---(14)]]>公式(13)中的这个概率序列就叫做亮度信号s(t)的直方图。
从上述推导中可以明显看出直方图是可以由信号亮度的连续分布概率密度函数直接得到的,反过来,亮度信号的连续分布概率密度函数也可以由直方图经过数据插值、拟合等处理后得到。
下面给出一个直方图的具体例子。当信号的亮度包括256亮度级时,概率序列和直方图中具体数值的对应关系为h(0)=0h(1)=0....
h(64)=0.005h(65)=0.006.....
h(190)=0.006h(191)=0.005h(192)=0.001h(193)=0…h(255)=0下面用函数y=g(x;p),p=[p1,p2,...,pM]T来表示Gamma特性参数未知的Gamma环节的Gamma特性,这里的Gamma环节包括单个Gamma环节或者多个Gamma环节的级联组合的情况。上述Gamma特性的表示方式中,p=[p1,p2,...,pM]T是一个参数向量,一般情况下,参数向量由M个参数组成。这些参数的全部或者部分是需要确定的。因此,按照这个很一般的形式,Gamma特性函数几乎可以是任何形式的函数,只要满足函数是连续的条件即可,而且,一般来说,Gamma特性函数是光滑可导的,至少是分段光滑可导的,因此,假设Gamma特性函数关于变量x的导数存在是合理的。Gamma特性函数的导数可以用如下符号表示d(x;p)=dg(x;p)dx---(15)]]>并且,Gamma特性函数,还应该满足g(1;p)=1 (16)一般来说,Gamma特性函数可以用如下三种常用的方式来表示方式一、幂函数y=g(x;p)=p1xp2+p3,p=[p1,p2,p3]T---(17)]]>方式二、多项式函数y=g(x;p)=p1xK+p2xK-1+....+pKx+PK+1,p=[p1,p2,p3,...,PK+1]T(18)方式三、分段线性函数
在公式(19)中,为了保证Gamma特性函数连续性这个约束条件,Gamma特性函数需满足约束条件akck+1+bk=ak+1ck+1+bk+1,k=1,2,...,K-1。
如果用e(t)和r(t)分别表示输入亮度信号和输出亮度信号,那么,e(t)和r(t)各自对应的分布概率密度函数是fe(x,t)和fr(x,t),并且,e(t)、r(t)和Gamma特性函数之间存在如下关系r(t)=g(e(t);p),p=[p1,p2,p3]T。
根据概率理论,可以推导出d(e;p)fr(r)=fe(e),其中r=g(e,p),e∈
(20)推导出公式(20)的具体过程可以参见常用的概率书籍,在本实施例中不再详细描述。
本发明可以根据公式(20)获得Gamma特性函数的导函数d(x;p),从而可以求出参数向量p。当然如果Gamma特性函数含有常数项,那么,对应的常数作为一个参数是不出现在d(x;p)中的,通过确定d(x;p)的参数之后,可以再通过关系(16)作为一个约束条件来确定常数项。其实不论如何,因为公式(16)的存在,使得所有的参数之间满足一个约束关系,这些参数如果是M个,那么,只有M-1个是独立的,我们只要确定其中的任意M-1个参数就可以了,剩下的一个参数用公式(16)得到的方程来求解。
从公式(20)可以看出,公式(20)和时间变量t无关。其实,Gamma特性函数本身和时间变量无关,因此,只要在一段相当长的时间内测定一次Gamma特性参数,在整个通信过程中就可以一直使用这组Gamma特性参数,如在IPTV中,一个节目的Gamma特性参数可以认为是相同的,这样,最多在每个节目开始的时候测量一次Gamma特性参数就可以了。
设定本发明获得的输入亮度信号和输出亮度信号的直方图为{he(k)|k=0,1,2,....,N-1},{hr(k)|k=0,1,2,....,N-1}。也就是说,每个直方图中包含有N个项,在直方图术语中,每一个项叫做“柱”(bin)。
根据公式(20)和公式(13),可以容易的得到Gamma特性函数的具体表现形式为d(2k+12N;p)(λhr(k′)+(1-λ)hr(k′+1))=he(k);---(21)]]>其中,k′=[g(2k+12N,p)],k=0,1,2,...,N-1;]]>λ=k′+1N-g(2k+12N,p)]]>将公式(21)写成一体的形式,则可以变为公式(22)d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))=he(k);]]>其中k=0,1,2,...,N-1,N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,he(k)为输入亮度信号的亮度直方图,he(k)为输出亮度信号的亮度直方图。
在公式(22)中,函数[.]表示取整函数,即对于实数x,[x]表示不大于x的最大整数,比如[2.4]=[2.6]=2。
通过上述描述获得了Gamma特性函数的形式,只是参数向量p需要确定,如果有N个方程,只要N≥M-1,这里的M是需要确定参数的个数,那么,一定可以通过这组方程求解出唯一的解作为其中M-1个参数的确定值,然后,再利用公式(16)推导出的方程确定剩下的一个参数向量,从而确定出了Gamma特性函数的所有参数。在一般情况下,条件N≥M-1都能满足。
根据公式(22)获得参数向量p的方法有多种,下面主要介绍两种确定参数向量p的方法一、直接求解方程法。
由公式(22)获得N个方程,从中任意选择M-1个方程,形成方程组联立求解。一般来说,这组方程是非线性的,而且是超越(Transcendental)的,如对于幂函数形式的Gamma特性函数等。因此不存在解析解(closed-form solution)。需要用数值解法(numerical solution method)。关于方程组的数值解法,属于常规技术,在本实施例中不在详细描述。
方法二、非线性函数优化方法。
根据公式(22)构造一个代价函数J(p)=Σk=0N-1(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))-he(k))2(d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+---(23)]]>其中N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,M为伽玛特性参数的个数。
显然,对于Gamma特性函数的真正的参数向量ptrue,在理论上应该使得J(ptrue)=0,由于参数向量ptrue满足公式(22)中的每一个方程,因此,公式(23)中的每个求和项都是零。但是在实际应用中,因为存在误差和近似如获得直方图过程中的近似等,所以,公式(23)中的每个求和项不会都等于零,但是,应该是一个很小的数值,并且应该满足如下条件对于任何p∈RM,都有J(p)≥J(ptrue)。也就是说,ptrue是函数J(p)的全局最小点(global minimal point)。
因此,在非线性函数优化方法中,确定参数向量p的问题就转化成对代价函数J(p)求其全局最小点的数学问题。
由于Gamma特性函数的函数形式是非线性并且是超越的,并且因为存在取整函数[.],因此,J(p)不一定是可导的,甚至不一定是连续的,这样,不能用涉及导数的优化方法如梯度方法,共轭梯度方法等来确定参数向量p。
在非线性函数优化方法中,可以采用如下三种方式来确定参数向量p。
方式(1)、组合优化方法。
通过组合优化方法来确定参数向量p的具体过程属于常规技术,在本实施例中不再详细描述。
方式(2)、神经网络方法。
通过神经网络方法来确定参数向量p的具体过程属于常规技术,在本实施例中不再详细描述。
方式(3)、蛮力搜索方法(Brutal Force Search Method)。所谓蛮力搜索,顾名思义,就是穷尽搜索所有的可能性。
对于参数取离散值的情况,则参数所有可能取值的集合是有限集合,那么,逐一搜索集合中的每个点,就能够找到使得J(p)最小的点,这样,就找到了全局最小点ptrue。但是,这种情况很少,在绝大多数情况下,参数是取连续值的,因此,参数所有可能取值的集合是无限集合,无法真正进行穷尽搜索。
对于参数取连续值的情况,蛮力搜索方法是将参数空间分成多个小的超立方体(Hypercube),然后,在每个超立方体中取一个点作为采样点,如超立方体的几何中心点等,最后,计算代价函数在各个超立方体的采样点上的函数值,找到使得代价函数最小的超立方体的采样点,将该采样点作为全局最小点。
下面针对参数取连续值的情况、结合附图对利用蛮力搜索方法确定Gamma特性函数的参数的实现过程进行详细描述。
蛮力搜索方法虽然有点“笨”,但是,蛮力搜索方法在工程实践中有着广泛用途,比如密码破译等,最有效的方法还是蛮力搜索方法。而且,对于数学优化问题,现有的技术除了蛮力搜索之外,都存在陷入局部极小(local minimalpoint)的问题,但是蛮力搜索法不存在该问题。
在本发明中,可以利用关于参数的先验知识,找到合适的起始搜索点,这样,可以大大降低需要搜索的次数,从而提高蛮力搜索方法的效率。
蛮力搜索方法的原理如附图10示。
图10中,M-1个参数所有可能取值的集合构成了参数空间(Parameter Space,简称PS),参数空间是M-1维欧式空间RM-1的一个子集。在确定了参数空间后,蛮力搜索的实现方法包括如下步骤步骤一、超立方体划分。
将PS划分成多个M-1超立方体(Hypercube),附图8中,ABCDEFGH,8个点组成一个超立方体。由于每个参数的取值范围大小不同,所以,超立方体的每个边长也不同。
设定第k(k=1,2,...,M-1)个参数pk的取值范围是[mink,maxk],对第k个维度进行Pk等分,从而该维度上每个超立方体的边长是Δk=maxk-minkPk---(24)]]>从而每个超立方体的体积为V=Πk=1M-1Δk=Πk=1M-1(maxk-minkPk)---(25)]]>因此,总的超立方体个数大于T=Πk=1M-1Pk.]]>这个数值是可能的最大值,如果PS的形状不是一个超立方体,那么其包含的长立方体个数可能远远小于T=Πk=1M-1Pk.]]>对于每个超立方体用指标向量I=[i1,i2,....,iM-1]T来表示,其中ik(k=1,2,...,M-1,ik=1,2,3,....,Pk)表示在第k个维度上,该超立方体是第ik个。
因此,对于第I=[i1,i2,....,iM-1]T个超立方体来说,其在各个维度上的坐标范围是[mink+(ik-1)Δk,mink+ikΔk] (26)该超立方体的几何中心QI的坐标是[ain1+(i1-1/2)Δ1,min2+(i2-1/2)Δ2,..............,minM-1+(iM-1-1/2)ΔM-1]T(27)步骤二、选取初始搜索点。
一般来说,对于某个参数pk,都存在一个比较合理的值,可以作为初始值。如Gamma特性函数采用幂函数形式时,对于视频输入设备来说,Gamma参数的取值一般在2.2左右,因为工业标准要求是2.2,由于制造技术和产品品质的原因,Gamma参数可能会正负偏离2.2,但是在多数情况下,比较接近2.2。这样的话,如果以2.2作为初始搜索点开始搜索,由于2.2比较接近真实值,则找到真实值需要尝试的次数就比较少。同样的道理,可以为每个参数pk(k=1,2,...,M-1)找到其合适的初始值pint k,那么这些初始值形成一个向量,就是参数向量p的初始值pint=[pintk,pintk,pintk,....,pint(M-1)]T.]]>步骤三,根据选取的初始搜索点开始搜索。
首先判断pint落在哪个超立方体中,通过比较坐标等方法可以判定出pint落在哪个超立方体中。比较坐标方法的具体实现过程为设定该超立方体的坐标是Iint=[iint1,iint2,....,iint(M-1)]T,判断条件是对于k=0,1,2.....M-1时,当且仅当下述公式(28)均成立,则确定出pint落在坐标为Iint=[iint1,iint2,....,iint(M-1)]T的超立方体中。
mink+(iintk-1)Δk≤pintk<mink+iintkΔk(28)在确定了初始搜索点所在的超立方体后,从该超立方体开始搜索。根据公式(22)计算J(pint),如果J(pint)能够使公式(29)成立,或者搜索过的超立方体满足预定条件,如搜索过的超立方体的数量达到预定数值等,那么,整个搜索过程结束,到步骤五。此时,ptrue=pint;其中,ptrue为全局最优点,即最终获得的Gamma特性函数的参数向量。
J(pint)≤Jthreshold(29)其中,Jthreshold是一个预先给定的门限值。
如果J(pint)不能够使公式(29)成立,则到步骤四。
步骤四、继续搜索。
在步骤四中的继续搜索可以是分层进行的。包围在初始超立方体外部的超立方体可以为一层或多层。
在二维情况下,初始超立方体和其外围多层超立方体如附图11所示。
图11中,中间灰色的正方形为初始超立方体,与初始超立方体的边邻接的立方体为第一层超立方体,与第一层超立方体的边邻接的立方体为第二层超立方体,与第二层超立方体的边廉价的立方体为第三层超立方体。
本发明的分层搜索方法为逐次搜索初始超立方体之外的每一层中的每个超立方体,在每一层超立方体的搜索中,应按照预定顺序遍历该层中的每一个超立方体。预定顺序可以是多种多样的,本发明不限制预定顺序的形式,只要能够遍历一层中的超立方体就可以。
在某一层超立方体搜索过程中,当根据预定顺序搜索到某个超立方体时,需要按照公式(27)计算该超立方体的几何中心Q的坐标,然后,计算函数值J(Q),如果J(Q)能够使公式(30)成立,或者搜索过的超立方体的数量达到预定数值,那么,整个搜索过程结束,到步骤五。此时,ptrue=Q;J(Q)≤Jthreshold(30)如果J(Q)不能够使公式(30)成立,而且,搜索过的超立方体的数量也没有达到预定数值,继续根据预定顺序在本层超立方体中搜索。如果本层的超立方体搜索完成、且本层各超立方体的J(Q)均不能够使公式(30)成立,则继续搜索下一层超立方体。
当搜索完第L层后,不论是否找到满足条件(30)的超立方体的几何中心Q,搜索过程都将结束,此时,应将搜索到的最小的J(Q)中的Q作为ptrue。到步骤五。这里的L是预先给定的一个门限值,表示最多需要搜索的层数。
步骤五,搜索过程结束。
上述搜索方法还可以应用于按照从粗到细的搜索过程中,达到最快最好的搜索效果。
从粗到细的搜索过程如附图12所示。
图12中,首先,按照上述步骤一到步骤五进行第一次搜索。第一次搜索可以看成是粗搜索,这样,可以将参数空间中的超立方体的边长设置的大一些,这样,参数空间中超立方体的数量较少。在从粗到细的搜索过程中,步骤三中搜索过的超立方体满足的预定条件可以为划分粒度是否粗于预定划分粒度。第一次搜索如果找到了满足条件(29)或者(30)的超立方体,则搜索过程结束。第一次搜索过程可以采用分层搜索的方法。
如果在第一次搜索过程中没有搜索到满足条件(29)或者(30)的超立方体,而且,超立方体的粒度划分还没有达到预定粒度,则在第一次搜索到的最小的J(Q)对应的超立方体中进行第二次较细的搜索,此时,应将第一次搜索到的最小的J(Q)对应的超立方体当成新的整个参数空间。此时的新的参数空间中的每个超立方体的边长变小了,重复上述步骤一至步骤五的搜索过程。同样,在第二次较细的搜索过程中,如果找到了满足条件(29)或者(30)的超立方体,第二次搜索过程结束。第二次搜索过程可以采用分层搜索的方法。
如果在第二次较细的搜索过程中,没有找到满足条件(29)或者(30)的超立方体,而且,超立方体的粒度划分还没有达到预定粒度,则在第二次搜索到的最小的J(Q)对应的超立方体中进行第三次更细的搜索,依此类推,将最后一次精细搜索找到的满足条件(29)或者(30)的超立方体的几何中心作为全局最优点ptrue,即Gamma特性函数的参数向量。
如果超立方体的粒度划分达到了预定粒度,但是仍然没有找到满足条件(29)或(30)的超立方体,则搜索过程结束。此时,应将查找到的最小的J(Q)对应的超立方体的集合中心作为全局最优点ptrue,即Gamma特性函数的参数向量。
在上述从粗到细的每一次的搜索过程中都可以采用分层搜索的方法。
在通过上述方法确定了Gamma特性函数的参数向量后,就能够对Gamma环节进行Gamma校正了。这里,需要进行Gamma校正的环节可以为视频数据源设备,也可以为视频通信网络中的中间设备,还可以为视频数据目的设备。
本发明为多媒体通信系统提供了一种易于实现的Gamma校正方法,本发明提供的Gamma校正方法具有很高的应用可行性,从而大大拓宽了Gamma校正的应用范围,特别能够针对IPTV、协作数据会议、广泛使用低端视频输入设备的公众视频通信提供了很好的Gamma校正功能,大大提高用户体验和服务质量,进一步提升上述业务的竞争力,为电信运营商、服务提供商和设备厂商带来巨大的经济效益。
本发明提供的视频通信伽玛特性的校正装置,包括获取直方图模块、存储模块、计算模块和伽玛校正模块。
获取直方图模块主要用于获取输入、输出图像的亮度直方图。获取直方图模块可以利用标准参考信号来获得输入亮度信号的亮度直方图,也可以利用具有相似统计特性的图像来获得输入亮度信号的亮度直方图。具体如上述方法中的描述。
存储模块主要用于存储输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系,如存储模块存储公式(21)或者公式(22)。
计算模块主要用于根据获取直方图模块获得的输入、输出图像的亮度直方图对存储模块中的所述数学关系进行求解,以获得伽玛特性参数。计算模块可以采用直接求解方程法、非线性函数优化方法如组合优化方法、神经网络方法、蛮力搜索方法等来对公式(21)或者(22)求解。在使用蛮力搜索方法时,计算模块首先构造代价函数(23),然后,采用逐层搜索、从粗到细搜索等方法来获得伽玛特性参数。计算模块将其获得的伽玛特性参数传输至伽玛校正模块。具体实现过程如上述方法中的描述。
伽玛校正模块主要用于根据其接收到的伽玛特性参数对伽玛环节进行伽玛校正。这里的伽玛环节可以为单个给定的伽玛环节,也可以为多个伽玛环节的级联。
本发明提供的装置位于视频设备中,如位于视频数据源设备中、位于视频通信网络的中间设备中,再如位于视频数据目的设备中。
虽然通过实施例描绘了本发明,本领域普通技术人员知道,本发明有许多变形和变化而不脱离本发明的精神,本发明的申请文件的权利要求包括这些变形和变化。
权利要求
1.一种视频通信伽玛特性的校正方法,其特征在于,包括a、获取输入、输出亮度信号的亮度直方图;b、建立输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系;c、根据所述输入、输出图像的亮度直方图对所述数学关系进行求解,以确定伽玛特性参数;d、根据所述伽玛特性参数对伽玛环节进行伽玛校正。
2.如权利要求1所述的方法,其特征在于,所述伽玛环节为单个给定伽玛环节、或者多个给定伽玛环节的级联组合。
3.如权利要求1所述的方法,其特征在于,所述步骤a中获取输入图像的亮度直方图的步骤包括根据已知统计特性的测试图像获取输入图像的亮度直方图;或者根据图像统计特性相似性从预先存储的直方图信息中选取输入图像的亮度直方图。
4.如权利要求1所述的方法,其特征在于,所述步骤b包括输入、输出亮度信号各自的亮度分布概率密度函数fe(x,t)、fr(x,t)、以及Gamma特性函数之间的数学关系为d(e;p)fr(r)=fe(e);其中r=g(e,p),e∈
;亮度分布概率密度函数和亮度直方图之间的数学关系为1Nfs(2k+12N,t)=Prob{kN≤s(t)≤k+1N},]]>其中k=0,1,2,...,N-1;将输入、输出信号各自的亮度分布概率密度函数和Gamma特性函数之间的数学关系转化为输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))=he(k);]]>其中k=0,1,2,...,N-1,N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,he(k)为输入亮度信号的亮度直方图,he(k)为输出亮度信号的亮度直方图,函数[.]表示取整函数。
5.如权利要求1至4中任一权利要求所述的方法,其特征在于,所述步骤c包括利用直接求解方程法、或者非线性函数优化方法进行求解,以确定伽玛特性参数。
6.如权利要求5所述的方法,其特征在于,非线性函数优化方法包括根据输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系构造代价函数J(p)=Σk=0N-1(d(2k+12N;p)((k′+1N-g(2k+12N,p))hr([g(2k+12N,p)])+(g(2k+12N,p)-k′+1N)hr([g(2k+12N,p)]+1))-he(k))2;]]>其中N为输入信号亮度直方图所包含的项数和输出信号亮度直方图所包含的项数之中的最小值,M为伽玛特性参数的个数;在M个参数中任意选取M-1个参数,将参数空间的维度降低为M-1维;确定M-1个参数的方法为确定参数向量ptrue,使对于任意p∈RM-1,关系J(p)>=J(pture)成立,即寻找到代价函数的全局最小点,参数向量ptrue即为M-1个参数;利用关系g(1;p)=1结合ptrue,确定剩余一个伽玛特性参数。
7.如权利要求6所述的方法,其特征在于,确定ptrue的方法包括组合优化方法、神经网络方法、蛮力搜索方法。
8.如权利要求7所述的方法,其特征在于,所述蛮力搜索方法包括步骤将维度降低一维的伽玛特性参数空间划分为多个超立方体;选取初始搜索点,并根据预定顺序从该初始搜索点所在的超立方体开始进行遍历搜索;计算搜索过程中进入的每个超立方体的几何中心坐标Q,并根据所述代价函数的表达式以及搜索过程中进入的超立方体计算J(Q);如果J(Q)小于等于预定门限,或者搜索过的超立方体满足预定条件,则将本次搜索进入的超立方体的几何中心坐标Q作为ptrue,搜索过程结束;否则,继续搜索过程。
9.如权利要求8所述的方法,其特征在于,所述初始搜索点根据实际应用中伽玛特性参数的经验数值来设置。
10.如权利要求8所述的方法,其特征在于,所述遍历搜索包括以初始搜索点所在的超立方体作为初始超立方体,将包围初始超立方体的超立方体按照距离初始超立方体的远近划分为依次包裹前一层的多层超立方体阵列,并逐层进行搜索。
11.如权利要求8所述的方法,其特征在于在蛮力搜索方法中,维度降低一维的伽玛特性参数空间被划分为多个粗粒度的超立方体,所述遍历搜索包括以初始搜索点所在的超立方体作为初始超立方体,将包围初始超立方体的超立方体按照距离初始超立方体的远近划分为依次包裹前一层的多层超立方体阵列,并逐层进行搜索;将搜索到的满足条件的超立方体作为新的伽玛特性参数空间,并将其划分为更细粒度的超立方体,依次类推,进行从粗到细的逐层搜索。
12.一种视频通信伽玛特性的校正装置,其特征在于,所述装置包括获取直方图模块、存储模块、计算模块和伽玛校正模块;获取直方图模块用于获取输入、输出图像的亮度直方图;存储模块用于存储输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系;计算模块用于根据输入、输出图像的亮度直方图对存储模块中的所述数学关系进行求解,以确定伽玛特性参数;伽玛校正模块用于根据所述伽玛特性参数对伽玛环节进行伽玛校正。
13.如权利要求12所述的装置,其特征在于,所述装置位于视频数据源设备中、和/或位于视频通信网络的中间设备中、和/或位于视频数据目的设备。
全文摘要
本发明提供一种伽玛特性的校正方法和装置,其核心为根据输入、输出亮度信号的亮度直方图对输入、输出信号各自的亮度直方图和Gamma特性函数之间的数学关系进行求解,以确定伽玛特性参数,并对伽玛环节进行伽玛校正。本发明提供的Gamma校正方法仅需要输入亮度信号的直方图和输出亮度信号的直方图即可,这样,本发明的Gamma参数确定方法可以称为半盲测量方法;由于本发明不需要输入亮度信号的全部知识,因此,本发明的Gamma校正方法具有很高的应用可行性;本发明的伽玛校正方法特别适用于IPTV、协作数据会议、以及广泛使用低端视频输入设备的公众视频通信等多媒体视频通信业务;从而达到了提高Gamma校正易用性,拓宽Gamma校正的应用范围,提高用户体验和服务质量的目的。
文档编号H04N9/69GK101052132SQ20061009290
公开日2007年10月10日 申请日期2006年6月9日 优先权日2006年6月9日
发明者罗忠 申请人:华为技术有限公司
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1