FT 在图形渲染中的应用：基于 FFT 的海浪模拟

2022-12-18 08:11 354人阅读评论(0)

接上文：FT 在图像处理中的应用

五、一个大型案例：基于 FFT 的海浪模拟

前置：

5.1 FFT 海洋公式：二维 IDFT

https://tore.tuhh.de/bitstream/11420/1439/1/GPGPU_FFT_Ocean_Simulation.pdf

https://www.youtube.com/watch?v=ClW3fo94KR4

既然任意连续且收敛的函数都可以分解为无数个不同频率、不同幅值的正、余弦信号，且这个性质可以扩展到二维，即任意一个二维图像都可以分解为无数个复平面波的叠加

那么海浪作为“波”的典型，是否也可以将其拆成任意方向不同频率强度的基波的叠加？当然没有问题：令为时刻，位置上海面的高度，其构建公式就为

可以看出，这正是二维离散傅里叶变换的逆变换(IDFT)版本，只不过多了一个量度即时间，考虑到 DFT 形式下任意将对应所有可能的复平面波，就能顺理成章的推出波矢量

其中常量为海平面大小，采样点满足，N 为采样点数量，当然也可以做个变换，即将范围映射到，此时

本质不变

代入上面的海洋方程，并将所有向量全部展开，为方便可以直接任就有

当然也可以任，这在计算方式上是相同的

5.1.1 二维 IFFT

前面已经详细介绍过了一维 FFT 的原理和计算过程，这里同样要扩展二维以得到最终的，其实看这个二维 IDFT 式子，你就会发现结论很明显了：

它其实是可以拆成两个一维的 IDFT 的：

那么最终我们计算海洋公式的步骤就应如下：

计算生成频谱，如果是交给 GPU 计算，可以将结果存储在一张 NxN 的纹理中（关于的具体内容及计算方式可以参考下一节 5.2 的内容），t 可以理解成第几帧，必然每一帧这个图像都会不同
执行 N 次水平 FFT：，上一步的纹理作为输入，输出结果作为下一步的输入
执行 N 次纵向 FFT：，得到的结果再进行最后的符号计算：即乘上，搞定

到此，我们就拿到了一张高度图，可以用于后续的顶点动画及渲染着色了

5.2 海浪统计模型

事实上，这种基于 FFT 或者说基于频谱逆向推出海洋水面高度的方法，本质上是一种统计模型方法，这里面涉及不少其他领域内的知识，例如海洋物理学，又或是对现实海洋的实时观测统计与实验，我们当然不用深入了解这部分的内容，只要拿对应的统计结果来用就可以了

上面的内容还有一个很重要的东西还未知，就是用于表达海洋波的频率的集合、IDFT 中的频域函数，搞定了这个之后，我们才能通过 IDFT 算法将频谱转化为时域的值，从而得到海洋随时间变化的高度图

~~公式都给你写好了，照着抄就行~~

5.2.1 菲利浦频谱(Phillips spectrum)

菲利浦频谱来自海洋动力学的成果，是一个反馈了真实海洋的经验模型，也是这里一切的核心，当然这个模型并不唯一

对于菲利浦频谱，其中

，是风速，为风向，为引力常数
为的幅值，即频率大小

从论文得到的信息是：菲利浦频谱主要由两个单独部分组成，分别为

非定向波谱，由海洋统计分析得出
方向扩展函数，体现的是风向对波最终幅值的贡献

5.2.2 海洋初始状态

回到前面的，可以先考虑它的初始状态，也就是 t = 0 的时刻的振幅有：

和为两个相互独立的服从均值为0，标准差为1的高斯随机数，实部表初始幅度，虚部表初始相位，而正是前面的方向波谱

这里再介绍另一个案例，其实后面的 FFT 实现也正是借鉴的这篇视频，这个案例中额外考虑了频谱到波的离散转化：

5.2.3 海洋公式最终形式

最后就是把时间参数加上：

，这里又有

角频率与波长的频散关系满足，其考虑了引力风波和重力场的关系，当然如果额外考虑水深 h，其频散关系
是的共轭复数

到此为止，整个 IDFT 公式就明朗了，然后就是套用上面 IFFT 的步骤

5.3 基于 GPGPU 的 Ocean IFFT 计算

通过上面的内容，可以预见的是生成高度图会有大量的数学计算，尽管这已经是算法加速后的结果了，因此大多数情况下都会把这部分计算任务交给更擅长并行计算的 GPU 而非 CPU

在 Unity 中，往往使用 ComputeShader 来实现

5.3.1 初始频谱计算

