diff --git a/src/main.cpp b/src/main.cpp
index d1c0b11..e8649a9 100755
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -124,7 +124,7 @@ void Update() {
     if (frameNumber == 31) {
       isMovePolygon = true;
     } else if (frameNumber == 56) {
-      isMovePolygon = false;
+      // isMovePolygon = false;
       // WaterSolver->polygons[0].omega = 0.f;
     } else if (frameNumber == 401) {
       // isMovePolygon = true;
@@ -292,9 +292,7 @@ int main(int argc, char **argv) {
       // if (frameNumber == 81)
       // if (WaterSolver->particles.size() > 0)
       WaterSolver->Draw();
-
       // SandSolver->Draw();
-
       // ElasticSolver->Draw();
 
       // for testing sdf
@@ -409,10 +407,12 @@ void drawSdf(GLFWwindow *wnd) {
   // for test
   // draw sdf gradient at a given point
   glm::vec2 start = glm::vec2(wx, wy);
-  float dist = SnowSolver->getDistance(start);
-  glm::vec2 grad = SnowSolver->getGradient(start);
+  float dist = WaterSolver->getDistance(start);
+  glm::vec2 grad = WaterSolver->getGradient(start);
   glm::vec2 end = start + dist * grad;
 
+  // std::cout << "dist = " << dist << '\n';
+
   glLineWidth(4);
   glColor3f(1.f, 0.f, 0.f);
 
diff --git a/src/solver.cpp b/src/solver.cpp
index cee0ec4..f188b40 100755
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -23,42 +23,122 @@ Solver<Type>::Solver(const std::vector<Border> &inBorders,
 |			          Functions for computing sdf
 |
 ----------------------------------------------------------------------- */
+// For each grid point, compute its sdf
+// template <typename Type> void Solver<Type>::computeSdfNaive() {
+//   // compute sdf for (sdfWidth * sdfHeight) world space points
+//   // the interval between two adjacent points is sdfCellSize
+//   for (int x = 0; x < sdfWidth; x++) {
+//     for (int y = 0; y < sdfHeight; y++) {
+//       // std::cout << x << ", " << y << '\n';
+//
+//       // world space point
+//       glm::vec2 p = worldOrigin + glm::vec2(sdfCellSize * x, sdfCellSize *
+//       y);
+//
+//       float fDist = 9999.f;
+//       float temp = 0.f;
+//
+//       // node index
+//       int idx = x + y * (sdfWidth - 1); // # of nodes in a row is (sdfWidth -
+//       1)
+//
+//       for (int i = 0; i < polygons.size(); i++) {
+//           temp = (inside_polygon(p, polygons[i])) ? -1.f : 1.f;
+//           temp *= nearest_distance(p, polygons[i]);
+//
+//           if (temp < fDist) {
+//             fDist = temp;
+//
+//             // the nearest polygon from this node
+//             nodes[idx].nearestPolygonPtr = &polygons[i];
+//           }
+//       }
+//
+//       // save sdf distance
+//       nodes[idx].sdfDist = fDist;
+//       // std::cout << "(" << x << ", " << y << ") distance = " << fDist <<
+//       '\n';
+//     }
+//   }
+// }
+
 template <typename Type> void Solver<Type>::computeSdf() {
-  // compute sdf for (sdfWidth * sdfHeight) world space points
-  // the interval between two adjacent points is sdfCellSize
+  // A simple strategy to reduce the cost of computing sdf:
+  // Compute an aabb which is a little bigger than the aabb of the polygon
+  // Compute the grid points (or nodes) occupied by this aabb
+  // Only compute sdf for these grid points
+  // For the orther, set the distance to 9999.f
+
+  // First, we need to initialize the sdf for all grid points
   for (int x = 0; x < sdfWidth; x++) {
     for (int y = 0; y < sdfHeight; y++) {
-      // std::cout << x << ", " << y << '\n';
+      // node index
+      int idx = x + y * (sdfWidth - 1); // # of nodes in a row is (sdfWidth - 1)
+
+      nodes[idx].sdfDist = 9999.f;
+    }
+  }
 
-      // world space point
-      glm::vec2 p = worldOrigin + glm::vec2(sdfCellSize * x, sdfCellSize * y);
+  // Then, iterate all the polygons,
+  // to find the necessary region in which we need to compute sdf
+  for (int i = 0; i < polygons.size(); i++) {
+    Polygon &poly = polygons[i];
 
-      float fDist = 9999.f;
-      float temp = 0.f;
+    // start and end world position of the region
+    glm::vec2 startWorld(poly.lb.x, poly.lb.y);
+    glm::vec2 endWorld(poly.rt.x, poly.rt.y);
 
-      // node index
-      int idx = x + y * (sdfWidth - 1); // # of nodes in a row is (sdfWidth - 1)
+    // world position to grid points
+    glm::ivec2 startGrid = world2Grid(startWorld);
+    glm::ivec2 endGrid = world2Grid(endWorld);
+
+    // extend the region a little bit
+    int alpha = 2;
+    startGrid.x -= alpha;
+    startGrid.y -= alpha;
+    endGrid.x += alpha;
+    endGrid.y += alpha;
+
+    // restrict the region
+    startGrid.x = (startGrid.x < 0) ? 0 : startGrid.x;
+    startGrid.y = (startGrid.y < 0) ? 0 : startGrid.y;
+    endGrid.x = (endGrid.x > X_GRID - 1) ? X_GRID - 1 : endGrid.x;
+    endGrid.y = (endGrid.x > Y_GRID - 1) ? Y_GRID - 1 : endGrid.y;
+
+    // compute the sdf for the necessary region
+    // compute sdf for (sdfWidth * sdfHeight) world space points
+    // the interval between two adjacent points is sdfCellSize
+    for (int x = startGrid.x; x < endGrid.x; x++) {
+      for (int y = startGrid.y; y < endGrid.y; y++) {
+        // std::cout << x << ", " << y << '\n';
 
-      for (int i = 0; i < polygons.size(); i++) {
-        // std::cout << i << '\n';
-        temp = (inside_polygon(p, polygons[i])) ? -1.f : 1.f;
-        temp *= nearest_distance(p, polygons[i]);
-        // std::cout << "here" << '\n';
+        // world space point
+        glm::vec2 p = worldOrigin + glm::vec2(sdfCellSize * x, sdfCellSize * y);
+
+        float fDist = getDistance(p);
+        float temp = 0.f;
+
+        // node index
+        // # of nodes in a row is (sdfWidth - 1)
+        int idx = x + y * (sdfWidth - 1);
+
+        temp = (inside_polygon(p, poly)) ? -1.f : 1.f;
+        temp *= nearest_distance(p, poly);
 
-        // fDist = glm::min(temp, fDist);
         if (temp < fDist) {
           fDist = temp;
 
           // the nearest polygon from this node
-          nodes[idx].nearestPolygonPtr = &polygons[i];
+          nodes[idx].nearestPolygonPtr = &poly;
         }
-      }
 
-      // save sdf distance
-      nodes[idx].sdfDist = fDist;
-      // std::cout << "(" << x << ", " << y << ") distance = " << fDist << '\n';
+        // save sdf distance
+        nodes[idx].sdfDist = fDist;
+        // std::cout << "(" << x << ", " << y << ") distance = " << fDist <<
+        // '\n';
+      }
     }
-  }
+  } // end iterate all polygons
 }
 
 template <typename Type>
@@ -126,17 +206,18 @@ float Solver<Type>::nearest_distance(glm::vec2 p, Polygon &poly) {
 template <typename Type> float Solver<Type>::getDistance(glm::vec2 p) {
 
   // world position to sdf index
-  int idx_x, idx_y;
-  idx_x = int(glm::floor(p.x / sdfCellSize));
-  idx_y = int(glm::floor(p.y / sdfCellSize));
+  glm::ivec2 gridPos = world2Grid(p);
+  // int idx_x, idx_y;
+  // idx_x = int(glm::floor(p.x / sdfCellSize));
+  // idx_y = int(glm::floor(p.y / sdfCellSize));
 
-  if (idx_x < 0 || idx_x > sdfWidth - 1) {
+  if (gridPos.x < 0 || gridPos.x > sdfWidth - 1) {
     return 9999.f;
-  } else if (idx_y < 0.f || idx_y > sdfHeight - 1) {
+  } else if (gridPos.y < 0.f || gridPos.y > sdfHeight - 1) {
     return 9999.f;
   } else {
-    int idx = idx_x + idx_y * (sdfWidth - 1);
-    // std::cout << idx << '\n';
+    int idx = gridPos.x + gridPos.y * (sdfWidth - 1);
+    // std::cout << gridPos << '\n';
     return nodes[idx].sdfDist;
   }
 }
@@ -174,21 +255,24 @@ template <typename Type> glm::vec2 Solver<Type>::getGradient(glm::vec2 p) {
 // p is world position
 template <typename Type> Polygon *Solver<Type>::getPolygon(glm::vec2 p) {
   // world position to sdf index
-  int idx_x, idx_y;
-  idx_x = int(glm::floor(p.x / sdfCellSize));
-  idx_y = int(glm::floor(p.y / sdfCellSize));
+  glm::ivec2 gridPos = world2Grid(p);
+  // int idx_x, idx_y;
+  // idx_x = int(glm::floor(p.x / sdfCellSize));
+  // idx_y = int(glm::floor(p.y / sdfCellSize));
 
-  if (idx_x < 0 || idx_x > sdfWidth - 1) {
+  if (gridPos.x < 0 || gridPos.x > sdfWidth - 1) {
     return NULL;
-  } else if (idx_y < 0.f || idx_y > sdfHeight - 1) {
+  } else if (gridPos.y < 0.f || gridPos.y > sdfHeight - 1) {
     return NULL;
   } else {
-    int idx = idx_x + idx_y * (sdfWidth - 1);
+    int idx = gridPos.x + gridPos.y * (sdfWidth - 1);
     return nodes[idx].nearestPolygonPtr;
   }
 }
 
-template <typename Type> void Solver<Type>::applySdfCollision(Type &p) {
+// if collision happens, return true
+// otherwise, return false
+template <typename Type> bool Solver<Type>::applySdfCollision(Type &p) {
   // particle world space position
   glm::vec2 pos(p.Xp[0], p.Xp[1]);
   glm::vec2 v(p.Vp[0], p.Vp[1]);
@@ -203,10 +287,6 @@ template <typename Type> void Solver<Type>::applySdfCollision(Type &p) {
 
   glm::vec2 vco(0.f, 0.f);
 
-  // 基于 sdf 的碰撞检测
-  // 和边界处的碰撞检测
-  // 唯一不同的是法向量 n 的计算方式
-  // 而碰撞后的速度，采用同样��处理
   // for sdf collision detection
   // narrow band threshold
   if (glm::abs(dist) < NARROW_BAND) {
@@ -293,6 +373,8 @@ template <typename Type> void Solver<Type>::applySdfCollision(Type &p) {
       //           << '\n';
     }
   }
+
+  return isCollisionOn;
 }
 
 /* -----------------------------------------------------------------------
diff --git a/src/solver.h b/src/solver.h
index ebc038e..5e625de 100755
--- a/src/solver.h
+++ b/src/solver.h
@@ -47,12 +47,21 @@ template <typename Type> class Solver {
   float getDistance(glm::vec2);
   glm::vec2 getGradient(glm::vec2);
   Polygon *getPolygon(glm::vec2);
-  void applySdfCollision(Type &);
+  bool applySdfCollision(Type &);
 
   void Draw(); // Draw particles, border and nodes (if selected)
   void WriteToFile(
       int frame); // Write point cloud coordinates to .ply file (Houdini)
 
+  static glm::ivec2 world2Grid(glm::vec2 worldPos) {
+    glm::ivec2 gridPos;
+
+    gridPos.x = int(glm::floor(worldPos.x / sdfCellSize));
+    gridPos.y = int(glm::floor(worldPos.y / sdfCellSize));
+
+    return gridPos;
+  }
+
 /* Static functions */
 #if INTERPOLATION == 1
   static float Bspline(float x) // Cubic Bspline
diff --git a/video/output.gif b/video/output.gif
index 0ef291b..ef55e07 100644
Binary files a/video/output.gif and b/video/output.gif differ