亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種并行加速改進的快速相位解包裹算法

        2021-01-08 01:37:54饒長輝高國慶周璐春
        光電工程 2020年12期
        關(guān)鍵詞:區(qū)塊像素圖像

        龍 瀟,鮑 華,饒長輝,高國慶,周璐春

        一種并行加速改進的快速相位解包裹算法

        龍 瀟1,2,3,鮑 華1,2*,饒長輝1,2,高國慶1,2,周璐春1,2

        1中國科學(xué)院自適應(yīng)光學(xué)重點實驗室,四川 成都 610209;2中國科學(xué)院光電技術(shù)研究所,四川 成都 610209;3中國科學(xué)院大學(xué),北京 100049

        針對Miguel等人提出的質(zhì)量圖引導(dǎo)相位解包裹算法中串行運算效率較低的缺點,構(gòu)造了一種多個低可靠度區(qū)塊并行合并的改進算法。在滿足原始算法設(shè)計思想的前提下,對解包裹路徑進行重新定義,并根據(jù)原始算法的解包裹路徑非連續(xù)的特性,構(gòu)建了一種低可靠度區(qū)塊亂序合并的策略,使得多個低可靠度區(qū)塊的合并任務(wù)可以同時進行。改進算法采用多線程軟件架構(gòu),主線程負責(zé)循環(huán)遍歷未處理的區(qū)塊,子線程接收待處理的區(qū)塊執(zhí)行合并任務(wù)。實驗結(jié)果表明,改進方法與原始算法的處理結(jié)果完全一致,而并行改進策略可有效利用計算機多核資源,使得相位解包裹算法的運行效率提高了50%以上。

        相位解包裹;質(zhì)量引導(dǎo);路徑相關(guān);并行計算;相位測量

        1 引 言

        在相位測量技術(shù)中,直接獲取的相位通常為被折疊進[-p,p)區(qū)間的包裹相位圖,當(dāng)相位的變化范圍超過2p時會出現(xiàn)跳變,形成條紋圖案[1]。將包裹相位恢復(fù)為連續(xù)相位的相位展開技術(shù)作為相位測量中不可缺少的關(guān)鍵技術(shù),已廣泛應(yīng)用于自適應(yīng)光學(xué)[2]、散斑干涉成像[3]、核磁共振成像[4]、光柵條紋投影測量[5]以及干涉合成孔徑雷達[6]等多個領(lǐng)域。

        近年來國內(nèi)外科研工作者針對解包裹算法做出了大量研究,已經(jīng)提出了許多相位展開算法。這些算法主要可以分為兩類[7]:1) 最小范數(shù)算法,2) 路徑跟蹤算法。

        最小范數(shù)算法是一種全局算法,將解包裹轉(zhuǎn)化為對全圖的目標函數(shù)求最小值的過程[8],當(dāng)目標函數(shù)使用1范數(shù)時,稱為最小費用流算法[9];當(dāng)目標函數(shù)使用2范數(shù)時,稱為最小二乘算法[10],對二維相位解包裹通常使用最小二乘算法。全局算法通常計算量較大,因此出現(xiàn)了對圖像分區(qū)域并行計算再對處理后的區(qū)域進行展開的方法[11],是計算速度與魯棒性的折衷。

        路徑跟蹤算法的思想是選擇合適的解包路徑,避免受噪聲影響大的區(qū)域產(chǎn)生的誤差隨著路徑不斷傳遞。路徑的選擇主要分為枝切法[6]、最小不連續(xù)算法[12]和質(zhì)量圖引導(dǎo)法[13],枝切法尋找包裹圖像中會出現(xiàn)解包錯誤的殘差點并形成枝切線,在最后單獨進行解包。枝切法計算速度較快,但枝切線出現(xiàn)封閉環(huán)時無法正確解包裹。最小不連續(xù)算法將差異大于p的相鄰像素稱為一個不連續(xù),在整個包裹相位中循環(huán)搜索不連續(xù)區(qū)域并展開,使圖像的不連續(xù)最小,循環(huán)搜索使最小不連續(xù)算法執(zhí)行效率較低。質(zhì)量引導(dǎo)路徑算法生成描述相位質(zhì)量高低的質(zhì)量圖作為引導(dǎo),優(yōu)先對高質(zhì)量像素點解包,使低質(zhì)量區(qū)域產(chǎn)生的誤差不會被傳播,因此具有較好的抗噪性能。

        質(zhì)量引導(dǎo)路徑算法中,Miguel等人[14]提出的快速相位解包裹算法高效準確且具有較強的抗噪魯棒性,但其串行的計算方式效率較低,對于如自適應(yīng)光學(xué)、干涉法測量等實時性強的任務(wù)需要進一步提速。針對上述問題,本文對Miguel等人提出的解包裹算法進行改進,提出了一種低可靠度區(qū)塊亂序并行合并的策略,改進算法與原始算法的處理結(jié)果完全一致,而解包裹算法的運行效率得到了顯著提升。

        2 原始相位解包裹算法原理

        式中:為整數(shù),對于的矩陣有,。

        質(zhì)量導(dǎo)引路徑算法的關(guān)鍵在于找到包裹相位圖中的平緩區(qū)域。將解包路徑設(shè)置在高可靠度像素之間,避免算法在噪聲較大和相位階躍處進行解包從而導(dǎo)致解包產(chǎn)生誤差并傳播至圖像其他區(qū)域。

        其中:

        定義像素的可靠度為

        則可靠度值越大表示像素越可靠。

        定義、為水平、豎直方向上兩相鄰像素間的連接(稱為邊),則邊的可靠度為

        初始化每個像素自成一組,對水平和豎直兩方向全部邊按照可靠度排序。Miguel算法的相位展開路徑從可靠度最大邊開始以降序遍歷所有邊,每條邊的解包規(guī)則為:若邊連接的兩像素屬于同一組則順序執(zhí)行下一邊;若邊連接的兩像素屬于不同組,計算它們間的2p跳變大小,像素數(shù)少的組內(nèi)所有像素點均加減該跳變值,使此邊連接的兩像素間大于p的差異被消去,此后像素數(shù)少的組整體并入像素數(shù)多的組,并繼續(xù)順序執(zhí)行下一邊。

        Miguel算法以串行的方式嚴格按照可靠度降序處理所有邊,即使下一邊與當(dāng)前邊的所屬組不包含重復(fù)的像素點,也需要將已完成的處理結(jié)果作為下一邊的處理方式的判斷依據(jù),這種設(shè)計方法使算法的執(zhí)行效率較低。

        3 多線程并行改進方法

        Miguel算法的核心思想是以最大可靠度的公共邊作為兩個相鄰組的解包路徑,當(dāng)相鄰組合并為一組后,它們原有的其他公共邊便不再執(zhí)行,避免了對低可靠度邊進行合并,從而降低解包誤差。分析Miguel算法的執(zhí)行流程可知,原始算法為串行運算,運行中同時最多只有一條邊解包,在處理高分辨率包裹相位時存在耗時較長的問題,對實時性要求較高的任務(wù),需要對此算法進行并行加速。另一方面,Miguel算法在解包過程中存在大量的邏輯分支判斷,不適用于GPU協(xié)處理器加速,因此本文選擇多核CPU硬件平臺,采用多線程任務(wù)協(xié)同的程序設(shè)計方法,使用主線程進行任務(wù)分發(fā)和邏輯判斷,子線程執(zhí)行合并任務(wù)的并行模式。

        本文改進算法的思想如圖2所示,在圖2(a)所示的4′4包裹圖像中,所有邊已按照可靠度排序為1至24。Miguel算法需要24次合并操作順序執(zhí)行,如圖2(a)所示,其中綠色邊為有效合并邊,白色邊表示相鄰像素已屬于同一組,數(shù)字代表合并時的循環(huán)次數(shù);而本文的改進策略是并行執(zhí)行低可靠度邊中互不影響的合并任務(wù),如圖2(b)所示,相同顏色的邊代表在一個循環(huán)內(nèi)完成判斷與合并,其中數(shù)字代表合并時的循環(huán)次數(shù),僅有顏色的邊表示相鄰像素已屬于同一組。改進后合并操作減少為15次,并且每次循環(huán)可同時執(zhí)行多個合并任務(wù),僅需要6次循環(huán)即可完成。

        根據(jù)Miguel算法的合并規(guī)則可知,對串行算法中的有效執(zhí)行邊亂序執(zhí)行,解包裹結(jié)果不變,僅存在因基準點不同而整體差異整數(shù)個2p,若對解包裹結(jié)果進行減均值處理則可得到相同結(jié)果。改進相位解包裹算法的關(guān)鍵在于主線程找到有效執(zhí)行邊中互不影響、可同時執(zhí)行解包的邊,并交由多子線程同時完成合并。

        圖2 改進算法的目的。(a) 順序執(zhí)行合并操作; (b)并行執(zhí)行合并操作

        圖3 并行改進算法流程圖

        圖4 鄰域模型示意圖

        Miguel算法從高可靠度到低可靠度的執(zhí)行順序使組與組之間通過最大可靠度的公共邊進行合并,若改變執(zhí)行順序,則執(zhí)行邊一定為兩相鄰組中至少一組的最大可靠度邊。根據(jù)此規(guī)律,本文提出的并行化流程圖如圖3所示,主線程循環(huán)遍歷未處理邊,檢查其所屬的兩個組,若其中一組在本循環(huán)中首次出現(xiàn),則此邊為該組的最大可靠度邊,屬于有效執(zhí)行邊,具體策略如下:對每個組增加一個“空閑/執(zhí)行/等待”標志,空閑狀態(tài)表示該組在此輪循環(huán)中可以與非執(zhí)行態(tài)的組發(fā)起合并任務(wù);執(zhí)行狀態(tài)表示該組已進入任務(wù)隊列,將在子線程內(nèi)執(zhí)行完成后置為等待態(tài);等待狀態(tài)表示該組等待空閑態(tài)的組發(fā)起合并任務(wù),主線程按可靠度從大到小循環(huán)遍歷每一條邊,并根據(jù)此標志判斷任務(wù)是否分發(fā)。

        主線程的遍歷規(guī)則為:

        1) 兩點屬于同一組,丟棄此邊;

        2) 兩個組都空閑,將這兩組的標志設(shè)為執(zhí)行,分發(fā)此邊;

        3) 組1等待,組2空閑,將這兩組的標志設(shè)為執(zhí)行,分發(fā)此邊;

        4) 組1執(zhí)行,組2空閑,將組2設(shè)為等待,跳過此邊;

        5) 組1執(zhí)行,組2執(zhí)行或等待,跳過此邊;

        6) 兩個組都等待,跳過此邊;

        7) 主線程進入下一循環(huán)時,處于“等待”態(tài)的組變?yōu)椤翱臻e”態(tài)。

        本文將組看成宏觀意義上的點,組與組之間的有效邊為公共邊中可靠度最大的一條邊,則包裹圖像可由如圖4中所示的三種鄰域模型進行描述,其中4條邊按照可靠度排序為1至4。鄰域模型中每個像素亦可表示一個鄰域模型。

        以圖4中的三種鄰域模型為例,解包順序如圖5所示。

        第一種模型,如圖4(a):第一次循環(huán),邊1執(zhí)行,邊2跳過,邊3執(zhí)行,邊4跳過;第二次循環(huán),邊2執(zhí)行,邊4跳過;第三次循環(huán)邊4成為內(nèi)部邊。執(zhí)行順序為(1,3)→(2)。

        第二種模型,如圖4(b):第一次循環(huán),邊1執(zhí)行,邊2、3、4跳過;第二次循環(huán),邊2執(zhí)行,邊3、4跳過。第三次循環(huán),邊3執(zhí)行,邊4跳過;第四次循環(huán)邊4成為內(nèi)部邊。執(zhí)行順序(1)→(2)→(3)。

        第三種模型,如圖4(c):第一次循環(huán),邊1、2執(zhí)行,邊3、4跳過;第二次循環(huán),邊3執(zhí)行,邊4跳過;第三次循環(huán)邊4成為內(nèi)部邊。執(zhí)行順序(1,2)→(3)。

        4 算法仿真

        為檢驗本文提出的低可靠度區(qū)塊亂序合并策略的有效性,實驗使用標準C++進行編程實現(xiàn)。實驗環(huán)境為Windows 10,Qt5.13.2,編譯器MSVC2017 64 bit,處理器Intel Core i5-8250U,8 GB RAM。

        實驗數(shù)據(jù)采用兩幅包裹圖像進行算法驗證。原始連續(xù)相位圖為隨機系數(shù)的前36階澤尼克組合像差和前105階澤尼克組合像差,疊加噪聲為均值為0,標準差0.5的高斯白噪聲,圖像分辨率為1024′1024,有效數(shù)據(jù)為圓域孔徑。

        圖6(a)為前36階澤尼克組合像差,經(jīng)過疊加高斯白噪聲與相位折疊得到圖6(b)的包裹圖像。圖6(c)為原算法解包裹結(jié)果,圖6(d)為改進算法解包裹結(jié)果,兩者面形完全一致。

        圖5 鄰域模型的解包流程示意圖

        圖6 對低階像差解包效果。(a) 組合像差;(b) 包裹圖像;(c) 原算法解包裹結(jié)果;(d) 改進算法解包裹結(jié)果

        圖7 對高階像差解包效果。(a) 組合像差;(b) 包裹圖像;(c) 原算法解包裹結(jié)果;(d) 改進算法解包裹結(jié)果

        圖8 解包裹過程中的分組示意圖。(a) 原算法;(b) 改進算法

        表1 兩種算法的耗時比較

        圖7(a)為前105階澤尼克組合像差,經(jīng)過疊加高斯白噪聲與相位折疊得到圖7(b)的包裹圖像。圖7(c)為原算法解包裹結(jié)果,圖7(d)為改進算法解包裹結(jié)果,兩者面形完全一致。

        本文對兩種1024′1024分辨率的組合像差進行20次實驗,統(tǒng)計原方法和改進算法的執(zhí)行時間。兩種方法的對比如表1所示。

        分析改進算法和原始算法的執(zhí)行過程,以不同顏色代表像素所屬的不同組,截取兩種方法在解包裹過程中的分組圖,可以看出原算法趨向于一個組吞并周圍的點,而改進算法適當(dāng)提前了可靠度較低的邊的執(zhí)行順序,均勻地使每個像素點與鄰近的像素合并成組,從而達到了并行加速的目的。

        實驗最后分析了在不同大小的包裹相位下,改進算法與原算法的平均耗時比較,比較結(jié)果如圖9所示。

        在不同的包裹圖像中,改進算法均達到了提速的效果,但與改進算法的理論加速存在一定差距。分析其原因有兩點:一是在算法后期,可同時執(zhí)行的邊太少,算法退化為串行算法;二是前期小區(qū)塊間的合并造成后續(xù)合并任務(wù)的運算量增大,降低了計算效率。通過大量仿真實驗,結(jié)合本文實驗平臺CPU調(diào)度開銷,本文提出的算法的最優(yōu)線程數(shù)為4。

        5 結(jié) 論

        本文針對Miguel等人提出的相位解包裹算法,提出了一種基于低可靠度區(qū)塊亂序合并的改進策略,使之可以并行計算。改進算法與原算法進行對比,仿真實驗表明,改進算法達到了降低原算法耗時并且不影響解包裹結(jié)果的目的,證實了改進策略正確有效。目前改進算法的加速比與計算機的核心數(shù)存在一定差距,是否可以進一步提高并行化效率或利用空余計算能力提高算法精度,是下一步研究工作的重點。

        [1] Fornaro G, Franceschetti G, Lanari R,. Robust phase-unwrapping techniques: a comparison[J]., 1996, 13(12): 2355–2366.

        [2] Bao H, Rao C H, Tian Y,. Research progress on adaptive optical image post reconstruction[J]., 2018, 45(3): 170730. 鮑華, 饒長輝, 田雨, 等. 自適應(yīng)光學(xué)圖像事后重建技術(shù)研究進展[J]. 光電工程, 2018, 45(3): 170730.

        [3] Wu S J, Yang J, Pan S Y,. Dynamic deformation measurement of discontinuous surfaces using digital speckle pattern interferometry and spatiotemporal three-dimensional phase unwrapping[J]., 2018, 47(2): 0212002. 吳思進, 楊靖, 潘思陽, 等. 數(shù)字散斑干涉術(shù)和時空三維相位解包裹用于非連續(xù)表面動態(tài)變形測量[J]. 光子學(xué)報, 2018, 47(2): 0212002.

        [4] Chavez S, Xiang Q S, An L. Understanding phase maps in MRI: a new cutline phase unwrapping method[J]., 2002, 21(8): 966–977.

        [5] Yuan H X, Li J L, Luo P,. Image restoration for blurred fringes of rail profile 3D online measurement based on PMP[J]., 2017, 44(7): 695?700. 袁宏翔, 李金龍, 羅鵬, 等. 基于PMP的鋼軌三維形貌在線測量模糊條紋復(fù)原[J]. 光電工程, 2017, 44(7): 695–700.

        [6] Goldstein R M, Zebker H A, Werner C L. Satellite radar interferometry: two-dimensional phase unwrapping[J]., 1988, 23(4): 713–720.

        [7] Wang Y H, Chen W J, Zhong S M,. Research progress in phase unwrapping technology and its applications[J]., 2018, 37(12): 1–7, 16. 王永紅, 陳維杰, 鐘詩民, 等. 相位解包裹技術(shù)及應(yīng)用研究進展[J]. 測控技術(shù), 2018, 37(12): 1–7, 16.

        [8] Qian XF, Zhang Y A, Li X Y,. Phase unwrapping algorithm based on mask and least-squares iteration[J]., 2010, 30(2): 440–444. 錢曉凡, 張永安, 李新宇, 等. 基于掩膜和最小二乘迭代的相位解包裹方法[J]. 光學(xué)學(xué)報, 2010, 30(2): 440–444.

        [9] Costantini M. A novel phase unwrapping method based on network programming[J]., 1998, 36(3): 813–821.

        [10] Guo Y, Chen X T, Zhang T. Robust phase unwrapping algorithm based on least squares[J]., 2014, 63: 25–29.

        [11] Huang K M, Yamada T. Phase unwrapping by regions using least-squares approach[J]., 1998, 37(11): 2965–2970.

        [12] Flynn T J. Two-dimensional phase unwrapping with minimum weighted discontinuity[J]., 1997, 14(10): 2692–2701.

        [13] Xu W, Cumming I. A region-growing algorithm for InSAR phase unwrapping[J]., 1999, 37(1): 124–134.

        [14] Herráez M A, Burton D R, Lalor M J,. Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path[J]., 2002, 41(35): 7437–7444.

        Improved fast phase unwrapping algorithm based on parallel acceleration

        Long Xiao1,2,3, Bao Hua1,2*, Rao Changhui1,2, Gao Guoqing1,2, Zhou Luchun1,2

        1Key Laboratory of Adaptive Optics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;2Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;3University of Chinese Academy of Sciences, Beijing 100049, China

        (a) Group diagram of original algorithm; (b) Group diagram of parallel algorithm

        Overview:In the phase measurement technology, the phase directly obtained is usually folded into a range of wavelengths, so that the phase pattern appears as a stripe pattern. Generally, this stripe pattern is not the final result required for phase measurement. There is a phase unwrapping operation needed to obtain a continuous phase map.

        The main problem facing by the phase unwrapping algorithms is how to balance the robustness and the computational efficiency. At present, there have been a lot of researches on unwrapping algorithms. They are mainly divided into two categories, including the minimum norm method and the path tracking method. The minimum norm algorithm is a global algorithm that transforms the process of phase expansion into a process of minimizing an objective function of the full graph. In the least norm method, the least squares algorithm is a commonly used unwrapping algorithm for two-dimensional phase wrapping images. Because global algorithms usually need a large amount of calculations, they require high computing power. The idea of the path tracking algorithm is to choose a suitable path for expansion, so as to avoid that the areas affected by noise appear prematurely in the path and cause errors and continue to be transmitted along the path. Quality map guidance algorithm is a common type of path tracking algorithm. This algorithm first generates a quality map describing the impact of noise. The quality map guides the unwrapping path through high-quality areas, so that the errors generated in low-quality areas will not be propagated. Thus, quality map guidance algorithm has good noise immunity.

        The quality map-guided unwrapping algorithm proposed by Miguel needs a small amount of calculations and has strong noise immunity, but its serial calculation method has low operating efficiency. To solve this problem, an improved algorithm for parallel merging of multiple low-reliability blocks is proposed. Under the condition that the original algorithm design idea is satisfied, the unwrapping path is redefined as the largest reliable edge of the block. In addition, based on the non-continuous characteristic of the unwrapping path of the original algorithm, a low-reliability block out-of-order merging strategy is proposed to make multiple merging tasks can be performed simultaneously. The improved algorithm uses a multi-threaded software architecture. The main thread is responsible for looping through the unprocessed blocks to check whether they meet the requirements of merging, and the child threads receive and perform the merge tasks. The experimental results show that the improved method is completely consistent with the processing results of the original algorithm, and the parallel improvement strategy can effectively use the computer's multi-core resources, so that the operational efficiency of the phase unwrapping algorithm is improved by more than 50%.

        Citation: Long X, Bao H, Rao C H,. Improved fast phase unwrapping algorithm based on parallel acceleration[J]., 2020,47(12): 200111

        * E-mail: hbao@ioe.ac.cn

        Improved fast phase unwrapping algorithm based on parallel acceleration

        Long Xiao1,2,3, Bao Hua1,2*, Rao Changhui1,2, Gao Guoqing1,2, Zhou Luchun1,2

        1Key Laboratory of Adaptive Optics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;2Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, Sichuan 610209, China;3University of Chinese Academy of Sciences, Beijing 100049, China

        Aiming at the shortcoming of low serial operational efficiency in the quality-map-guided phase-unwrapping algorithm proposed by Miguel, an improved algorithm for parallel merging of multiple low-reliability blocks is proposed. Under the condition that the original algorithm design idea is satisfied, the unwrapping path is redefined as the largest reliable edge of the block. In addition, based on the non-continuous characteristic of the unwrapping path of the original algorithm, a low-reliability block out-of-order merging strategy is proposed to make multiple merging tasks can be performed simultaneously. The improved algorithm uses a multi-threaded software architecture. The main thread is responsible for looping through the unprocessed blocks to check whether they meet the requirements of merging, and the child threads receive and perform the merge tasks. The experimental results show that the improved method is completely consistent with the processing results of the original algorithm, and the parallel improvement strategy can effectively use the computer's multi-core resources, so that the operational efficiency of the phase unwrapping algorithm is improved by more than 50%.

        phase unwrapping; quality guidance; path dependent; parallel computing; phase measurement

        National Natural Science Foundation of China (11727805)

        10.12086/oee.2020.200111

        TP391;TN29

        A

        龍瀟,鮑華,饒長輝,等. 一種并行加速改進的快速相位解包裹算法[J]. 光電工程,2020,47(12): 200111

        : Long X, Bao H, Rao C H,Improved fast phase unwrapping algorithm based on parallel acceleration[J]., 2020, 47(12): 200111

        2020-04-02;

        2020-05-27

        國家自然科學(xué)基金資助項目(11727805)

        龍瀟(1994-),男,碩士,主要從事相位差波前測量技術(shù)的研究。E-mail:836868760@qq.com

        鮑華(1981-),男,博士,副研究員,主要研究自適應(yīng)光學(xué)圖像事后處理及相位差法波前探測。E-mail:hbao@ioe.ac.cn

        猜你喜歡
        區(qū)塊像素圖像
        趙運哲作品
        藝術(shù)家(2023年8期)2023-11-02 02:05:28
        改進的LapSRN遙感圖像超分辨重建
        像素前線之“幻影”2000
        區(qū)塊鏈:一個改變未來的幽靈
        科學(xué)(2020年5期)2020-11-26 08:19:12
        有趣的圖像詩
        區(qū)塊鏈:主要角色和衍生應(yīng)用
        科學(xué)(2020年6期)2020-02-06 08:59:56
        “像素”仙人掌
        區(qū)塊鏈+媒體業(yè)的N種可能
        傳媒評論(2018年4期)2018-06-27 08:20:12
        讀懂區(qū)塊鏈
        高像素不是全部
        CHIP新電腦(2016年3期)2016-03-10 14:22:03
        日本熟妇美熟bbw| 亚洲最大在线精品| 亚洲日韩一区二区一无码| 亚洲欧美成人久久综合中文网| 亚洲成人av大片在线观看| 女同同志熟女人妻二区| 青青草原综合久久大伊人精品| 亚洲一区 日韩精品 中文字幕| 久久久久国产一区二区三区| 国产一区二区三区四区五区vm| 手机看片1024精品国产| 亚洲伦理一区二区三区| 国产日产免费在线视频| 日本国产一区二区在线| 人妖av手机在线观看| 丰满人妻一区二区三区免费视频| 日本japanese丰满多毛| 在线视频制服丝袜中文字幕| 黑人一区二区三区啪啪网站| 亚洲av成人av三上悠亚| 国产欧美性成人精品午夜| 精品人人妻人人澡人人爽牛牛| 国产成人无码精品久久99| 亚洲美女av二区在线观看| 国产91极品身材白皙| 无码gogo大胆啪啪艺术| 国产精品美女久久久浪潮av| 2021av在线| 亚洲精品视频一区二区三区四区| 国产日产精品_国产精品毛片| 亚洲精品国产美女久久久| 无码手机线免费观看| 国产久视频| 少妇激情一区二区三区| 国产精品亚洲一级av第二区| 人妻丰满熟妇av无码区app| 午夜性无码专区| 国产精品偷伦免费观看的| 日本女u久久精品视频| 久久久久国色av免费观看性色| 1000部精品久久久久久久久|