先是计算菲利浦频谱以生成对应纹理，看上去最陌生的公式其实直接抄作业就行了，其中作为纹理坐标，其余除了风向都为常量。很明显，游戏中我们不太可能每帧去改变风向，因此风向也可以视作常量（由美术或策划去配置），在这种情况下的计算理论可以离线或者只在风向改变时做一次，不需要每帧去重复计算

然后就是高斯随机数生成，这个只是为了得到不同的初始状态，因此同样可以离线计算：


  
   
    
     
    
    
     
      //如果本地已经存在 GuassTexture，则不需要生成了，直接读取，否则在 Awake 时生成一份
     
    
   
    
     
    
    
     
      private void Awake()
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
     
          gaussianNoise = GetNoiseTexture(size);
     
    
   
    
     
    
    
     
      }
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      Texture2D GetNoiseTexture(int size)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
         
      string filename = 
      "GaussianNoiseTexture" + size.ToString() + 
      "x" + size.ToString();
     
    
   
    
     
    
    
     
          Texture2D noise = Resources.Load<Texture2D>(
      "GaussianNoiseTextures/" + filename);
     
    
   
    
     
    
    
         
      return noise ? noise : GenerateNoiseTexture(size, 
      true);
     
    
   
    
     
    
    
     
      }
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
     
      float NormalRandom()
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
         
      return Mathf.Cos(
      2 * Mathf.PI * Random.
      value) * Mathf.Sqrt(
      -2 * Mathf.Log(Random.
      value));
     
    
   
    
     
    
    
     
      }
     
    
   
    
     
    
    
         
     
    
   
    
     
    
    
     
      Texture2D GenerateNoiseTexture(int size, bool saveIntoAssetFile)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
     
          Texture2D noise = 
      new Texture2D(size, size, TextureFormat.RGFloat, 
      false, 
      true);
     
    
   
    
     
    
    
     
          noise.filterMode = FilterMode.Point;
     
    
   
    
     
    
    
         
      for (
      int i = 
      0; i < size; i++)
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
             
      for (
      int j = 
      0; j < size; j++)
     
    
   
    
     
    
    
     
              {
     
    
   
    
     
    
    
     
                  noise.SetPixel(i, j, 
      new Vector4(NormalRandom(), NormalRandom()));
     
    
   
    
     
    
    
     
              }
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
     
          noise.Apply();
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
         
      //此处尝试将 GuassTexture 保存到本地，代码略
     
    
   
    
     
    
    
         
      return noise;
     
    
   
    
     
    
    
     
      }

GPU Gems3:Chapter 37 中提到过如何生成满足服从均值为0，标准差为1的高斯随机数，公式如下：

其中为均匀分布的随机数

然后两者相组合，就能得到：


  
   
    
     
    
    
     
      public void CalculateInitials(WavesSettings wavesSettings, float lengthScale)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
         
      //设置计算频谱的各项参数，例如风速、风向、重力加速度常量等
     
    
   
    
     
    
    
     
          wavesSettings.SetParametersToShader(initialSpectrumShader, KERNEL_INITIAL_SPECTRUM, paramsBuffer);
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
         
      //将计算的 H0k 保存到一张纹理中
     
    
   
    
     
    
    
     
          initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, H0K_PROP, buffer);
     
    
   
    
     
    
    
         
      //计算角频率与波长的频散关系 w 保存到另一张纹理中
     
    
   
    
     
    
    
     
          initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, PRECOMPUTED_DATA_PROP, precomputedData);
     
    
   
    
     
    
    
         
      //传入高斯随机数纹理，用于计算 H0K
     
    
   
    
     
    
    
     
          initialSpectrumShader.SetTexture(KERNEL_INITIAL_SPECTRUM, NOISE_PROP, gaussianNoise);
     
    
   
    
     
    
    
         
      //交给 GPGPU 开始计算
     
    
   
    
     
    
    
     
          initialSpectrumShader.Dispatch(KERNEL_INITIAL_SPECTRUM, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 
      1);
     
    
   
    
     
    
    
      
     
    
   
    
     
    
    
         
      //计算得到 H0k 后，把共轭结果 H*0K 也存下来
     
    
   
    
     
    
    
     
          initialSpectrumShader.SetTexture(KERNEL_CONJUGATE_SPECTRUM, H0_PROP, initialSpectrum);
     
    
   
    
     
    
    
     
          initialSpectrumShader.SetTexture(KERNEL_CONJUGATE_SPECTRUM, H0K_PROP, buffer);
     
    
   
    
     
    
    
     
          initialSpectrumShader.Dispatch(KERNEL_CONJUGATE_SPECTRUM, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 
      1);
     
    
   
    
     
    
    
     
      }

