shader编程-三维场景下SDF建模，对模型进行扭曲、弯曲、裁剪、掏空操作（WebGL-Shader开发基础12）_三维 sdf

1. demo效果

如上图所示，第一张是对方形圆柱进行扭曲操作的效果图，第二张图是对扁平立方体面进行弯曲操作的效果图，第三张图是从内到外依次是原物体，裁剪后的物体，裁剪后掏空物体的效果，从左到右依次对圆柱、球体、甜圈圈进行同样变换的效果

2. 实现要点

先回顾一下之前学习的SDF建模相关的一些操作，首先学习了基本变换旋转、缩放、平移，之后学习了通过布尔运算(交集、并集、差集)组合模型 ，进而拓展了平滑布尔运算组合模型，今天解锁一点点新的技能，对模型进行扭曲、弯曲、裁剪、掏空操作

2.1 模型扭曲

通过前面的demo效果展示知道，我们首先绘制了一个长条状的立方体，然后又绘制了一个以同样立方体扭曲后的物体，实现的核心思想就是变换用来绘制扭曲立方体的坐标，即依据y分量的变化使x分量和z分量转圈圈，具体代码如下

vec3 opTwistY( in vec3 p ,float k)
{

  float c = cos(k*p.y);
  float s = sin(k*p.y);
  mat2  m = mat2(c,-s,s,c);
  p.xz *=m;
  return p;
}
1
2
3
4
5
6
7
8
9

如果你想物体绕z轴或x轴扭曲，两种办法，一是先按y轴扭曲然后将物体旋转，二是将上面的函数中的分量替换为对应的坐标轴

扭曲的坐标处理完接着就是像之前一样绘制图形，为了方便使用我们将绘制过程封装成了一个函数，具体如下

