zoukankan      html  css  js  c++  java
  • Computer Graphics note(4):Shading

    目录

    一.Shading Model(Blinn-Phong Reflectance Model)

    Games101 lecture7-8-9-10
    Shading(着色)定义为对不同对象应用不同材质的过程。不同的材质也就是不同的着色方法。有许多着色模型,例如Blinn-Phong Reflectance Model(经验模型)。如下图着色模型:

    从上面来看,分为三部分,高光(Specular highlights)漫反射(Diffuse reflection)环境光照(Ambient lighting)

    1.相关定义(Shading is Local)

    对于着色,现阶段光照是正对某一点的,该点称为Shading Point(着色点),虽然点可以位于任何一个曲面上,但在一个足够小的局部范围下,总是认为该点处于一个平面上(如上图)。从而可以定义该平面法线(pmb{n})(Surface normal),观测方向(pmb{v})(Viewer direction)(Shading Point看向观测点方,观测点减去着色点然后归一化),以及光照方向(pmb{l})(Light direciton,光源位置减去点的位置然后归一化即可)。以上向量表示方向,故为单位向量。除此之外还有着色点的颜色等参数。

    同时着色具有局部性(Local),即考虑着色时只考虑Shading Point本身,但是不考虑其他物体的存在(比如该点是否在阴影内,被遮挡情况)

    2.Diffuse Reflection(Blinn-Phong)

    所谓漫反射,也就是光线在反射的时候时往在各个方向上均匀散射的。

    (1)Shading Point接收的能量强度

    当Shading Point表面法线和光线有一定夹角( heta)的时候,夹角越小,亮度越大,如下图所示。因为角度不同时Shading Point周围的单位面积(固定)所接受到的能量情况是不同的,即有(cos heta=pmb{l}·pmb{n})

    (2)到达Shading Point的能量强度

    这里的光源假设认为是点光源(三维空间),它会往四周辐射光线,能量会逐步衰减(每个圈是球壳)。如下图所示:

    假设距离点光源距离为1(单位距离)的时候,其光强度为(I),现在如果要算出距离点光源距离为(r)处某点的强度(x),因为能量守恒,则有如下式子(强度x球壳面积),则(x=frac{I}{r^2})

    [I*4*pi*1 = x*4*pi*r^2 ]

    (3)漫反射表示 & Shading Point的颜色

    至此,如果知道一个Shading Point距离点光源的距离(r),即可知道有多少强度的光传播到该点附近的单位面积上以及单位面积上接收到的能量情况,就能表示漫反射了,式子如下图:

    式子如下:

    [L_d = k_d (frac{I}{r^2}) max(0,pmb{n}cdot pmb{l}) ]

    由于漫反射是均匀辐射的,所以对于Shading Point而言,其结果与观测点(v)无关。
    假设有一个点光源,Shading Point和光源距离为(r),光源单位距离强度为(I),则到该点的能量强度为(frac{I}{r^2});该点单位面积接收的能量为(cos heta=pmb{n}·pmb{l})(单位向量点乘等于其夹角余弦),但是余弦可能为负数(也就是光线从点的下方射入,无物理意义),所以取值为(max(0,pmb{n}·pmb{l}))。

    如果其材质有颜色,则考虑(k_d)为漫反射系数,如果其值为1,则表示该点完全不吸收能量(完全反射),即为最亮;如果其值为0,则吸收所有能量,即为最暗(例子如下图)。如果将系数(k_d)表示一个三通道的RGB值,每个通道取值为([0,1]),那么就可以在Shading Point上定义某一种颜色。

    3.Specular Term(Blinn-Phong)

    由观察可得,当观测方向(pmb{v})和镜面反射方向(pmb{R})足够接近的时候,就能看到高光,如下图。

    (1)高光条件

    Blinn-Phong模型基于此,认为(pmb{v})(pmb{R})接近也就相当于法线方向(pmb{n})和半程向量(bisector)(pmb{h})接近。两者夹角为(alpha),用其余弦(两者点乘)衡量是否接近,取值为(0,1),越接近则值趋于1,否则趋于0。如下图。

    (2)半程向量(pmb{h})

    这里的半程向量(h)也就是光照方向(pmb{l})观测方向(pmb{v})角平分线方向,这里只要用平行四边形法则然后归一化(长度变为1)处理即为所求,式子如下:

    [egin{aligned} underset{ ext{(半程向量)}}{pmb{h}} &=bisector(pmb{v},pmb{l}) \ &= frac{pmb{v}+pmb{l}}{||pmb{v}+pmb{l}||} end{aligned} ]

    (3)高光表示

    则高光表示如下,(I)为光照强度,其中(k_s)为镜面反射系数,系数越大,亮度越大。同时因为高光通常是白色,所以(k_s)表示颜色通常也为白色。

    [egin{aligned} L_s &=k_s (frac{I}{r^2})max(0,cosalpha)^p\ &=k_s (frac{I}{r^2})max(0,pmb{n}cdot pmb{h})^p end{aligned} ]

    (4)Q1:为什么Blinn-Phong中的高光项表示中没有Shading Point接收到的能量?

    同上面的漫反射比较,会发现这里省略了Shading Point接收到的能量(max(0,pmb{n}·pmb{l}))(实际任何反射都应该考虑),但是因为Blinn-Phong模型是经验模型,更加关注是否能看到高光,即(pmb{h})(pmb{n})是否接近。

    (5)Q2:为什么采用的是半程向量(pmb{h})和法线(pmb{n})的夹角而不是观测方向(pmb{v})和反射方向(pmb{R})的夹角?

    如果采用的观测方向(pmb{v})和镜面反射方向(pmb{R})的夹角的话,其模型称为Phong Reflectance Model,而Blinn-Phong Reflectance Model是对其的改进,因为半程向量容易计算且反射方向不易计算。

    (6)Q3:式子中指数(p)怎么来的?

    虽然向量之间夹角余弦可以表示两者是否接近,但是其容忍度较高,也就是其能看到高光的区间(余弦正值区域)较大,而实际高光是两者方向十分接近的时候才会出现高光,否则是没有高光,所以需要缩小高光区间。下图可以看出(p)增加时,能看到高光的区间逐渐缩小。在Blinn-Phong模型中,(p)的值约为100~200。

    在实际例子中,当(p)增加时,高光会越来越小,如下图:

    4.Ambitent Term(Blinn-Phong))

    此时暂时假设任何一点接收到环境光的强度相同,如下图:

    (1)环境光表示

    式子如下,其中(I_a)为接收到的环境光强度,同理(k_a)为环境光系数,也可以在点上定义一种颜色。

    [L_a=k_aI_a ]

    同时环境光不依赖观测方向(pmb{v})、光照方向(pmb{l})和法线方向(pmb{n}),因为各个方向上强度相同,所以环境光是一个常数,而这就保证物体任何地方至少有一个常数颜色(物体本身不存在一个地方是完全黑的)。

    但是实际的环境光计算,需要用到全局光照,此时用的常数环境光只是便于理解。

    5.模型总结

    环境光是个常数;漫反射与观测方向无关,与法线方向和关照方向有关;以及高光,将三项相加,结果如下图:

    其式子如下:

    [egin{aligned} L&=L_a+L_d+L_s\ &=k_a I_a + k_d (frac{I}{r^2}) max(0,pmb{n}cdot pmb{l})+ k_s (frac{I}{r^2})max(0,pmb{n}cdot pmb{h})^p end{aligned} ]

    二.Shading Frequencies

    有了模型之后,就是对Shading Point进行着色,这里涉及到着色频率的概念。
    着色频率(Shading Frequencies)是指要对哪些点进行着色,例子如下图:

    三个球采用同样的模型,但是着色频率不同结果。

    1.Flat shading(Shade each triangle)

    如第一个球体而言,对每一个面进行着色,对每个三角形面求出其法线(任意两边做叉积),然后根据公式求出结果,即为整个三角形面的颜色。

    2.Grouraud shading(Shade each vertex)

    如第二个球体而言,对每个三角形顶点着色,求出其法线,三角形内容颜色通过插值算出。

    3.Phong shading(Shade each pixel)

    如第三个球体而言,对每个像素进行着色,求出其法线。注意区别Phong shading是个着色频率,Blinn-Phong是个着色模型。

    三种模型结果需要看物体模型的复杂度,顶点数目等,如果物体足够复杂(面足够多,足够小),Flat shading的效果也会很好,如下图:

    4.逐顶点法线计算

    如上图所示,对于一个顶点,会被(N)个三角形所共用,则其法线等于与之相邻的三角形面的法线的平均或者加权平均(权重为三角形的面积),法线最后都需要归一化(化为单位向量),式子如下:

    [N_v=frac{sum^i N_i}{||sum^i N_i||} ]

    5.逐像素法线计算

    如上图,当已知顶点(上图左右两个黑点)的法线之后,中间的法线插值得出,法线最后都需要归一化(化为单位向量)

    三.Graphics(Real-time Rendering)Pipeline

    1.简介

    渲染管线(Graphics Pipeline)或实时渲染管线(Real-time Rendering Pipeline),指的是从某个场景到最后输出图像的整个过程,GPU从硬件上实现了图形管线,如下所示:

    首先输入的空间中的顶点,经过Vertex Processing(包含MVP变换,着色等)之后,将点投影到屏幕空间上;接着Triangle Processing,这些点会在屏幕空间上形成三角形;之后经过Rasterization光栅化(包含采样过程)将三角形离散成为Fragment(片元,类比像素,如果使用MSAA,即为sample);然后在Fragment Processing阶段(包含Z-Buffer测试,着色等),对Fragment着色,此阶段也可以包含在上一步光栅化阶段里;最后在Framebuffer Operations阶段输出图像。

    需要注意的是上面所说的Shading发生在第一步和第四步里面,这是因为对于不同的着色频率发生阶段不同。如果是Grouraud shading对顶点着色就可以在Vertex Processing阶段处理,但是如果是Phong shading对像素着色,就需要等待像素产生,即Fragment Processing阶段处理。

    2.Shader

    图像管线中存在可编程的部分,所以可以人为控制顶点和像色着色部分,而决定顶点和像素如何处理运作的代码即为Shader。Shader会对每个顶点或者像素执行一次(通用执行,即shader只处理一个顶点或者像素 )。如果shader是处理顶点,即为vertex shader(顶点着色器);处理像素(或者fragment),则为fragment shader(片段着色器)或者pixel shader(像素着色器)。

    比如使用OpenGL的着色语言GLSL写的fragment shader例子如下:

    uniform sampler2D myTexture; // 全局变量,表示纹理
    uniform vec3 lightDir; // 全局变量,表示光照方向(认为每个像素都有一个固定的光照方向)
    varying vec2 uv; // per fragment value (interp. by rasterizer)
    varying vec3 norm; //法线(插值)per fragment value (interp. by rasterizer)
    void diffuseShader()
    {
     	vec3 kd;//漫反射系数
     	kd = texture2d(myTexture, uv); //(纹理替代kd)material color from texture
     	kd *= clamp(dot(–lightDir, norm), 0.0, 1.0); // 漫反射部分,Lambertian shading model(OpenGL认为关照方向射向点,所以需要有个负号)
     	gl_FragColor = vec4(kd, 1.0); // output fragment color,gl_FragColor是个固定值,表示fragment的颜色
    } 
    

    四.Text Mapping(纹理映射)

    纹理映射的根本作用是定义任何一点的基本属性,在Blinn-Phong Reflectance Mode中我们使用纹理来替代漫反射系数(k_d)
    前提,任何物体表面都是二维的,也就是说纹理是一张应用于物体表面可变化的图,从而物体表面和纹理上的点就有一一对应关系,如下图:

    前提已知空间中的三角形到纹理的映射关系(美工解决),我们所作的工作是将纹理应用到物体表面。

    1.纹理坐标系((u,v))

    这需要在纹理上定义一个坐标系((u,v)),通常默认规定(u,v)取值为((0,1))(与纹理图片大小无关),如下图:

    2.纹理复用(Tileable Texture)

    同时纹理可以应用在不同的物体表面,并且可以进行复用,如下图(左边是渲染图,右边是((u,v))图)。会发现虽然纹理复用在坐标上是有边界的,但是在物理表面是没有边界的,这种纹理被称为Tileable Texture

    五.Interpolation Across Triangles:Barycentric Coordinates(重⼼坐标)

    1.插值定义

    前面在讲着色频率以及纹理时提到了插值,三角形的顶点有各自不同的属性,插值的目的就是让该属性在三角形内进行平滑的过渡。属性包括纹理映射,顶点颜色,Phong shading中的顶点法线等等。插值通过重心坐标完成。

    2.重心坐标

    重心坐标是针对三角形的,不同的三角形有着不同的重心坐标系统。它可以用三角形(ABC)三个顶点的线性组合来表示三角形平面内任何一点((x,y)),只要三者的系数和为1且为非负数即可,式子如下,从中可以很容易看出三个顶点的重心坐标为,(A(1,0,0))(B(0,1,0))(C(0,0,1))

    [egin{cases} (x,y)=alpha A+eta B+gamma C \ alpha+eta+gamma=1\ alpha >0 ,eta>0,gamma>0 end{cases} ]

    (1)三角形内任意点重心坐标

    对于三角形内任意一点的重心坐标,可以通过面积比算出来,如下图中三角形内部黑点。

    先将其与三个顶点连接,可以得到三个内部三角形,若(A_A)表示A点所对的三角形面积,同理知道(A_B)(A_C),则其重心坐标如下:

    [egin{aligned} alpha&=frac{A_A}{A_A+A_B+A_C}\ eta&=frac{A_B}{A_A+A_B+A_C}\ gamma&=frac{A_C}{A_A+A_B+A_C}\ end{aligned} ]

    如果只知道三顶点的坐标(A(x_A,y_A))(B(x_B,y_B))(C(x_C,y_C)),则可以用上面的公式退出下面的一般表达式:

    [egin{aligned} alpha&=frac{-(x-x_B)(y_C - y_B)+(y - y_B)(x_C - x_B)}{-(x_A - x_B)(y_C - y_B)+(y_A - y_B)(x_C - x_B)}\ eta&= frac{-(x -x_C )(y_A- yC )+(y - y_C )(x_A - x_C )}{-(x_B - x_C )(y_A - y_C )+(y_B - y_C )(x_A - x_C )}\ gamma&= 1- alpha - eta end{aligned} ]

    (2)三角形重心的重心坐标

    而对于三角形的重心而言,三角形的重心是三角形三条中线的交点,重心是三点坐标的平均值,并且重心和三角形3个顶点组成的3个三角形面积相等,可以得到其重心坐标如下:

    [egin{aligned} (alpha,eta,gamma)&=(frac{1}{3},frac{1}{3},frac{1}{3})\ (x,y)&=frac{1}{3}A+frac{1}{3}B+frac{1}{3}C end{aligned} ]

    (3)使用重心坐标进行插值

    使用重心坐标对三角形内的点进行插值,对于需要插值的属性也用重心坐标进行线性组合,如下图所示,三个顶点的属性为(V_A)(V_B)(V_C),这些属性可以是位置,纹理坐标,颜色,法线,深度,材质属性等等。则对于三角形内一点,其重心坐标为((alpha,eta,gamma)),则其属性即为(V=alpha V_A+eta V_B+gamma V_C)

    需要注意的是重心在投影下不能保证不变的。也就是说如果想要对三维空间的某种属性进行插值的话,就应该在三维空间下计算重心坐标系统,不能在投影之后的三角形上计算。比如光栅化阶段三角形已经被投影到屏幕空间上时,此时如果要对深度进行插值,不能直接在该三角形内计算,而是应该在三维空间下的三角形先计算好重心坐标,插值完成后再投影到屏幕空间上。

    六.纹理应用

    至此,可以进行简单的纹理应用,其应用过程如下:

    for each rasterized screen sample (x,y): //(x,y)通常为像素中心
     	(u,v) = evaluate texture coordinate at (x,y)//将点的u,v坐标通过重心坐标进行插值得到
     	texcolor = texture.sample(u,v);//得到u,v坐标对应的纹理颜色
     	set sample’s color to texcolor;//用texcolor代替原本颜色(通常是kd)
    

    1.Texture Magnification(纹理过小情况)

    但是通过上述的过程进行应用之后,会造成Texture Magnification问题,也就是纹理的分辨率过低时被应用到高分辨率物体上(查询纹理坐标时会得到非整数的值)。解决方式如下:

    (1)Nearest

    可以是让一定范围内(u,v)坐标查找距离最近的(比如四舍五入为整数值)的纹理元素(纹素,纹理上的像素,texel),如下图中Nearest所示,会造成块状。

    (2).Bilinear interpolation(双线性插值)

    但是其效果不好,更好的方式是双线性插值法(Bilinear interpolation),具体做法如下:

    对于一个texel,(f(x,y)),如上图红点所示,此时我们想要知道的是纹理在红点处的值。其中黑点表示纹理坐标。首先找到与之相邻的四个点(u_{01})(u_{00})(u_{11})(u_{10}),再找到红点和左下角(u_{00})的水平距离(s)以及竖直距离(t),此时以两个像素之间为单位1的话,可以认为(s,vin(0,1))。如下图所示。

    接下来需要定义线性插值(Linear interpolation)操作(lerp(x,v_0,v_1)),其中(v_0)(v_1)分别表示定义在位置0和1上的值,(xin(0,1))表示需要线性插值属性,即有

    [lerp(x,v_0,v_!)=v_0+x(v_1-v_0) ]

    这样就可以使用(s)对两条水平的边进行插值,即有

    [egin{aligned} u_0&=lerp(s,u_{00},u_{10})\ u_1&=lerp(s,u_{01},u_{11}) end{aligned} ]

    接下来在竖直方向上用(t)(u1)(u_0)进行插值,即可得到红点处的值,即有

    [f(x,y)=lerp(t,u_0,u_1) ]

    以上即为双线性插值方法,即水平和竖直方向上都做插值(顺序不限)。其效果如下图:

    2.Texture Magnification(纹理过大情况)

    当纹理过大时,比如下图中纹理是格子,此时会出现严重的走样问题。

    这是因为屏幕上远处和近处的像素覆盖的纹理范围是各不相同的,如下图。当覆盖区域过大但仍旧采用像素重心进行应用时就会走样。

    之前说过走样产生的本质原因是信号(函数)变换过快(高频率),但是采样速度跟不上。在纹理过大时问题上体现为在一个像素内部可能包含很大的纹理(高频),纹理会发生变化,但是只用一个像素对其进行采样,就会发生走样现象。考虑之前使用MSAA/超采样来反走样,这里同样可以使用,但是开销过大。

    所以考虑避免采样,这就需要知道如何得到一个区域(像素覆盖区域)内的平均值。图形学中使用Mipmap/image pyramid来完成该操作。

    (1)Mipmap定义

    Mipmap的特点是只能用于方形区域,属于近似查询,查询速度快。
    Mipmap就是用一张图生成更多层的图,然后去查询像素所在层数,再进行计算,具体操作如下图:

    其中原始纹理称为第0层(Level 0)纹理,可以通过Mipmap生成更高层纹理,使得第(i)层的分辨率都是第(i-1)层缩小到一半。在这里假设原始纹理的存储为1,则一共使用了(frac{4}{3})存储量,即使用了额外(frac{1}{3})的存储量,如下式子:

    [egin{aligned} 存储量&=lim_{n o infty} left( 1+frac{1}{4}+frac{1}{16}+cdots+frac{1}{4^n} ight)\ &=lim_{n o infty}left(frac{1-left(frac{1}{4} ight)^{left(n+1 ight)}}{1-frac{1}{4}} ight)\ &=frac{4}{3} end{aligned} ]

    (2)使用Mipmap计算Level D

    首先如下图(左边时屏幕空间,右边是纹理坐标),假设此时需要得到四个红点中的左下角红点所代表像素的覆盖面积,可以将与之相邻的像素中心同样投影到纹理坐标上,如上边右图所示。这样就可以算出纹理坐标上该点与相邻点的距离,简单起见从中取一个最大值作为区域的边长,从而近似得到该像素在纹理上覆盖区域的正方形大小,如下图所示。

    计算公式如下,其中(D)为Mipmap的纹理层数。

    [egin{aligned} D&=log_2 L \ \ L&=max left(sqrt{left(frac{du}{dx} ight)^2+left(frac{dv}{dx} ight)^2},sqrt{left(frac{du}{dy} ight)^2+left(frac{dv}{dy} ight)^2} ight) end{aligned} ]

    得到了近似区域的大小(L^2),即可知道区域对应的层数(log_2),比如区域为4x4大小,则其必定在第二层上,则在第二层上会变成一个像素。

    (3)Trilinear interpolation(三线性插值)

    但是将需要查询的Mipmap层数进行可视化如下图,会发现此时层层之间的查询不连续,因为此时层数都是整数值,所以需要使用插值来进行平滑的过度。

    考虑在层与层之间进行插值,首先在第(D)层和第(D+1)层上分别做双线性插值得到对应结果,再利用这两个值在层与层之间进行插值(非水平非竖直方向,层与层),此为Trilinear interpolation(三线性插值)。

    其结果图如下,层与层之间有很好的过度:

    (4)Anisotropic Filter(各向异性过滤)

    现在回到最开始的网格图,应用Mipmap结果如下图,远处出现了严重的Overblur(过度模糊)现象。

    这是因为屏幕空间上的像素对应到纹理上覆盖的区域不一定都是正方形,如下图所示。采用Mipmap,这些长条状的区域都会使用一个正方形区域近似,就会求得一个更大区域的平均值,从而过度模糊。

    Anisotropic Filter(各向异性过滤)可以部分解决这个问题,它和Mipmap相比,它计算了不同长宽比的纹理层数,如下图所示,Mipmap计算的是对角线上,而Anisotropic Filter还计算了周围的对象,这就可以解决矩形状区域的快速查询,但是总共的存储量是原本的3倍(收敛极限)。

    3.其他的纹理应用

    (1)Environment Map/Lighting

    使用纹理记录环境光,然后再进行渲染。还可以将环境光记录在球面上(Spherical Environment Map),如下图。

    但是这样将其展开的时候就会发现有扭曲现象,如下图:

    改进方式是Cube Map,也就是给球面一个包围盒,如下图,则原来存储在球面上的关照信息衍生到包围盒上,这样就会得到一个六张图片组合的环境光,如下图。

    (2)凹凸/法线贴图(Bump Mapping)

    纹理除了可以定义颜色之外,还可以定义其他不同的属性,比如定义在一个表面上任意一点的相对于基础表面上沿着法线方向的相对高度,从而避免使用大量的三角形来定义部分复杂(凹凸不平)的几何形体,如下图:

    通过凹凸贴图这样就可以在不改变几何形体的情况下来改变着色结果,将像素的法线进行(Perturb)扰动(仅仅为了着色计算而使用),也就是相对于平面的高度变化,相当于改变了法线(实际没有改变物体原本的法线),比如下图中(p)点通过凹凸贴图会被认为移动到了上面的点。

    <1>接着考虑凹凸贴图上的法线如何变化?(先在贴图上定义切线,通过切线得到对应法线)

    简单起见先考虑一维贴图/Flatland上的变换情况,如下图情况,原本平面是平的,蓝色线是由凹凸贴图定义得到的,原本的表面法线在(p)点是(n(p)=(0,1)),这里是一个临时的坐标系,然后用差分方式(用相邻两点的高度差进行计算)近似求出贴图上(p)点的切线(也就是该点的导数),如下式子,其中(c)为凹凸贴出的影响因子。

    [d_p=c*[h_{p+1}-h_p] ]

    导数相当于水平方向上移动一个单位距离,竖直方向上移动(d_p),则此处的切线用向量表示为((1,d_p)),则扰动后的法线应该和切线垂直,即将其逆时针旋转(90^o),得到对应法线为((-d_p,1).normalized()),最后需要归一化。

    将上述结果推广到实际情况,对于二维贴图,假设原本法线为(n(p)=(0,0,1)),再求出(u,v)坐标上的对应的导数,如下:

    [frac{d_p}{d_u}=c_1*[h_{(u+1)}-h_u]\ frac{d_p}{d_v}=c_2*[h_{(v+1)}-h_v] ]

    类比flatland情况,推导略,则得到法线如下:

    [(-frac{d_p}{d_u},-frac{d_p}{d_v},1).normalized() ]

    (3)位移贴图(Displacement mapping)

    于凹凸贴图相比较,位移贴图实际上通过顶点的移动改变了物体的几何形状,这对模型三角形要求精细,两者结果比较如下图:

  • 相关阅读:
    Win10 ntoskrnl.exe蓝屏解决
    Log POST Data in Nginx
    MACOS关闭指定端口
    获取Skype用户IP地址
    禁止windows10带来的三大隐患问题
    各种语言一句话反弹shell
    2015阿里巴巴安全峰会PPT
    HTTPS反向代理嗅探
    利用arpspoof和urlsnarf 进行ARP嗅探
    收集的几个存在漏洞的程序
  • 原文地址:https://www.cnblogs.com/FlyerBird/p/13245807.html
Copyright © 2011-2022 走看看