CN105976337B - 一种基于中值引导滤波的图像去雾方法 - Google Patents

一种基于中值引导滤波的图像去雾方法 Download PDF

Info

Publication number
CN105976337B
CN105976337B CN201610307360.5A CN201610307360A CN105976337B CN 105976337 B CN105976337 B CN 105976337B CN 201610307360 A CN201610307360 A CN 201610307360A CN 105976337 B CN105976337 B CN 105976337B
Authority
CN
China
Prior art keywords
image
transmissivity
value
intermediate value
guiding filtering
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610307360.5A
Other languages
English (en)
Other versions
CN105976337A (zh
Inventor
黄鹤
宋京
王会峰
汪贵平
张弢
雷旭
郭璐
黄莺
许哲
易盟
袁东亮
陈志强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhilian cloud big data technology Nanjing Co.,Ltd.
Original Assignee
Changan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changan University filed Critical Changan University
Priority to CN201610307360.5A priority Critical patent/CN105976337B/zh
Publication of CN105976337A publication Critical patent/CN105976337A/zh
Application granted granted Critical
Publication of CN105976337B publication Critical patent/CN105976337B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/73—Deblurring; Sharpening
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10004—Still image; Photographic image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G06T2207/20032—Median filtering

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于中值引导滤波的图像去雾方法，获取雾霾天气下的图像；将获得的含雾图像转换到RGB颜色空间，并且获得各个像素点的R、G、B三个通道的值；对各个像素点的R、G、B三个通道值进行比较，获得最小值，从而得到图像的初始暗通道图像；根据获得的初始暗通道图像，计算大气光值A；根据透射率估计公式，得到粗估计透射率；通过中值引导滤波对粗估计透射率进行滤波，获得精准的透射率；将获得的图像以及获得的大气光值A和获得的透射率代入大气光成像模型方程反向求解即获得去雾图像。本发明在实现对透射率进行精确优化的同时还对暗通道先验对大气光估计的过程进行了优化，同时具有更快的运算速度，更加具有实用性。

Description

一种基于中值引导滤波的图像去雾方法
技术领域
本发明属于图像处理技术领域，具体涉及一种基于中值引导滤波的图像去雾方法。
背景技术
雾霾天气下，光线受空气中大量的小水珠、固态微粒的影响而产生散射和折射，使得图像采集设备采集到的图像对比度降低、色彩饱和度下降，甚至色调偏移。这使得我们提取图像中的有用信息的难度大大提高，而且这样的图像往往难以满足机器视觉或者其他设备的要求。在实际中，图像去雾在交通信息的监控、民用的安防、军用的遥感侦查等方面都有着不可或缺的作用。因此，对雾霾天气下采集到的图像进行去雾处理，提高图像的质量是十分有意义的。
目前，数字图像去雾研究可大致分为两大类：一类是基于图像处理的图像增强方法。这类方法是为了根据需要，不考虑图像退化的过程，有针对的、有选择性的突出某些需要的信息，抑制不需要的或者干扰信息，从而改善图像的视觉效果或者使图像更能满足图像处理设备的要求，实现去雾的目的。例如彩色增强、对比度增强、白平衡，但这类方法往往计算量较大、实时性较差。其本质上并没有实现去雾，没有考虑图像退化的过程，处理后的图像容易失真，而且算法是针对某一图像而编写的普适性较差。
另一类是基于物理模型的图像复原方法。这类方法考虑到空气中介质对成像的影响，即入射光的衰减和大气光成像过程来建立图像退化模型，利用先验知识或假设某些条件，来反演图像的退化过程，从而实现图像的复原。与前一类方法相比，后者对图像的退化过程进行了建模，通过求解退化过程，反解出原始图像，得到的复原图像更加符合实际，具有极大的优越性。但是实际中得到精确的图像退化模型是十分困难的，这给基于物理模型的图像复原带来了极大的难度。
近年来，基于物理模型的图像去雾算法，由于其内在的优越性，取得了极大的进展。何凯明提出的“基于暗通道优先的单幅图像去雾”方法取得了很好的去雾效果，被广泛认可。但是该方法在优化透射率时采用了软抠图的方法，使得该算法的复杂度很高。后来有人改进为使用联合双边滤波来代替，使得算法的效率有所提升，但效果依然十分有限。最后，何凯明等人又提出了“引导滤波”，极大的提高了“基于暗通道优先的单幅图像去雾”算法的效率。但是由于暗通道先验的对连续浅色区域不适用的缺点及算法仍具有较高的复杂度，所以仍有改进的空间。
在对雾霾图像进行去雾处理时，基于图像增强的方法由于其并没有实质的去雾过程，因此得到的去雾图像往往失真严重，而且这类方法对不同的图像很难取得同样好的效果，不具有普适性。基于物理模型的去雾方法，以传统的暗通道先验去雾的引导滤波方法为例，由于传统的暗通道先验对于明亮区域无效，使得处理得到的图像存在图像失真，此外已有的引导滤波算法计算效率不是很高。
发明内容
本发明的目的在于提供一种基于中值引导滤波的图像去雾方法，以克服上述现有技术存在的缺陷，本发明在实现对透射率进行精确优化的同时还对暗通道先验对大气光估计的过程进行了优化，使得获得的大气光更加精准，同时具有更快的运算速度，算法更加简单，更加具有实用性。
为达到上述目的，本发明采用如下技术方案：
一种基于中值引导滤波的图像去雾方法，包括以下步骤：
步骤1：获取雾霾天气下的图像I；
步骤2：将步骤1中获得的含雾图像I转换到RGB颜色空间，并且获得各个像素点的R、G、B三个通道的值；
步骤3：对步骤2中的各个像素点的R、G、B三个通道值进行比较，获得最小值，从而得到图像I的初始暗通道图像；
步骤4：根据步骤3中获得的初始暗通道图像，取亮度前0.1％的像素，计算大气光值A；
步骤5：根据透射率估计公式，得到粗估计透射率t；
步骤6：通过中值引导滤波对步骤5中的粗估计透射率t进行滤波，获得精准透射率t_d；
步骤7：将步骤1获得的图像I以及步骤4中获得的大气光值A和步骤 6获得的透射率t_d代入大气光成像模型方程反向求解即获得去雾图像J。
进一步地，步骤4中计算大气光值A的方法为：取步骤3中获得的初始暗通道图像中亮度前0.1％的像素，并将这些像素的亮度和位置保存，然后从含雾图像I中取出这些像素点，求出这些像素点在R、G、B三个通道中的中值，作为大气光值A，所述大气光值A为1×3数组。
进一步地，步骤5中的透射率估计公式为：
其中，y为局部区域Ω(x)内的一像素点，Ic(y)表示这一像素点的RGB三个通道的值，Ac为步骤4得到的RGB三个通道的大气光值，ω为透射率调整系数。
进一步地，ω＝0.8。
进一步地，步骤6中中值引导滤波公式为：
其中，Ii为引导图像中的某一像素点，且引导图像为由含雾图像转换的灰度图像，qi为引导滤波器输出图像中的某一像素点，且引导滤波器输出图像为精准透射率，ak和bk为滤波窗口wk中的一组线性系数， 其中为引导图像在窗口wk中的中值，为引导滤波器输入图像即粗估计透射率在窗口wk中的中值，ε为滤波调整系数。
进一步地，ε＝10-6。
进一步地，wk的大小为3×3。
进一步地，步骤7中大气光成像模型方程为：
其中，I(x,y)为输入的含雾图像中的某一像素，A为大气光值，J(x，y) 为去雾后的该像素点的值，t_d为引导滤波后的精准透射率，t0为防止t_d过小的限幅常数。
进一步地，t0＝0.1。
与现有技术相比，本发明具有以下有益的技术效果：
本发明提出的中值引导滤波算法使用了含雾图像的灰度图做为引导图像，在计算引导滤波器的系数ak、bk时使用了新的计算方式，避免了已有方法需要对引导图像和滤波图像求方差、求协方差，同时对滤波框内临近像素点采用同一ak、bk，减少了引导滤波器的计算量，另外本发明中提出的基于中值引导滤波的快速图像去雾方法对暗通道大气光估计进行了优化，计算的大气光更加准确，中值引导滤波计算方法更加简单，计算复杂度更低、速度更快。
附图说明
图1是本发明的流程示意图；
图2是本发明与传统的引导滤波算法的对比图，其中，其中图(a)、(d)、 (g)为含雾的图像，图(b)、(e)、(h)为采用传统的引导滤波算法去雾后的图像，图(c)、(f)、(i)为采用本发明提出的中值引导滤波算法去雾后的图像。
具体实施方式
下面结合附图对本发明作进一步详细描述：
参见图1和图2，本发明提供了一种基于中值引导滤波的快速图像去雾方法，采用暗通道先验原理对图像进行大气光估计，再对图像的透射率进行估计，而后用中值引导滤波对透射率进行优化，最后根据大气光模型求解出去雾图像。具体思路是，针对传统的去雾算法计算复杂度高，去雾后的图像存在失真的问题，在优化了大气光估计的同时提出了中值引导滤波算法，新的算法在估计大气光时提取了图像中亮度前0.1％的点，并取得了这些点的R、 G、B三个通道的均值，对大气光值A的三个颜色通道做出了不同的估计。中值引导滤波算法重新定义了引导滤波器方程的系数a、b。此外求解了图像中处理框中像素的中值，用这一中值来代替处理框中的其他像素，在取得与原有算法同样的效果的同时简化了计算，提高了算法的时效性。
具体步骤如下：
步骤1、获得雾霾天气下的图像：利用图像采集设备，获得要处理的雾霾天气下的图像I。
步骤2、计算步骤1中获得的雾霾图像I的尺寸m×n，并计算雾霾图像的各个像素点上各个通道的值，即R(i，j)、G(i，j)、B(i，j)。
步骤3、对获得的各个像素点R、G、B值进行比较，取最小值作为该点像素的暗通道亮度值，即
步骤4、根据步骤2中获得的雾霾图像的暗通道图像I_dark，提取亮度前0.1％的像素，将这些像素存入向量DarkVec中，同时将这些像素点在原图像中的位置存入Location中。而后根据DarkVec和Location向量，从含雾图像中取出这些像素点，分别计算出这些像素点在R、G、B三个通道中的中值，将得到的R、G、B三个通道的中值，作为估计的大气光值A在三个通道中估计值。
步骤5、按照透射率粗估计公式此处ω取ω＝0.8，将步骤3中获得的大气光值A及含雾图像I代入上式中，获得粗估计的透射率t(x)，其中，y为局部区域Ω(x)内的一像素点，Ic(y)表示这一像素点的RGB三个通道的值，Ac为步骤4得到的RGB三个通道的大气光值，ω为透射率调整系数。
步骤6、根据qi＝akIi+bk，取ε为10-6，将含雾图像转换为灰度图作为引导图像，将步骤4中计算的粗估计透射率t作为引导滤波器输入图像，在计算ak,bk时，根据中值引导滤波器定义，选择的像素窗口wk的大小为3×3。从而得到中值引导滤波器的输出图像：精准透射率t_d，其中，Ii为引导图像中的某一像素点，且引导图像为由含雾图像转换的灰度图像，qi为引导滤波器输出图像中的某一像素点，且引导滤波器输出图像为精准透射率，ak和bk为滤波窗口wk中的一组线性系数，其中为引导图像在窗口wk中的中值，为引导滤波器输入图像即粗估计透射率在窗口wk中的中值，ε为滤波调整系数。
步骤7、根据步骤3得到的大气光估计值和步骤5得到的精细化后的透射率t_d，代入根据大气光成像求解得到的去雾图像方程为了防止t_d过小，造成去雾图像失真，取t0＝0.1，其中，I(x,y)为输入的含雾图像中的某一像素，A为大气光值，J(x，y)为去雾后的该像素点的值，t_d为引导滤波后的精准透射率，t0为防止t_d过小的限幅常数。
图1为算法的流程图。图2为不同算法的处理效果图。从图2可以看出，图2中(a)、(d)、(g)为原始雾霾天气下的图像，我们分别使用已有的传统的引导滤波去雾方法和本文中提出的中值引导滤波去雾方法对图像进行去雾处理，其中如(b)、(e)、(h)为使用已有的引导滤波去雾的结果，可以看到处理后的图像有着不错的去雾效果，但是在天空等区域出现了颜色失真的情况，图(c)、(f)、(i)为使用本发明方法处理后的结果，可以看到图像色彩比较符合原图，同时有着相当或者由于已有的引导滤波算法的效果。
表2是不同算法对不同尺寸的图像处理用时的对比表。从表2中可以看出，对于尺寸比较大的图像，本文提出的基于中值引导滤波的快速图像去雾方法的时效性优于已有的引导滤波去雾算法。
表2 不同算法对不同尺寸的图像处理用时的对比表

Claims (8)

1.一种基于中值引导滤波的图像去雾方法，其特征在于，包括以下步骤：
步骤1：获取雾霾天气下的图像I；
步骤2：将步骤1中获得的含雾图像I转换到RGB颜色空间，并且获得各个像素点的R、G、B三个通道的值；
步骤3：对步骤2中的各个像素点的R、G、B三个通道值进行比较，获得最小值，从而得到图像I的初始暗通道图像；
步骤4：根据步骤3中获得的初始暗通道图像，取亮度前0.1％的像素，计算大气光值A；
步骤5：根据透射率估计公式，得到粗估计透射率t；
步骤6：通过中值引导滤波对步骤5中的粗估计透射率t进行滤波，获得精准透射率t_d；其中，中值引导滤波公式为：
其中，Ii为引导图像中的某一像素点，且引导图像为由含雾图像转换的灰度图像，qi为引导滤波器输出图像中的某一像素点，且引导滤波器输出图像为精准透射率，ak和bk为滤波窗口wk中的一组线性系数， 其中为引导图像在窗口wk中的中值，为引导滤波器输入图像即粗估计透射率在窗口wk中的中值，ε为滤波调整系数；
步骤7：将步骤1获得的图像I以及步骤4中获得的大气光值A和步骤6获得的透射率t_d代入大气光成像模型方程反向求解即获得去雾图像J。
2.根据权利要求1所述的一种基于中值引导滤波的图像去雾方法，其特征在于，步骤4中计算大气光值A的方法为：取步骤3中获得的初始暗通道图像中亮度前0.1％的像素，并将这些像素的亮度和位置保存，然后从含雾图像I中取出这些像素点，求出这些像素点在R、G、B三个通道中的中值，作为大气光值A，所述大气光值A为1×3数组。
3.根据权利要求2所述的一种基于中值引导滤波的图像去雾方法，其特征在于，步骤5中的透射率估计公式为：
其中，y为局部区域Ω(x)内的一像素点，Ic(y)表示这一像素点的RGB三个通道的值，Ac为步骤4得到的RGB三个通道的大气光值，ω为透射率调整系数。
4.根据权利要求3所述的一种基于中值引导滤波的图像去雾方法，其特征在于，ω＝0.8。
5.根据权利要求1所述的一种基于中值引导滤波的图像去雾方法，其特征在于，ε＝10-6。
6.根据权利要求1所述的一种基于中值引导滤波的图像去雾方法，其特征在于，wk的大小为3×3。
7.根据权利要求1所述的一种基于中值引导滤波的图像去雾方法，其特征在于，步骤7中大气光成像模型方程为：
其中，I(x,y)为输入的含雾图像中的某一像素，A为大气光值，J(x，y)为去雾后的该像素点的值，t_d为引导滤波后的精准透射率，t0为防止t_d过小的限幅常数。
8.根据权利要求7所述的一种基于中值引导滤波的图像去雾方法，其特征在于，t0＝0.1。
CN201610307360.5A 2016-05-10 2016-05-10 一种基于中值引导滤波的图像去雾方法 Active CN105976337B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610307360.5A CN105976337B (zh) 2016-05-10 2016-05-10 一种基于中值引导滤波的图像去雾方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610307360.5A CN105976337B (zh) 2016-05-10 2016-05-10 一种基于中值引导滤波的图像去雾方法

Publications (2)

Publication Number Publication Date
CN105976337A CN105976337A (zh) 2016-09-28
CN105976337B true CN105976337B (zh) 2018-12-18

Family

ID=56992134

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610307360.5A Active CN105976337B (zh) 2016-05-10 2016-05-10 一种基于中值引导滤波的图像去雾方法

Country Status (1)

Country Link
CN (1) CN105976337B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106846259A (zh) * 2016-12-15 2017-06-13 上海交通大学 一种腹腔镜手术视频帧率自适应的实时去烟雾方法
CN109993704A (zh) * 2017-12-29 2019-07-09 展讯通信（上海）有限公司 一种去雾图像处理方法及系统
CN108230275A (zh) * 2018-02-05 2018-06-29 电子科技大学 图像去雾的方法
CN109767407B (zh) * 2019-02-27 2022-12-06 西安汇智信息科技有限公司 一种去雾过程中大气透射率图像的二次估计方法
CN109961413B (zh) * 2019-03-21 2023-05-12 平顶山学院 大气光方向优化估计的图像去雾迭代算法
CN110473155B (zh) * 2019-08-02 2023-07-18 西安工程大学 一种基于视网膜颜色感知暗通道原理的图像去雾方法
CN110660026B (zh) * 2019-08-08 2023-04-18 西安电子科技大学 一种基于Retinex理论和高饱和度先验的图像去雾方法
CN110910319B (zh) * 2019-10-30 2022-10-21 中国医学科学院北京协和医院 一种基于大气散射模型的手术视频实时去雾增强方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104281998A (zh) * 2013-07-03 2015-01-14 中山大学深圳研究院 一种基于引导滤波的快速单幅彩色图像去雾方法
CN104318535A (zh) * 2014-11-20 2015-01-28 广东欧珀移动通信有限公司 图像去雾的方法、装置及移动终端
CN104881879A (zh) * 2015-06-15 2015-09-02 北京航空航天大学 一种基于暗通道先验的遥感图像雾霾仿真方法
CN105046666A (zh) * 2015-07-24 2015-11-11 中国科学技术大学 一种基于暗原色先验的交通视频实时去雾的方法
CN105118027A (zh) * 2015-07-28 2015-12-02 北京航空航天大学 一种图像的去雾方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104281998A (zh) * 2013-07-03 2015-01-14 中山大学深圳研究院 一种基于引导滤波的快速单幅彩色图像去雾方法
CN104318535A (zh) * 2014-11-20 2015-01-28 广东欧珀移动通信有限公司 图像去雾的方法、装置及移动终端
CN104881879A (zh) * 2015-06-15 2015-09-02 北京航空航天大学 一种基于暗通道先验的遥感图像雾霾仿真方法
CN105046666A (zh) * 2015-07-24 2015-11-11 中国科学技术大学 一种基于暗原色先验的交通视频实时去雾的方法
CN105118027A (zh) * 2015-07-28 2015-12-02 北京航空航天大学 一种图像的去雾方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
"基于引导滤波的单幅图像去雾算法研究;唐鉴波 等;《科学技术与工程》;20130430;第13卷(第11期);第1-2节 *
Investigating Haze-Relevant Features in a Learning Framework for Image Dehazing;Ketan Tang et al.;《Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on》;20140628;第4.2节第1段 *
基于引导滤波的单幅图像自适应去雾算法;杨燕 等;《计算机工程》;20160131;第42卷(第1期);第265-271页 *
视频抠图算法研究;詹肇楷;《中国优秀硕士学位论文全文数据库 信息科技辑》;20160315;第1.3节 *

Also Published As

Publication number Publication date
CN105976337A (zh) 2016-09-28

Similar Documents

Publication Publication Date Title
CN105976337B (zh) 一种基于中值引导滤波的图像去雾方法
CN106296612B (zh) 一种图像质量评价和天气状况引导的分级式监控视频清晰化系统及方法
CN104253930B (zh) 一种实时视频去雾方法
CN102750674B (zh) 基于自适应容差的视频图像去雾方法
CN106157267B (zh) 一种基于暗通道先验的图像去雾透射率优化方法
CN106530257A (zh) 一种基于暗通道先验模型的遥感图像去雾方法
CN105046666B (zh) 一种基于暗原色先验的交通视频实时去雾的方法
CN106251300A (zh) 一种基于Retinex的快速夜间雾天图像复原方法
CN102063706B (zh) 一种快速去雾方法
CN108596853A (zh) 基于背景光统计模型和传输地图优化的水下图像增强方法
CN107767354A (zh) 一种基于暗原色先验的图像去雾算法
CN108230264A (zh) 一种基于ResNet神经网络的单幅图像去雾方法
CN106548463A (zh) 基于暗通道与Retinex的海雾图像自动去雾方法及系统
CN104318524A (zh) 一种基于YCbCr色彩空间的图像增强方法、装置及系统
CN111292258A (zh) 一种基于暗通道和亮通道先验的图像去雾方法
Li et al. Single image haze removal based on haze physical characteristics and adaptive sky region detection
CN104794688A (zh) 基于深度信息分离天空区域的单幅图像去雾方法及装置
CN111553862B (zh) 一种海天背景图像去雾和双目立体视觉定位方法
Wang et al. Variational single nighttime image haze removal with a gray haze-line prior
CN103996178A (zh) 一种沙尘天气彩色图像增强方法
CN104794697A (zh) 一种基于暗原色先验的图像去雾方法
CN109087254A (zh) 无人机航拍图像雾霾天空和白色区域自适应处理方法
CN110827221A (zh) 一种基于双通道先验和侧窗导向滤波的单幅图像去雾方法
CN110136079A (zh) 基于场景深度分割的图像去雾方法
CN111325688A (zh) 融合形态学聚类优化大气光的无人机图像去雾方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20210128

Address after: 210012-015, 7th floor, building 11, yunmi City, No.19, ningshuang Road, Yuhuatai District, Nanjing City, Jiangsu Province

Patentee after: Zhilian cloud big data technology Nanjing Co.,Ltd.

Address before: 710064 Shaanxi province Xi'an Beilin District South Second Ring Road No. 33 in the middle

Patentee before: CHANG'AN University