vec2 twistBox (vec4 pos){

  //绘制原立方体
  pos.x += 2.5; 
  float box1 = sdBox(pos.xyz,vec3(0.4,1.2,0.4),0.06);//方块
  vec2 mBox1 = vec2(box1,1.0);

  //绘制扭曲立方体
  pos.x -= 5.0; 
  pos.xyz=opTwistY(pos.xyz,3.0);//扭曲
  float box2 = sdBox(pos.xyz,vec3(0.4,1.2,0.4),0.06);//方块
  vec2 mBox2 = vec2(box2,1.0);

  vec2 res = opU(mBox1,mBox2);

  return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

调用过程如下

vec2 res =  vec2(p.y,0.0);//地面

vec3 pos = p-vec3(0,2,5);//确定模型的中心
vec4 oriPos = vec4(pos,1.0);//转为其次坐标

//扭曲
vec2 twistBox = twistBox(oriPos);

res = opU(res,twistBox);//扭曲立方体
1
2
3
4
5
6
7
8
9

2.2 模型弯曲

弯曲也是同样的思想，通过改变绘制弯曲模型的坐标来实现，返回弯曲坐标的函数如下

vec3 opCheapBend(in vec3 p ,float k)
{
  float c = cos(k*p.x);
  float s = sin(k*p.x);
  mat2  m = mat2(c,-s,s,c);
  vec3  q = vec3(m*p.xy,p.z);
  return q;
}
1
2
3
4
5
6
7
8

像上一次一样，使用处理好的坐标绘制图形即可，这次也将绘制过程封装成了一个函数，具体如下

vec2 cheapBendBox (vec4 pos){

  //绘制原立方体面
  pos.x += 2.5; 
  float box1 = sdBox(pos.xyz,vec3(1.0,0.2,0.6),0.06);//方块
  vec2 mBox1 = vec2(box1,1.0);

  //绘制弯曲立方体面
  pos.x -= 5.0; 
  pos.xyz=opCheapBend(pos.xyz,0.5);//弯曲
  float box2 = sdBox(pos.xyz,vec3(1.0,0.2,0.6),0.06);//方块
  vec2 mBox2 = vec2(box2,1.0);

  vec2 res = opU(mBox1,mBox2);

  return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

调用过程如下

vec2 res =  vec2(p.y,0.0);//地面

vec3 pos = p-vec3(0,2,5);//确定模型的中心
vec4 oriPos = vec4(pos,1.0);//转为其次坐标

//扭曲
vec2 twistBox = twistBox(oriPos);

//弯曲
vec2 cheapBendBox = cheapBendBox(oriPos);
//res = opU(res,twistBox);//扭曲立方体
res = opU(res,cheapBendBox); //弯曲立方体
1
2
3
4
5
6
7
8
9
10
11
12

2.3 模型裁剪与掏空

2.3.1 球体裁剪与掏空

掏空函数

//掏空
float opOnion(float sdf, in float thickness )
{
  return abs(sdf)-thickness;
}
1
2
3
4
5

绘制过程是，我们首先绘制了一个完整的球体，然后绘制了一个裁剪的球体，最后绘制一个掏空且裁剪的球体，之所以也裁剪是因为不裁剪的画我们看不到内部结构，具体如下

vec2 onionSphere (vec4 pos){

  //绘制球体
  vec4 pos1 = pos;
  pos1.z -= 6.0; 
  float sphere1 = sdSphere(pos1.xyz,0.6);
  vec2 mSphere1 = vec2(sphere1,2.0);

  //绘制裁剪球体
  vec4 pos2 = pos;
  pos2.z -= 2.5; 
  float sphere2 = sdSphere(pos2.xyz,0.6);
  sphere2 = max( sphere2, pos2.y );//裁剪

  vec2 mSphere2 = vec2(sphere2,2.0);

  //绘制裁剪掏空球体
  vec4 pos3 = pos;
 
  float sphere3 = sdSphere(pos3.xyz,0.6);

  sphere3 = opOnion(sphere3,0.01);//掏空
  sphere3 = max( sphere3, pos3.y );//裁剪
  vec2 mSphere3 = vec2(sphere3,2.0);

  vec2 res = opU(mSphere1,mSphere2);
  res = opU(res,mSphere3);

  return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

关于裁剪，不知你有没有联想到，这里的裁剪效果是裁掉了y轴大于0 的部分，如果想裁掉x轴或z上的大于0的部分，只需把max( sphere3, pos3.y );//裁剪此处的分量y替换,如果你想要留下y轴上部分，即demo中裁掉和留下的部分互换，那么在y分量前添加负号即可，max( sphere3, -pos3.y );//裁剪，其他轴中的裁剪同理

2.3.2 圆柱裁剪与掏空

圆柱的实现与球体的思路一样，只不过把球体换成圆柱，如下

vec2 onionCylinder (vec4 pos){

  //绘制圆柱
  vec4 pos1 = pos;
  pos1.z -= 6.0; 
  float cylinder1 = sdCylinder(pos1.xyz,0.4,0.6);
  vec2 mCylinder1 = vec2(cylinder1,1.0);

  //绘制裁剪圆柱
  vec4 pos2 = pos;
  pos2.z -= 2.5; 
  float cylinder2 = sdCylinder(pos2.xyz,0.4,0.6);
  cylinder2 = max( cylinder2, pos2.y );//裁剪

  vec2 mCylinder2 = vec2(cylinder2,1.0);

  //绘制裁剪掏空圆柱
  vec4 pos3 = pos;

  float cylinder3 = sdCylinder(pos3.xyz,0.4,0.6);

  cylinder3 = opOnion(cylinder3,0.01);//掏空
  cylinder3 = max( cylinder3, pos3.y );//裁剪
  vec2 mCylinder3 = vec2(cylinder3,1.0);

  vec2 res = opU(mCylinder1,mCylinder2);
  res = opU(res,mCylinder3);

  return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

2.3.3 甜圈圈裁剪与掏空

甜圈圈裁剪与掏空也是同样的思路，绘制过程如下

vec2 onionTorus (vec4 pos){

  //绘制甜圈圈
  vec4 pos1 = pos;
  pos1.z -= 6.0; 
  float torus1 = sdTorus(pos1.xyz,vec2(0.6,0.2));
  vec2 mTorus1 = vec2(torus1,3.0);

  //绘制裁剪甜圈圈
  vec4 pos2 = pos;
  pos2.z -= 2.5; 
  float torus2 = sdTorus(pos2.xyz,vec2(0.6,0.2));
  torus2 = max( torus2, pos2.y );//裁剪

  vec2 mTorus2 = vec2(torus2,3.0);

  //绘制裁剪掏空甜圈圈
  vec4 pos3 = pos;

  float torus3 = sdTorus(pos3.xyz,vec2(0.6,0.2));

  torus3 = opOnion(torus3,0.01);//掏空
  torus3 = max( torus3, pos3.y );//裁剪
  vec2 mTorus3 = vec2(torus3,3.0);

  vec2 res = opU(mTorus1,mTorus2);
  res = opU(res,mTorus3);

  return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

3. demo代码

继续，可直接跑起来的代码

<body>
  <div id="container"></div>
  <script src="http://www.yanhuangxueyuan.com/versions/threejsR92/build/three.js"></script>
  <script>
    var container;
    var camera, scene, renderer;
    var uniforms;
    var vertexShader = `
      void main() {
        gl_Position = vec4( position, 1.0 );
      } 
    `
    var fragmentShader = `
    #ifdef GL_ES
    precision mediump float;
    #endif
    uniform float u_time;
    uniform vec2 u_mouse;
    uniform vec2 u_resolution;

    const int MAX_STEPS = 100;//最大步进步数
    const float MAX_DIST = 100.0;//最大步进距离
    const float SURF_DIST = 0.01;//相交检测临近表面距离

    //绕z轴旋转矩阵
    mat4 rotZ(float a) {
      return mat4(cos(a),-sin(a),0.0,0.0,
                  sin(a),cos(a),0.0,0.0,
                  0.0,0.0,1.0,0.0,
                  0.0,0.0,0.0,1.0
                  );
    }
    
    //绕x轴旋转矩阵
    mat4 rotX(float a) {
      return mat4(1.0,0.0,0.0,0.0,
                  0.0,cos(a),-sin(a),0.0,
                  0.0,sin(a),cos(a),0.0,
                  0.0,0.0,0.0,1.0
                );
    }
    
    //绕y轴旋转矩阵
    mat4 rotY(float a) {
      return mat4(cos(a),0.0,sin(a),0.0,
                  0.0,1.0,0.0,0.0,
                  -sin(a),0.0,cos(a),0.0,
                  0.0,0.0,0.0,1.0
                );
    } 

    //平滑交集
    float opSmoothI( float d1, float d2, float k )
    {
      float h = max(k-abs(d1-d2),0.0);
      return max(d1, d2) + h*h*0.25/k;
    }

    //平滑并集
    float opSmoothU( float d1, float d2, float k )
    {
      float h = max(k-abs(d1-d2),0.0);
      return min(d1, d2) - h*h*0.25/k;
    }

    //平滑差集
    float opSmoothS( float d1, float d2, float k )
    {
      float h = max(k-abs(-d1-d2),0.0);
      return max(-d1, d2) + h*h*0.25/k;
    }

    //掏空
    float opOnion(float sdf, in float thickness )
    {
      return abs(sdf)-thickness;
    }

    //并集
    vec2 opU( vec2 d1, vec2 d2 )
    {
      return (d1.x<d2.x) ? d1 : d2;
    }


    //球体
    float sdSphere( vec3 p, float s )
    {
      return length(p)-s;
    }

    //立方体
    float sdBox( vec3 p, vec3 b,float rad )
    {
      vec3 d = abs(p) - b;
      return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0)) - rad;
    }

    //圆柱
    float sdCylinder( vec3 p, float h, float r )
    {
      vec2 d = abs(vec2(length(p.xz),p.y)) - vec2(r,h);
      return min(max(d.x,d.y),0.0) + length(max(d,0.0));
    }

    //甜圈圈
    float sdTorus( vec3 p, vec2 t )
    {
      vec2 q = vec2(length(p.xz)-t.x,p.y);
      return length(q)-t.y;
    }

    vec3 opTwistY( in vec3 p ,float k)
    {

      float c = cos(k*p.y);
      float s = sin(k*p.y);
      mat2  m = mat2(c,-s,s,c);
      p.xz *=m;
      return p;
    }

    vec3 opCheapBend(in vec3 p ,float k)
    {

      float c = cos(k*p.x);
      float s = sin(k*p.x);
      mat2  m = mat2(c,-s,s,c);
      vec3  q = vec3(m*p.xy,p.z);
      return q;
    }


    vec2 twistBox (vec4 pos){

      //绘制原立方体
      pos.x += 2.5; 
      float box1 = sdBox(pos.xyz,vec3(0.4,1.2,0.4),0.06);//方块
      vec2 mBox1 = vec2(box1,1.0);

      //绘制扭曲立方体
      pos.x -= 5.0; 
      pos.xyz=opTwistY(pos.xyz,3.0);//扭曲
      float box2 = sdBox(pos.xyz,vec3(0.4,1.2,0.4),0.06);//方块
      vec2 mBox2 = vec2(box2,1.0);

      vec2 res = opU(mBox1,mBox2);

      return res;
    }

    vec2 onionCylinder (vec4 pos){

      //绘制圆柱
      vec4 pos1 = pos;
      pos1.z -= 6.0; 
      float cylinder1 = sdCylinder(pos1.xyz,0.4,0.6);
      vec2 mCylinder1 = vec2(cylinder1,1.0);

      //绘制裁剪圆柱
      vec4 pos2 = pos;
      pos2.z -= 2.5; 
      float cylinder2 = sdCylinder(pos2.xyz,0.4,0.6);
      cylinder2 = max( cylinder2, pos2.y );//裁剪

      vec2 mCylinder2 = vec2(cylinder2,1.0);

      //绘制裁剪掏空圆柱
      vec4 pos3 = pos;

      float cylinder3 = sdCylinder(pos3.xyz,0.4,0.6);

      cylinder3 = opOnion(cylinder3,0.01);//掏空
      cylinder3 = max( cylinder3, pos3.y );//裁剪
      vec2 mCylinder3 = vec2(cylinder3,1.0);

      vec2 res = opU(mCylinder1,mCylinder2);
      res = opU(res,mCylinder3);

      return res;
    }
    vec2 onionSphere (vec4 pos){

      //绘制球体
      vec4 pos1 = pos;
      pos1.z -= 6.0; 
      float sphere1 = sdSphere(pos1.xyz,0.6);
      vec2 mSphere1 = vec2(sphere1,2.0);

      //绘制裁剪球体
      vec4 pos2 = pos;
      pos2.z -= 2.5; 
      float sphere2 = sdSphere(pos2.xyz,0.6);
      sphere2 = max( sphere2, pos2.y );//裁剪

      vec2 mSphere2 = vec2(sphere2,2.0);

      //绘制裁剪掏空球体
      vec4 pos3 = pos;
     
      float sphere3 = sdSphere(pos3.xyz,0.6);
 
      sphere3 = opOnion(sphere3,0.01);//掏空
      sphere3 = max( sphere3, pos3.y );//裁剪
      vec2 mSphere3 = vec2(sphere3,2.0);

      vec2 res = opU(mSphere1,mSphere2);
      res = opU(res,mSphere3);

      return res;
    }

    vec2 onionTorus (vec4 pos){

      //绘制甜圈圈
      vec4 pos1 = pos;
      pos1.z -= 6.0; 
      float torus1 = sdTorus(pos1.xyz,vec2(0.6,0.2));
      vec2 mTorus1 = vec2(torus1,3.0);

      //绘制裁剪甜圈圈
      vec4 pos2 = pos;
      pos2.z -= 2.5; 
      float torus2 = sdTorus(pos2.xyz,vec2(0.6,0.2));
      torus2 = max( torus2, pos2.y );//裁剪

      vec2 mTorus2 = vec2(torus2,3.0);

      //绘制裁剪掏空甜圈圈
      vec4 pos3 = pos;

      float torus3 = sdTorus(pos3.xyz,vec2(0.6,0.2));

      torus3 = opOnion(torus3,0.01);//掏空
      torus3 = max( torus3, pos3.y );//裁剪
      vec2 mTorus3 = vec2(torus3,3.0);

      vec2 res = opU(mTorus1,mTorus2);
      res = opU(res,mTorus3);

      return res;
    }

    vec2 cheapBendBox (vec4 pos){

      //绘制原立方体面
      pos.x += 2.5; 
      float box1 = sdBox(pos.xyz,vec3(1.0,0.2,0.6),0.06);//方块
      vec2 mBox1 = vec2(box1,1.0);

      //绘制弯曲立方体面
      pos.x -= 5.0; 
      pos.xyz=opCheapBend(pos.xyz,0.5);//弯曲
      float box2 = sdBox(pos.xyz,vec3(1.0,0.2,0.6),0.06);//方块
      vec2 mBox2 = vec2(box2,1.0);

      vec2 res = opU(mBox1,mBox2);

      return res;
    }

    vec2 getDistandMaterial(vec3 p){

      vec2 res =  vec2(p.y,0.0);//地面
 
      vec3 pos = p-vec3(0,2,5);//确定模型的中心
      vec4 oriPos = vec4(pos,1.0);//转为其次坐标
 
      //扭曲
      vec2 twistBox = twistBox(oriPos);

      //弯曲
      vec2 cheapBendBox = cheapBendBox(oriPos);


      //裁剪掏空
      vec2 onionSphere = onionSphere(oriPos);

      vec2 onionCylinder = onionCylinder(oriPos+vec4(2.0,0.0,0.0,1.0));

      vec2 onionTorus = onionTorus(oriPos-vec4(2.0,0.0,0.0,1.0));

   
      //res = opU(res,twistBox);//扭曲立方体
      //res = opU(res,cheapBendBox); //弯曲立方体
      res = opU(res,onionSphere);//裁剪掏空球体
      res = opU(res,onionCylinder);//裁剪掏空圆柱
      res = opU(res,onionTorus);//裁剪掏空甜圈圈

      return res;
    }

    vec2 rayMarch(vec3 rayStart, vec3 rayDirection) {
      float depth=0.;
      float material=0.;
      for(int i=0; i<MAX_STEPS; i++) {
        vec3 p = rayStart + rayDirection*depth;//上一次步进结束后的坐标也就是这一次步进出发点

        vec2 dm = getDistandMaterial(p);
        float dist = dm.x;//获取当前步进出发点与物体相交时距离
        material = dm.y;
        depth += dist; //步进长度累加

        if(depth>MAX_DIST || dist<SURF_DIST) break;//步进距离大于最大步进距离或与物体表面距离小于最小表面距离(光线进入物体)停止前进
      }
      return vec2(depth,material);
    }

    vec3 getNormal(vec3 p){
      return normalize(vec3(
        getDistandMaterial(vec3(p.x + SURF_DIST, p.y, p.z)).x - getDistandMaterial(vec3(p.x - SURF_DIST, p.y, p.z)).x,
        getDistandMaterial(vec3(p.x, p.y + SURF_DIST, p.z)).x - getDistandMaterial(vec3(p.x, p.y - SURF_DIST, p.z)).x,
        getDistandMaterial(vec3(p.x, p.y, p.z  + SURF_DIST)).x - getDistandMaterial(vec3(p.x, p.y, p.z - SURF_DIST)).x
      ));
    }


    //Blinn-Phong模型光照计算
    vec3 calcBlinnPhongLight( vec3 materialColor, vec3 p, vec3 ro) {

      vec3 lightPos = vec3(5.0 * sin(u_time), 20.0, 10.0*cos(u_time)-18.);//光源坐标

      //计算环境光
      float k_a = 0.3;//环境光反射系数
      vec3 ambientLight = 0.6 * vec3(1.0, 1.0, 1.0);
      vec3 ambient = k_a*ambientLight;
      
      vec3 N = getNormal(p); //法线
      vec3 L = normalize(lightPos - p); //光照方向
      vec3 V = normalize(ro - p); //视线
      vec3 H = normalize(V+L); //半程向量

      float r = length(lightPos - p);


      //计算漫反射光
      float k_d = 0.6;//漫反射系数
      float dotLN = clamp(dot(L, N),0.0,1.0);//点乘，并将结果限定在0~1
      vec3 diffuse = k_d * (materialColor/r*r) * dotLN;


      //计算高光反射光
      float k_s = 0.8;//镜面反射系数
      float shininess = 160.0;
      vec3 specularColor = vec3(1.0, 1.0, 1.0);
      vec3 specular = k_s * (specularColor/r*r)* pow(clamp(dot(N, H), 0.0, 1.0), shininess);//计算高光
  
      
      //计算阴影
      vec2 res = rayMarch(p + N*SURF_DIST*2.0,L); 
      if(res.x<length(lightPos-p)-0.001){
        diffuse*=0.1;
      }
        
      //颜色 =  环境光 + 漫反射光 + 镜面反射光
      return ambient +diffuse + specular;
    }

    void main( void ) {

      //窗口坐标调整为[-1,1],坐标原点在屏幕中心
      vec2 st = (gl_FragCoord.xy * 2. - u_resolution) / u_resolution.y;


      vec3 ro = vec3(0.0,4.0,0.0);//视点
      vec3 rd = normalize(vec3(st.x,st.y,1.0));//视线方向

      vec2 res = rayMarch(ro,rd);//反向光线追踪求交点距离与材质ID

      float d = res.x;//物体与视点的距离
      float m = res.y;//材质ID
  
      vec3 p = ro + rd * d;

      vec3 materialColor = vec3(1.0, 0.0, 1.0);//默认材质色，使用差集计算出来的内壁会使用该色填充
      
      
      //为不同物体设置不同的材质颜色
      if(m==0.0){
        materialColor = vec3(.2, 0.0, 0.0);
      }
      if(m==1.0){
        materialColor = vec3(.2, 0.0, 1.0);
      }
      if(m==2.0){
        materialColor = vec3(.7, 0.2, 0.0);
      }
      if(m==3.0){
        materialColor = vec3(.8, .9, 0.0);
      }

      vec3 color = vec3(1.0,1.0,1.0);

      //使用Blinn-Phong模型计算光照
      color *= calcBlinnPhongLight( materialColor, p, ro);

      gl_FragColor = vec4(color, 1.0);

    }
    `

    init();
    animate();

    function init() {
      container = document.getElementById('container');

      camera = new THREE.Camera();
      camera.position.z = 1;

      scene = new THREE.Scene();

      var geometry = new THREE.PlaneBufferGeometry(2, 2);

      uniforms = {
        u_time: {
          type: "f",
          value: 1.0
        },
        u_resolution: {
          type: "v2",
          value: new THREE.Vector2()
        },
        u_mouse: {
          type: "v2",
          value: new THREE.Vector2()
        }
      };

      var material = new THREE.ShaderMaterial({
        uniforms: uniforms,
        vertexShader: vertexShader,
        fragmentShader: fragmentShader
      });

      var mesh = new THREE.Mesh(geometry, material);
      scene.add(mesh);

      renderer = new THREE.WebGLRenderer();
      //renderer.setPixelRatio(window.devicePixelRatio);

      container.appendChild(renderer.domElement);

      onWindowResize();
      window.addEventListener('resize', onWindowResize, false);

      document.onmousemove = function (e) {
        uniforms.u_mouse.value.x = e.pageX
        uniforms.u_mouse.value.y = e.pageY
      }
    }

    function onWindowResize(event) {
      renderer.setSize(800, 800);
      uniforms.u_resolution.value.x = renderer.domElement.width;
      uniforms.u_resolution.value.y = renderer.domElement.height;
    }

    function animate() {
      requestAnimationFrame(animate);
      render();
    }

    function render() {
      uniforms.u_time.value += 0.02;
      renderer.render(scene, camera);
    }
  </script>
</body>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469