中间用到了一个专门用于计算频谱函数的 ComputeShader，里面就是纯计算，并且由于计算模型并不唯一，具体代码就不贴了

5.3.2 任意时刻频谱计算

关于的代码实现更简单，毕竟前面已经得到了以及（该图片来源于 youtube《Ocean waves simulation with Fast Fourier transform》，后同）：


  
   
    
     
    
    
     
      float4 wave = WavesData[
      id.xy];
     
    
   
    
     
    
    
     
      float phase = wave.w * Time;
     
    
   
    
     
    
    
     
      float2 exponent = float2(cos(phase), sin(phase));
     
    
   
    
     
    
    
     
      float2 h = ComplexMult(H0[
      id.xy].xy, exponent)
     
    
   
    
     
    
    
     
          + ComplexMult(H0[
      id.xy].zw, float2(exponent.x, -exponent.y));
     
    
   
    
     
    
    
     
      float2 ih = float2(-h.y, h.x);

不过考虑到海面并不只有高度变化，还有水平方向的流动，后续只有一张 y 轴的高度图还是不够的，还需要 x 轴和 z 轴的偏移图，这里可以先根据获取水平方向的频谱

Gerstner 波为了展现波峰尖角，也同样对 xz 平面做了变换

前面代码将一起塞到了一张4通道纹理中，这里继续拿来用：


  
   
    
     
    
    
     
      float2 displacementX = ih * wave.x * wave.y;
     
    
   
    
     
    
    
     
      float2 displacementY = h;
     
    
   
    
     
    
    
     
      float2 displacementZ = ih * wave.z * wave.y;

5.3.3 海洋法线

其实到这这一步应该就可以准备 IFFT 计算了，不过中间还有一个东西没有考虑，那就是法线图的生成，因为你要计算光照，法线图是必须的，实时的海浪必然意味着法线图也需要实时计算

这里关于法线计算有两种方案

得到高度图后，可以通过差分方法来求法线，即对变换后的新顶点，求出与其相邻的两组顶点和，交叉连线后得到两个向量，叉乘得出法线，
通过偏导计算求出曲面梯度，然后再根据梯度得到法线

方案①逻辑非常简单，可以直接在着色器中实时计算得到法线，缺点也很明显，那就是仅能得到近似法线，这个结果是不会准确的，所以后续若想要得到一个比较好的光照效果，就只能采取方案②

梯度的本意是一个向量，它的方向与取得最大方向导数的方向一致，而它的模为方向导数的最大值，也可以理解为函数在该点处沿着该方向变化最快，变化率最大

设函数在平面区域 D 内具有一阶连续偏导数，则对于每一点，都可定出一个向量，该向量称为函数在点的梯度，记作

可以证明：

对于高度图，其梯度满足

，其中偏导


  
   
    
     
    
    
     
      float2 displacementY_dx = ih * wave.x;
     
    
   
    
     
    
    
     
      float2 displacementY_dz = ih * wave.z;

不过还不够，我们还有两张水平方向的偏移图：

此时不止高度会发生变化，延 x 方向和 z 方向也有一段偏移，这段偏移是

和，对应的偏导


  
   
    
     
    
    
     
      float2 displacementX_dx = -h * wave.x * wave.x * wave.y;
     
    
   
    
     
    
    
     
      float2 displacementZ_dz = -h * wave.z * wave.z * wave.y;

因此，考虑了水平偏移后，真正的海浪梯度就为：

看上去这个公式很复杂，其实本质上就是 (高度图对 x 偏导 / 扭曲挤压后的新海面对 x 偏导, 高度图对 z 偏导 / 扭曲挤压后的新海面对 z 偏导)

高度图很好理解，但是平面延 x 和 z 方向的偏移过程可能不是特别好理解：这里可以沿用一张 wiki 上的图片来做参考：

可以看到，在应用水平偏移后，整个海平面就不再是一个线性的平面了，其坐标会被扭曲，甚至会出现挤压（有向面元在变换后面积为负数），

这样再看上面那个极其复杂的海浪梯度向量，其实分子部分就是

分母部分即是对上图右部分的新曲面求偏导，即

也可以以实际应用为例子，直接看最终海洋的顶点变换表现，当你没有应用水平偏移时，最终顶点动画的效果如左：

可以看出顶点网格水平方向上其实是静止的，只有高度在发生变化，而如果应用了水平方向的偏移，最终顶点动画效果如右

差别明显，如果我们把视角调成垂直往下，就可以直接观察到变换后海面的顶点水平分布情况：

