久草青青在线,亚洲久久日婷婷AV

在我的上一篇博客Unity custom shader中調(diào)用內(nèi)置Lightmap和Light Probes中有提到LightProbes以及如何在custom shader中獲取LightProbes中所存儲的數(shù)據(jù)。然而，LightProbes中的數(shù)據(jù)是如何來的我卻一無所知，于是我去谷歌了一下，不查不知道，一查嚇一跳，這背后有著復(fù)雜的理論依據(jù)（包括圖形學(xué)、概率論、信號分析、微積分等），以及令人望而生畏的數(shù)學(xué)計算，才最終實現(xiàn)了官方文檔所說的：

Light Probes provide a way to capture and use information about light that is passing through the empty space in your scene

在此，本篇博客會依據(jù)論文Spherical Harmonic Lighting: The Gritty Details盡可能多的介紹球諧光照，當(dāng)然，有很多地方我也不是很明白，如果有說錯的地方望大家指正。

首先，球諧光照在現(xiàn)代游戲圖形渲染領(lǐng)域有著廣泛的應(yīng)用，可快速模擬復(fù)雜的實時光照，unity中的LightProbes就用了此技術(shù)，當(dāng)然unity中那些不重要的實時光源也用了。所以對于LightProbes來說，其背后就是球諧光照技術(shù)的應(yīng)用。

那么，什么是球諧光照呢？球諧光照（Spherical Harmonic Lighting）從英文字面上來看，是一種光照算法，這種算法的內(nèi)核是定義在球面上的特殊函數(shù)，它可以對光照進行捕捉（capture）并在之后進行重新光照（relight）并且可以實時展示全局光照（global illumination）風(fēng)格的區(qū)域光源（area light）與軟陰影（soft shadow）。

