smooth isolines

1 files changed, 109 insertions, 170 deletions

diff --git a/Smoke/renderer_gl.c b/Smoke/renderer_gl.c
index 1d87640..59489c5 100644
--- a/Smoke/renderer_gl.c
+++ b/Smoke/renderer_gl.c
@@ -26,9 +26,7 @@
 #include "streamlines.h"
 #include "heightplots.h"
 #include "flowvis.h"
-
 #include "seedpoint.h"
-
 #include "renderer_gl.h"
 
 #define DEFAULT_X_POS 0.0f
@@ -541,59 +539,51 @@ static void render_glyphs(void)
   gluDeleteQuadric(qobj);
 }
 
-//#define percentage(h, l, t) ( 1 - ( (t - min(h, l) ) / (max(h, l) - min(h, l) ) ) )
-
-float percentage(float h, float l, float t)
+struct point get_intersection(struct point pi, struct point pj, float fi, float fj, float t)
 {
-  float perc, delp, deltp;
-  float min = min(h, l);
-  float max = max(h, l);
+  struct point rp;
+  rp.z = 8.0f;
 
-  delp  = max - min;
-  deltp = t - min;
-  perc  = deltp / delp;
+  if (fi == min(fi, fj))
+  {
+    rp.x = (pi.x * (fj - t) + pj.x * (t - fi)) / (fj - fi);
+    rp.y = (pi.y * (fj - t) + pj.y * (t - fi)) / (fj - fi);
+  }
+  else
+  {
+    rp.x = (pj.x * (fi - t) + pi.x * (t - fj)) / (fi - fj);
+    rp.y = (pj.y * (fi - t) + pi.y * (t - fj)) / (fi - fj);
+  }
 
-  return perc;
+  return rp;
 }
 
 static void render_isolines(void)
 {
-  fftw_real wn, hn;
-  int count, DIM;
-  float iso_scale;
+  int idx, count, i, j, state, DIM;
+  float wn, hn, v0, v1, v2, v3, threshold, iso_scale;
+  float x0, y0, z0, x1, y1, z1, x_offset, y_offset;
+  struct color4f color;
+  struct point p0, p1, p2, p3, e1, e2;
 
   DIM = fluids_get_dim();
 
   wn = (fftw_real)winWidth  / (fftw_real)(DIM + 1); // Grid cell width
   hn = (fftw_real)winHeight / (fftw_real)(DIM + 1); // Grid cell height
 
-  iso_scale = isolines_get_threshold_min();
-  iso_scale = isolines_get_threshold_max();
-  iso_scale = isolines_get_nr();
+  z0 = z1 = 8.0f;
 
-  if (isolines_get_nr()) {
-    iso_scale = (float)(fabs(isolines_get_threshold_min() - isolines_get_threshold_max()) / isolines_get_nr());
-  } else {
+  if (isolines_get_nr())
+  {
+    iso_scale = fabs(isolines_get_threshold_min() - isolines_get_threshold_max()) / isolines_get_nr();
+  }
+  else
+  {
     iso_scale = 0.0f;
   }
 
   for (count = 0; count < isolines_get_nr(); count++)
   {
-    int idx;
-    int i, j;
-    float v0, v1, v2, v3;
-    float x0, y0, x1, y1, x_offset, y_offset;
-    int state = 0;
-    static int prev_state = 0;
-    struct color4f color;
-    float threshold;
-
-    static int low = 0;
-    static int hig = 0;
-
-    v0 = v1 = v2 = v3 = 0.0f;
-    x0 = y0 = x1 = y1 = 0.0f;
-    
     threshold = min(isolines_get_threshold_min(), isolines_get_threshold_max()) + count * iso_scale;
 
     glDisable(GL_LIGHTING);
@@ -608,7 +598,6 @@ static void render_isolines(void)
         idx = (j * DIM) + i;
         color = isolines_get_color(get_dataset(idx));
         glColor4f(color.r, color.g, color.b, color.a);
-        //glColor3f(1.0f, 0.0f, 0.0f);
 
         v0 = get_dataset(idx + DIM);
         v1 = get_dataset(idx + 1 + DIM);
@@ -620,145 +609,93 @@ static void render_isolines(void)
         if (v2 >= threshold) { state += 4; }
         if (v3 >= threshold) { state += 8; }
 
-        if (state > 0 && state < 15)
+        x_offset = wn + (fftw_real)i * wn;
+        y_offset = (2 * hn) + (fftw_real)j * hn;
+
+        p0.x = x_offset;
+        p0.y = y_offset + hn;
+        p0.z = z0;
+
+        p1.x = x_offset + wn;
+        p1.y = y_offset + hn;
+        p1.z = z0;
+
+        p2.x = x_offset + wn;
+        p2.y = y_offset;
+        p2.z = z0;
+
+        p3.x = x_offset;
+        p3.y = y_offset;
+        p3.z = z0;
+
+        switch(state)
         {
-          x_offset = wn + (fftw_real)i * wn;
-          y_offset = hn + (fftw_real)j * hn;
-
-          switch (state)
-          {
-            case 1:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-              break;
-
-            case 2:
-              x0 = percentage(v1, v0, threshold) * wn;
-              y0 = hn;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 3:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 4:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 5:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-
-              glVertex3i(x_offset + x0, y_offset + y0, 5.0f);
-              glVertex3i(x_offset + x1, y_offset + y1, 5.0f);
-              
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 6:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-              break;
-
-            case 7:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = 0;
-              y1 = percentage(v3, v0, threshold) * hn;
-              break;
-
-            case 8:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = 0;
-              y1 = percentage(v3, v0, threshold) * hn;
-              break;
-
-            case 9:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-              break;
-
-            case 10:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-
-              glVertex3i(x_offset + x0, y_offset + y0, 5.0f);
-              glVertex3i(x_offset + x1, y_offset + y1, 5.0f);
-              
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 11:
-              x0 = percentage(v3, v2, threshold) * wn;
-              y0 = 0;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 12:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 13:
-              x0 = percentage(v1, v0, threshold) * wn;
-              y0 = hn;
-              x1 = wn;
-              y1 = percentage(v2, v1, threshold) * hn;
-              break;
-
-            case 14:
-              x0 = 0;
-              y0 = percentage(v3, v0, threshold) * hn;
-              x1 = percentage(v1, v0, threshold) * wn;
-              y1 = hn;
-              break;
-
-            case 0: case 15: default:
-              x0 = y0 = x1 = y1 = 0.0f;
-              break;
-          }
-
-           // draw line
-          glVertex3i(x_offset + x0, y_offset + y0, 5.0f);
-          glVertex3i(x_offset + x1, y_offset + y1, 5.0f);
-        }
+          default:
+            state = 0;
+          case 0:
+          case 15:
+            e1.x = e1.y = e1.z = 0.0f;
+            e2.x = e2.y = e2.z = 0.0f;
+            break;
+
+          case 1:
+          case 14:
+            e1 = get_intersection(p0, p3, v0, v3, threshold);
+            e2 = get_intersection(p0, p1, v0, v1, threshold);
+            break;
+
+          case 2:
+          case 13:
+            e1 = get_intersection(p1, p0, v1, v0, threshold);
+            e2 = get_intersection(p1, p2, v1, v2, threshold);
+            break;
+
+          case 3:
+          case 12:
+            e1 = get_intersection(p0, p3, v0, v3, threshold);
+            e2 = get_intersection(p1, p2, v1, v2, threshold);
+            break;
+
+          case 4:
+          case 11:
+            e1 = get_intersection(p2, p3, v2, v3, threshold);
+            e2 = get_intersection(p2, p1, v2, v1, threshold);
+            break;
+
+          case 5:
+          case 10:
+            e1 = get_intersection(p0, p3, v0, v3, threshold);
+            e2 = get_intersection(p0, p1, v0, v1, threshold);
+            glVertex3i(e1.x, e1.y, e1.z);
+            glVertex3i(e2.x, e2.y, e2.z);
+
+            e1 = get_intersection(p2, p3, v2, v3, threshold);
+            e2 = get_intersection(p2, p1, v2, v1, threshold);
+            break;
+
+          case 6:
+          case 9:
+            e1 = get_intersection(p2, p3, v2, v3, threshold);
+            e2 = get_intersection(p1, p0, v1, v0, threshold);
+            break;
+
+          case 7:
+          case 8:
+            e1 = get_intersection(p2, p3, v2, v3, threshold);
+            e2 = get_intersection(p0, p3, v0, v3, threshold);
+            break;
+        } // end switch state
+
+        // draw line
+        glVertex3i(e1.x, e1.y, e1.z);
+        glVertex3i(e2.x, e2.y, e2.z);
       }
-    }
-
+    } // end for count 
     glEnd();
   }
   glEnable(GL_LIGHTING);
 }
 
-
 static void render_streamlines(void)
 {
   int i, j, k, l, idx, DIM;
@@ -785,7 +722,7 @@ static void render_streamlines(void)
       cell_y = round(p.y / cell_y);
 
       idx = cell_x + (cell_y * DIM) + (j * DIM * DIM);
-      //v = frame_history[idx];
+      v = frame_history[idx];
 
       j++;
     }
@@ -994,5 +931,7 @@ void visualize(struct vis_data_arrays *vis_data)
     render_normals();
   }
 
+  render_streamlines();
+
 	render_legend();
 }