搞定梯度后，想要计算法线就很简单了，只需要拿上向量减去梯度向量，得到的就是法向量，这样下来对于空间中海面上的一点

，即

其精确法线（未考虑归一化）就为

搞定

5.3.4 IFFT 蝶形图生成

第三章已经对 FFT 蝶形算法做过了详细的解析，一个 N = 8 的蝶形算法如下图：

对于 IDFT ，，令，其满足性质：

其中，


  
   
    
     
    
    
     
      public void IFFT2D(RenderTexture input, RenderTexture buffer, bool outputToInput = false, bool scale = true, bool permute = false)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
         
      int logSize = (
      int)Mathf.Log(size, 
      2);
     
    
   
    
     
    
    
         
      bool pingPong = 
      false;
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_PRECOMPUTED_DATA, precomputedData);
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_BUFFER0, input);
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_HORIZONTAL_STEP_IFFT, PROP_ID_BUFFER1, buffer);
     
    
   
    
     
    
    
         
      for (
      int i = 
      0; i < logSize; i++)
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
     
              pingPong = !pingPong;
     
    
   
    
     
    
    
     
              fftShader.SetInt(PROP_ID_STEP, i);
     
    
   
    
     
    
    
     
              fftShader.SetBool(PROP_ID_PINGPONG, pingPong);
     
    
   
    
     
    
    
     
              fftShader.Dispatch(KERNEL_HORIZONTAL_STEP_IFFT, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 
      1);
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_PRECOMPUTED_DATA, precomputedData);
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_BUFFER0, input);
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_VERTICAL_STEP_IFFT, PROP_ID_BUFFER1, buffer);
     
    
   
    
     
    
    
         
      for (
      int i = 
      0; i < logSize; i++)
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
     
              pingPong = !pingPong;
     
    
   
    
     
    
    
     
              fftShader.SetInt(PROP_ID_STEP, i);
     
    
   
    
     
    
    
     
              fftShader.SetBool(PROP_ID_PINGPONG, pingPong);
     
    
   
    
     
    
    
     
              fftShader.Dispatch(KERNEL_VERTICAL_STEP_IFFT, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 
      1);
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
         
     
    
   
    
     
    
    
         
      if (pingPong)
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
     
              Graphics.Blit(buffer, input);
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
     
          fftShader.SetInt(PROP_ID_SIZE, size);
     
    
   
    
     
    
    
     
          fftShader.SetTexture(KERNEL_PERMUTE, PROP_ID_BUFFER0, outputToInput ? input : buffer);
     
    
   
    
     
    
    
     
          fftShader.Dispatch(KERNEL_PERMUTE, size / LOCAL_WORK_GROUPS_X, size / LOCAL_WORK_GROUPS_Y, 
      1);
     
    
   
    
     
    
    
     
      }

如果递归计算上面按所有的 X, G, H，总体复杂度为 O(NlogN)，例如 N = 8 时，从上图可以看出其递归深度为3，每次计算量都为8，即总共 3x8=24 个蝶形计算单元：

一次蝶形计算单元如上，由于我们可能需要进行多次 IFFT，但事实上由于每次 IFFT 的大小（采样次数）N 是固定的，因此可以预计算所有的 ，并将其复数结果放入一张 logN x N 大小的纹理的 R 和 G 通道中（纹理坐标 x, y 取值从0开始），令 ，则对于纹理 位置的值 满足：

还有两个通道 B 和 A，用于放入图中的上一层 IFFT 的索引 a 和 b，a 和 b 的计算公式如下：

对于 N = 8 的情况下，a, b, k 和 d 的取值即如下，纵轴对应 y = 0~8，横轴对应 x = 0~2

x = 0：a = 0 1 2 3 0 1 2 3 b = 4 5 6 7 4 5 6 7 k = 0 0 0 0 4 4 4 4 d = 4
x = 1：a = 0 1 4 5 0 1 4 5 b = 2 3 6 7 2 3 6 7 k = 0 0 2 2 4 4 8 8 d = 2
x = 2：a = 0 2 4 6 0 2 4 6 b = 1 3 5 7 1 3 5 7 k = 0 1 2 3 4 5 6 7 d = 1