既然是一種光照算法，那么論文一開始就是從一個光照公式開始談起的:
$L(X, \vec\omega_o) = L_e(X, \vec\omega_o) + \int_sf_r(X, \vec\omega_i\rightarrow\vec\omega_o)L(X', \vec\omega_i)G(X,X')V(X,X'){\rm d}\omega_i$
其中：
$L(X, \vec\omega_o)$ = $\omega_o$ 方向的光照在點 $X$ 處所反射出來的光的強度
$L_e(X, \vec\omega_o)$ = 點 $X$ 處的自發(fā)光
$f_r(X, \vec\omega_i\rightarrow\vec\omega_o)$ = 點 $X$ 處的BRDF(Bidirectional Reflectance Distribution Function)
$L(X', \vec\omega_i)$ = 從另一個物體上的點 $X’$ 傳遞過來對點 $X$ 有影響的光，并且這道光的入射角度為 $\omega_i$
$G(X,X')$ = $X$ 與 $X’$ 之間的幾何關(guān)系，一般就是 $\vec N \cdot \vec L$ ，也就是 $cosine \ term$
$V(X,X')$ = 可見性測試， $X$ 能看見 $X’$ 就返回1，否則返回0

對于 $f_r(X, \vec\omega_i\rightarrow\vec\omega_o)$ 、 $L(X', \vec\omega_i)$ 、 $G(X,X')$ 、 $V(X,X')$ 這四項進行實時積分是很困難而且昂貴的，所以我們需要一些數(shù)學(xué)工具來簡化這一積分計算，變不可能為可能！

蒙特卡洛積分
首先我們來看第一個數(shù)學(xué)工具，蒙特卡洛積分。公式是這樣子的：
$\int f(x){\rm d}x = \int \frac {f(x)}{p(x)} p(x) {\rm d}x \approx \frac {1}{N} \sum^N_{i=1}\frac {f(x)}{p(x)}$

這個東西的核心思想就是對于一個連續(xù)函數(shù)f(x)在其定義域隨機挑選N個值 $x_i$ （這個操作用個詞概括就叫采樣）并求出N個 $f(x_i)$ ，除以每個值被挑選中的概率 $p(x_i)$ ，將他們累加起來除以N，就近似得到了原函數(shù) $f(x)$ 的積分結(jié)果。那么對于這個式子而言，采樣越多則結(jié)果越逼近真實的積分結(jié)果。
不過由于這個積分是個球面積分，定義域是在球面上，要在定義域上進行均勻的采樣我們還需要一個公式：
$(2arccos(\sqrt{1-\xi_x}),2\pi\xi_y) \rightarrow (\theta,\psi)$

其中 $\xi_x$ 和 $\xi_y$ 是0-1之間的均勻隨機變量，他們倆是獨立的。 $\theta,\psi$ 是采樣方向，在球面上可以理解為維度經(jīng)度。

用上面這個公式就可以均勻隨機的在球面上進行采樣，

我們知道單位球面積為 ${4\pi}$ ，那么隨機采樣一個點的概率為 $\frac{1}{4\pi}$ ，

再利用蒙特卡洛積分來近似計算結(jié)果。如果有朋友對這個公式怎么來的有興趣，可以看看Generating uniformly distributed numbers on a sphere這篇博客。

基底函數(shù)
信號處理中有個叫傅里葉變換的東西（我沒學(xué)過信號處理，沒法專業(yè)的解釋這東西，需要專業(yè)解釋還是請教谷歌吧），簡單來說就是對于任意一個函數(shù) $f(x)$ ，可以把它拆成一組三角函數(shù)乘以系數(shù)之后的和，即 $F(x) = \sum c * f(x)$ 。這里的 $f(x)$ 為一組正交基底函數(shù)，函數(shù)正交意味著有兩個函數(shù) $f(x)$ 和 $g(x)$ $\int f(x)g(x){\rm d}x = 0$ 。
我們有了公式 $F(x) = \sum c * f(x)$ 后，就要開始找出 $c$ 和 $f(x)$ 。 $f(x)$ 是一組正交基底函數(shù)，古往今來已經(jīng)有很多前輩探索發(fā)現(xiàn)了許多的正交基底函數(shù)，我們在這里特別關(guān)注的是一個叫做伴隨勒讓德多項式的東西，不過首先我們先把這個放一邊，認為 $f(x)$ 是已知的，那么現(xiàn)在唯一要求的是 $c$ 。 $c$ 的求法是在 $F(x)$ 上采樣一堆 $x$ 繪制出 $F(x)$ 的曲線，乘以 $f(x)$ 再積分（積分相當(dāng)于把 $F(x)$ 投影到 $f(x)$ 上來求 $f(x)$ 對 $F(x)$ 的貢獻度），得到的值就是我們要求的 $c$ 。這一過程在原文中被稱為投影（projection），所以本篇博客中所有出現(xiàn)投影一詞時都是這個意思?，F(xiàn)在 $c$ 有了 $f(x)$ 有了，我們就可以大致求出原函數(shù) $F(x)$ ，過程如下圖。

只要我們采樣足夠多，選定的基底函數(shù)足夠多，計算出來的函數(shù)就越接近原函數(shù)。另外基底函數(shù)會通過權(quán)重進行排序，一般排在前面的基底函數(shù)比較重要，也稱為低頻部分，它為定下了基調(diào)，通常第一個基底函數(shù)模擬的是一個平均值，越往后面的基底函數(shù)重要性越低，被稱為高頻部分，展現(xiàn)了的一些噪聲部分。

現(xiàn)在可以來看 $f(x)$ 了。之前我們看到了傅里葉變換可以把任意一個函數(shù) $F(x)$ 拆分成一組正交基底函數(shù)， $F(x)$ 是基于x軸的。而我們之前看到的渲染函數(shù)卻是基于球面的，基于球面的函數(shù)一般表達為 $F(\theta,\psi)$ ，兩個參數(shù)是單位球在兩個軸的夾角。在球面上也有類似于傅里葉變換的東西，叫勒讓德多項式，特別的，我們會關(guān)注一個叫伴隨勒讓德多項式的東西。根據(jù)Wikipedia，伴隨勒讓德多項式與勒讓德多項式的關(guān)系是：

伴隨勒讓德多項式可以由勒讓德多項式求 m 次導(dǎo)得到：
$P^m_l(x) = (1-x^2)^{m/2}P^{(m)}_l(x)$
等號右邊的上標 (m) 表示求 m 次導(dǎo)。

公式中 $l$ 為band index， $m$ 為degree，什么意思呢？伴隨勒讓德多項式分了很多個band，每個band里面又用 $m$ 作為index標注多項式，像這樣：
$P^0_0(x)$
$P^0_1(x)P^1_1(x)$
$P^0_2(x)P^1_2(x)P^2_2(x)$
$\ldots\ldots\ldots\ldots$
其中要注意 $m \in [0,l], l \in R^+$ ，定義域為[-1,1]以及只返回實數(shù)（勒讓德多項式返回的是復(fù)數(shù)）
一般我們用下面三條規(guī)則來推導(dǎo)出所有伴隨勒讓德多項式，而不是直接一個個算。
1. $(l-m)P^m_l = x(2l-l)P^m_{l-1} - (l+m-1)P^m_{l-2}$
2. $P^m_m = (-1)^m(2m-1)!!(1-x^2)^{m/2}$
3. $P^m_{m+1} = x(2m+1)P^m_m$
其中2可以算出第一項，3可以推出下一項，然后用1不斷往下推導(dǎo)。
下面獻上知乎專欄雞哥的靈魂p圖來幫助大家理解這三條遞推規(guī)則。

另外，!!是雙階乘的意思，這是定義，不要自己亂來！
到此，我以為這就是whole story了。。?？上?，我還是太年輕了=_=!
事實上，以上所講的伴隨勒讓德多項式依舊是一維的情況，要把定義在球面上的一個函數(shù)拆分成一組正交基底函數(shù)（我們把這組函數(shù)稱之為球諧函數(shù)，下文都會用球諧函數(shù)表達），伴隨勒讓德多項式是核心，真正的公式是這樣的：

$y^m_l(\theta,\psi) = \begin{cases} \sqrt{2}K^m_lcos(m\psi)P^m_l(cos\theta), & m>0 \\ \sqrt{2}K^m_lsin(-m\psi)P^{-m}_l(cos\theta), & m<0 \\ K^0_lP^0_l(cos\theta), & m=0 \\ \end{cases}$

其中 $P$ 就代表伴隨勒讓德多項式，這里多出來的 $K$ 是一個縮放系數(shù)
$K^m_l = \sqrt{\frac{2l+1}{4\pi} \cdot \frac{(l-|m|)!}{(l+|m|)!}}$

而 $m$ 和伴隨勒讓德多項式中的定義有所不同，這里m可以取負數(shù)，即： $l \in R^+， -l\leq m\leq l$

有時候把球諧函數(shù)變成一維形式也挺有用的，所以我們定義 $y_i$
$y^m_l(\theta,\psi) = y_i(\theta,\psi)$ 其中 $i=l(l+1)+m$

原文中說大部分的論文都是只扔給讀者一堆令人疑惑的多項式，但把公式轉(zhuǎn)成圖會更有趣，所以原文作者給了一張前5個band的可視化圖:

是不是覺得這么一堆公式看著頭都大了，遑論寫代碼了。Wikipedia
貼心的把前幾個band的公式都解析出來了，照抄寫死在代碼都行！

還記得之前的公式 $F(x) = \sum c * f(x)$ 以及 $c$ 求法么？（不記得回上去看下） $f(x)$ 就是上面所說的球諧函數(shù)，原函數(shù) $F(x)$ 通過隨機在定義域上采樣再通過蒙特卡洛積分得到，那么 $c$ 就可以得到這樣一個公式
$c^m_l=\int_sf(s)y^m_l(s){\rm d}s$
那么近似原函數(shù) $\tilde{f(x)}$ 即是：
$\tilde{f(x)} = \sum^{n-1}_{l=0}\sum^l_{m=-l}c^m_ly^m_l(s)=\sum^{n^2}_{i=0}c_i y_i(s)$
以上公式中 $s$ 可以理解為采樣點。那為什么最后是從 $0$ 累加到 $n^2$ 個參數(shù)呢？由于m取 $[-l,l]$ ，那么伴隨勒讓德多項式就變成了
$P^0_0(x)$
$P^{-1}_1(x)P^0_1(x)P^1_1(x)$
$P^{-1}_2(x)P^{-2}_2(x)P^0_2(x)P^1_2(x)P^2_2(x)$
$\ldots\ldots\ldots\ldots$
那么取一個band時（ $l=0$ ）有1個多項式（即一個參數(shù) $c$ ），取兩個band時（ $l=1$ ）有4個多項式（即4個參數(shù) $c$ ），取三個band時（ $l=2$ ）有9個多項式（即9個參數(shù) $c$ ），以此類推，取n個band時則有 $n^2$ 個參數(shù) $c$ 。
同樣，原文作者也提供了一張可視化的圖供大家加深印象

可以看出，n越大越能還原，即球諧函數(shù)band展開的越多，越能還原最初的信號。

上面所說的蒙特卡洛積分以及基底函數(shù)這些數(shù)學(xué)的東西其實已經(jīng)可以還原原來的光信號了，然而。。。一看論文怎么才到一半左右，后面還有一半是什么鬼?。吭瓉?，接下來論文中闡述了一些球諧函數(shù)的特性，基于這些特性，我們認為球諧函數(shù)是優(yōu)秀的！

1.標準正交性
這一特性表明對于球諧函數(shù)，使其變成一維 $y_i$ （上文有說如何變成一維）則有：
$\begin{cases} \int y_iy_j=1 & i=j \\ \int y_iy_j = 0 & i\not =j \\ \end{cases}$

2.旋轉(zhuǎn)不變性
我們定義在球面上的原函數(shù) $f(s)$ ，將其旋轉(zhuǎn)，變成另一個函數(shù) $g(s)$ ，那么 $g(s) = R(f(s))$ ，也就是說旋轉(zhuǎn)函數(shù)本身的操作與將自變量旋轉(zhuǎn)后交給 $f(s)$ 是一樣的，這個特性有什么用呢？原文的答案是我們在做動畫、移動光源、或旋轉(zhuǎn)物體的時候，光強不會產(chǎn)生漲落（fluctuate）、挪移（crawl）、抖動（pulse）等不良現(xiàn)象（objectionable artifacts）。

3.殺手級特性
不要驚訝，原文對于這個特性的表述就是the killer one。我們在做光照計算的時候，要考慮入射光的強度乘以表面反射項（這個東西稱為傳輸函數(shù)（transfer function）），以此來獲得最終的反射光。而這些都是在單位圓的表面做的，所以會有以下公式：
$\int _sL(s)t(s){\rm d}s$
其中 $L$ 為入射光而 $t$ 為傳輸函數(shù)
如果我們把這兩個函數(shù)做投影得到一堆系數(shù) $c$ 時，那么球諧函數(shù)的標準正交性可以保證有：
$\int_s\tilde L(s)\tilde t(s){\rm d}s = \sum^{n^2}_{i=0}L_it_i$
即兩個函數(shù)乘積的積分等于他們兩個球諧系數(shù)的乘積的和
于是，我們一開始無法實時計算的 $\int_sf_r(X, \vec\omega_i\rightarrow\vec\omega_o)L(X', \vec\omega_i)G(X,X')V(X,X'){\rm d}\omega_i$ 這個東西就可以快速計算了。公式中 $L(X', \vec\omega_i)$ 這項是在描述入射光函數(shù)即 $\tilde L(s)$ ，剩下的則是傳輸函數(shù) $\tilde t(s)$ ，分別用蒙特卡洛積分法近似得出原函數(shù)表達式后做投影，得出一堆 $c$ 后兩兩相乘再求和，結(jié)果就近似是原來那個難以積分的式子的結(jié)果。考慮到實際運行中 $n$ （代表多少基底函數(shù)系數(shù)）一般都取3-6之間的數(shù)（也就是說 $n^2$ 最多是36，36對基底函數(shù)系數(shù)相乘再相加），計算量并不大，所以可以快速得出結(jié)果。

旋轉(zhuǎn)球諧函數(shù)（Rotating Spherical Harmonics）
原文在這一小節(jié)里闡述了一下如何旋轉(zhuǎn)那一堆我們得到的 $n^2$ 個參數(shù)，畢竟一旦物體旋轉(zhuǎn)了（比如做動畫），我們不可能去重新算一遍參數(shù)。原文這里看得我想死（數(shù)學(xué)太差T_T），大概就是構(gòu)造了一個 $n^2\times n^2$ 的旋轉(zhuǎn)矩陣來把那 $n^2$ 個參數(shù)進行旋轉(zhuǎn)，不用重新算。至于詳細。。?？丛幕蛘哌@里球諧光照與PRT學(xué)習(xí)筆記（四）：球諧函數(shù)的性質(zhì)與球諧旋轉(zhuǎn)。

接下來，原文中討論了好幾種傳輸函數(shù)，我這邊簡單概括下。
1.無陰影漫反射傳輸（Diffuse Unshadowed Transfer）

回到最初給出的公式，把不必要的東西都剝離，只剩光源和某平面上受光的點，不考慮陰影，則：
$L(x,\omega_o) = \int_sf_r(x,\omega_o,\omega_i)L_i(x,\omega_i)H(x,\omega_i){\rm d}\omega_i$
其中：
$L(x,\omega_o)$ = 點 $x$ 處沿 $\omega_o$ 方向的出射光 $\int_sf_r(x,\omega_o,\omega_i)$ = 點 $x$ 處的BRDF
$L_i(x,\omega_i)$ = 點 $x$ 處沿 $\omega_i$ 方向的入射光
$H(x,\omega_i)$ = $cosine \ term$

由于漫反射的BRDF在每個方向都反射相同的光，即 $view \ independent$ (與觀察方向無關(guān))，那么公式中的 $\omega_o$ 可以忽略，它就變成了一個常量（原文稱為 $constant \ scalar$ ），就可以提到積分外，公式就變?yōu)椋?br> $L_{DU}(x) = \frac{\rho_x}{\pi}\int_\Omega L_i(x,\omega_i)max(N_x \cdot \omega_i){\rm d}\omega_i$
其中：
$\rho_x =$ 點 $x$ 處的表面反照率（ $surface \ albedo$ ）
$N_x =$ 點 $x$ 處的法線

公式中的這個 $\pi$ 是怎么來的呢？我只想說這個不是瞎湊湊出來的，而是由數(shù)學(xué)推導(dǎo)出來的，詳見如何看懂這些"該死的"圖形學(xué)公式。

這時，傳輸函數(shù)就變?yōu)楹唵蔚?img class="math-inline" src="https://math.jianshu.com/math?formula=cosine%20%5C%20term" alt="cosine \ term" mathimg="1">，我們記為 $M_{DU}$
$M_{DU} = max(N_x \cdot \omega_i)$

Diffuse Unshadowed Transfer效果圖

2.有陰影漫反射傳輸（Shadowed Diffuse Transfer）

光有可能被擋到的情況

我們在 $L_{DU}(x) = \frac{\rho_x}{\pi}\int_\Omega L_i(x,\omega_i)max(N_x \cdot \omega_i){\rm d}\omega_i$ 這個公式基礎(chǔ)上添加一個可見性項 $V(\omega_i)$ 來獲得自陰影，則
$L_{DU}(x) = \frac{\rho_x}{\pi}\int_\Omega L_i(x,\omega_i)V(\omega_i)max(N_x \cdot \omega_i){\rm d}\omega_i$
其中來自 $\omega_i$ 方向可見返回1，否則返回0
那么:
$M_{SD} = V(\omega_i)max(N\cdot \omega_i,0)$

Shadowed Diffuse Transfer效果圖

3.有相互反射的漫反射傳輸（Diffuse Interreflected Transfer）
這種方法最具有挑戰(zhàn)性。方程最有意思的部分是它不僅將光源的直射光部分加進來，還將從其他可見點反射過來的光也考慮進來，原文中用了一個詞是recursively，我的理解是由于光的反射不止一次，會考慮多次反射直至光被完全吸收或者超出視線范圍，這么多次反射的光通過遞歸計算得到。于是：
$L_{DI}(x) = L_{DS}(x) + \frac{\rho_x}{\pi}\int_\Omega \overline L(x',\omega_i)(1-V(\omega_i))max(N\cdot \omega_i){\rm d}\omega_i$
其中：
$L_{DS}(x) =$ 有陰影漫反射傳輸（Shadowed Diffuse Transfer）所述公式
$V(\omega_i) =$ 可見性測試
$\overline L(x',\omega_i) =$ 同一模型下由點 $x'$ 反射向點 $x$ 的光，具體方向為 $\omega_i$
從數(shù)學(xué)上來講，這個公式很難簡潔的表述，但從算法角度卻比較容易表述，總共分四步走：

對于模型上的每個著色點 $x$ ，計算這個點的直接光傳輸函數(shù)（例如：用 有陰影漫反射傳輸（Shadowed Diffuse Transfer）的公式計算）。
然后從現(xiàn)在的點發(fā)一條光線直到它遇到物體上面的另一個三角形，獲得了一個撞擊點，那么在三角形重心坐標體系下，對三角形每個頂點上的傳輸函數(shù)進行線性插值求出撞擊點的傳輸函數(shù)，這個傳輸函數(shù)就包含了返射回 $x$ 點的光的信息。
將上一步獲得的傳輸函數(shù)乘以 $cosine \ term$ 然后累加起來，再用蒙特卡洛積分法算出點 $x$ 的傳輸函數(shù)參數(shù)。
一旦所有的點都算好了，這組新的參數(shù)就是一次反射所對應(yīng)的參數(shù)，要獲得N次反射所對應(yīng)的參數(shù)，就用第一組參數(shù)作為開始設(shè)置不斷做光線追蹤，做到N次上限后停止。最后把直接光部分加進來就大功告成了。

示意圖

幾何上來講這個點子很簡單（呵呵。。。簡單尼瑪呢），模型上的每個點的直接光信息都是已知并以傳輸函數(shù)的形式進行編碼，然后我們發(fā)射處=出一些射線來找到一些能反射過來光的點，把這些點的傳輸函數(shù)乘以 $cosine \ term$ 再加到原來發(fā)出點的傳輸函數(shù)中，舉個例子，上圖點A發(fā)出一條射線，撞擊點為B，那么B點的傳輸函數(shù)加到A是這樣子的：

注意所有球諧函數(shù)都在同一坐標體系下，所以相加是種可行的操作。這樣點A在被球諧光源所照射時，即時不能直接看見，也能被照亮。這個假設(shè)很重要的一點是，光照在整個模型上沒有變化（例如點B處要和點A處有同樣的光照函數(shù)）。這也是球諧光照的關(guān)鍵：低頻光源以及物體上光源變化很小。
最終實現(xiàn)出來的效果如下：

最后，球諧漫反射表面的渲染通過那個殺手級特性計算，雖然上面說過了，不過原文還是把公式再列了一遍。
$\int_s\tilde L(s)\tilde t(s){\rm d}s = \sum^{n^2}_{i=0}L_it_i = L_0t_0 + L_1t_1+L_2t_2+L_3t_3\ldots$

最最后論文寫了一些照明與顏色標準的東西，這個由CIE（International Commission on Illumination，為什么縮寫是CIE呢？原來這個組織的法語是 Commission internationale de l'éclairage，CIE是該組織法語的縮寫=_=）組織定制，這些東西和球諧函數(shù)有關(guān)系，但更和Real Time Rendering那本書里的第八章有關(guān)系（我指的是第四版，第三版不是第八章，具體第幾章我也不知道。。。），所以還是擱置到以后再說吧。

PS. 這個月就寫了這么一篇博客，一則因為公司事多，二則因為這個月公司去了三亞年會旅游，搞掉一周，所以心思不在寫博客上，希望下個月可以高產(chǎn)一點，多學(xué)點！

參考
簡單理解spherical harmonic lighting（球諧光照）
球諧光照與PRT學(xué)習(xí)筆記（一）：引入
 球諧光照與PRT學(xué)習(xí)筆記（二）：蒙特卡洛積分與球面上的均勻采樣
 球諧光照與PRT學(xué)習(xí)筆記（三）：球諧函數(shù)
球諧光照與PRT學(xué)習(xí)筆記（四）：球諧函數(shù)的性質(zhì)與球諧旋轉(zhuǎn)
球諧光照與PRT學(xué)習(xí)筆記（五）：預(yù)計算傳輸與著色

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

球諧光照入門簡介

球諧光照入門簡介

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

球諧光照入門簡介

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av