对于实际应用的例子中 N = 256 的情况，生成的蝶形纹理如下：


  
   
    
     
    
    
     
      void PrecomputeTwiddleFactorsAndInputIndices(uint3 
      id : SV_DispatchThreadID)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
     
          uint b = Size >> (
      id.x + 
      1);
     
    
   
    
     
    
    
     
          float2 mult = 
      2 * PI * float2(
      0, 
      1) / Size;
     
    
   
    
     
    
    
     
          uint i = (
      2 * b * (
      id.y / b) + 
      id.y % b) % Size;
     
    
   
    
     
    
    
     
          float2 twiddle = ComplexExp(-mult * ((
      id.y / b) * b));
     
    
   
    
     
    
    
     
          PrecomputeBuffer[
      id.xy] = float4(twiddle.x, twiddle.y, i, i + b);
     
    
   
    
     
    
    
     
          PrecomputeBuffer[uint2(
      id.x, 
      id.y + Size / 
      2)] = float4(-twiddle.x, -twiddle.y, i, i + b);
     
    
   
    
     
    
    
     
      }

这个无论是 k 的计算，还是最后索引 a 的计算，可能单看数字不是特别明白，但是没有关系，只要代入到后面的计算场景就会立刻清晰：

如果采样递归的手段计算 FFT，那么逻辑写起来会相对简单一点，因为它是顺着思路的，之前就有举过一道多项式乘法的例子，并且给出了代码，然而由于我们想要 GPU 帮我们并行计算，就不好递归去做，这个时候只能逆过来算，也就是按照蝶形图从左至右（stage 1~3 的顺序）的计算每一个单元：这个时候再看我们关于索引 a, b 以及 k 的计算就好理解多了

当然每一层计算完毕后，索引都会按照该层的计算顺序更新排列

不仅如此，由于我们计算每一个蝶形单元时，只依赖上一个层(stage)的结果，因此每一层的运算都可以并行处理，效率极高：


  
   
    
     
    
    
     
      void HorizontalStepInverseFFT(uint3 
      id : SV_DispatchThreadID)
     
    
   
    
     
    
    
     
      {
     
    
   
    
     
    
    
     
          float4 data = PrecomputedData[uint2(Step, 
      id.x)];
     
    
   
    
     
    
    
     
          uint2 inputsIndices = (uint2)data.ba;
     
    
   
    
     
    
    
         
      if (PingPong)
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
     
              Buffer1[
      id.xy] = Buffer0[uint2(inputsIndices.x, 
      id.y)]
     
    
   
    
     
    
    
     
                  + ComplexMult(float2(data.r, -data.g), Buffer0[uint2(inputsIndices.y, 
      id.y)]);
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
         
      else
     
    
   
    
     
    
    
     
          {
     
    
   
    
     
    
    
     
              Buffer0[
      id.xy] = Buffer1[uint2(inputsIndices.x, 
      id.y)]
     
    
   
    
     
    
    
     
                  + ComplexMult(float2(data.r, -data.g), Buffer1[uint2(inputsIndices.y, 
      id.y)]);
     
    
   
    
     
    
    
     
          }
     
    
   
    
     
    
    
     
      }

对所有行来一次横向 FFT 后，就是所有列的纵向 FFT，由于对应的计算方式几乎完全一致，就不再做详细介绍了

5.3.5 收尾

到此整个 FFT Ocean 的物理计算流程就要到尾声了，经过 IFFT 计算后，我们就能通过后续纹理合并，以及向量计算拿到：

一张海浪高度图
水平偏移图
法线图

其中①和②可以合到一张纹理中，然后运用到顶点变换：


  
   
    
     
    
    
     
      float3 displacement = 
      0;
     
    
   
    
     
    
    
     
      displacement += tex2Dlod(_Displacement_c0, worldUV / LengthScale0) * lod_c0;
     
    
   
    
     
    
    
     
      v.vertex.xyz += mul(unity_WorldToObject, displacement);

法线图的计算则如下


  
   
    
     
    
    
     
      float4 derivatives = tex2D(_Derivatives_c0, IN.worldUV / LengthScale0);
     
    
   
    
     
    
    
     
      float2 slope = float2(derivatives.x / (
      1 + derivatives.z),
     
    
   
    
     
    
    
     
          derivatives.y / (
      1 + derivatives.w));
     
    
   
    
     
    
    
     
      float3 worldNormal = normalize(float3(-slope.x, 
      1, -slope.y));

搞定，效果放在了文章最前面，代码的算法部分实际上是直接借鉴的 gasgiant/FFT-Ocean，当然为了能让高配手机上也能运行，以及要做非真实感的海洋，还是要做不少工作的，例如着色，水体交互、海浪泡沫、性能优化等等。考虑到本篇的原意只是介绍 FT 及其在游戏开发领域的应用，就先在这画上句号吧

转载：https://blog.csdn.net/Jaihk662/article/details/128284952

查看评论

小言_互联网的博客

小言_互联网的博客

个人资料

文章分类

文章存档

阅读排行

评论排行

推荐文章