Use rendering at startup only for perlin ground

But when the ground is not textured, this looks just as if the ball was
not moving… :(

Ball.cpp

@@ -5,10 +5,11 @@
 bool bouncing = true;
 double stop_bouncing;
 
-Ball::Ball(const Point& _center, double _min_height, double _v_x, double _v_z,
-        double _p, double _q) :
+Ball::Ball(const Point& _center, const Ground* _ground, double _min_height,
+        double _v_x, double _v_z, double _p, double _q) :
     Center(_center),
     surface(_center, _min_height, _p, _q),
+    ground(_ground),
     radius(_p),
     init_h(_center.y),
     min_height(_min_height),
@@ -24,22 +25,25 @@ Ball::Ball(const Point& _center, double _min_height, double _v_x, double _v_z,
 }
 
 void Ball::_compute_pos(double dt) {
+    double height = (*ground)(Center.x, Center.z);
     if (bouncing) {
         Center.y = fmax(
-            min_height,
-            A + B * crt_time - (G_CTE / 2) * crt_time * crt_time + min_height
+            min_height + height,
+            A + B * crt_time - (G_CTE / 2) * crt_time * crt_time + min_height 
         );
     } else {
         double n_time = crt_time - stop_bouncing;
         double min_rad = radius - min_height;
-        Center.y = min_height + ((min_rad * (1.0 - exp(-n_time)) + min_rad )/2.) +
+        v_x *= 0.8;
+        v_z *= 0.8;
+        Center.y = height + min_height + ((min_rad * (1.0 - exp(-n_time)) + min_rad )/2.) +
             (min_rad - ((min_rad* (1- exp(-n_time)) + min_rad)/2.)) * sin(5. *
                     n_time);
     }
     Center.x += dt * v_x;
     Center.z += dt * v_z;
     surface.update_center_pos(Center);
-    surface.check_horizontal_plan_collision(0.0);
+    surface.check_ground_collision(ground);
 }
 
 void Ball::_compute_T_n() {
@@ -59,21 +63,23 @@ void Ball::_compute_U_n() {
 }
 
 void Ball::_compute_v_x(Point normal) {
-    double factor = normal.x / (4*(normal.x + normal.z));
+    double factor = 0;
+    if (normal.x + normal.z != 0) {
+        factor = normal.x / (4*(normal.x + normal.z));
+    }
     v_x *= (0.5 + factor);
 }
 
 void Ball::_compute_v_z(Point normal) {
-    double factor = normal.z / (4*(normal.x + normal.z));
+    double factor = 0;
+    if (normal.x + normal.z != 0) {
+        factor = normal.z / (4*(normal.x + normal.z));
+    }
     v_z *= (0.5 + factor);
 }
 
-void Ball::update_pos(double dt) {
-    crt_time += dt;
-    if ((bouncing) && (crt_time >= T)) {
-        double old_t = T;
-        double max_t;
-        Point normal = surface.getNormalVector();
+void Ball::_compute_parameters() {
+        Point normal = (*ground).get_normal(Center.x, Center.z);
         bounce_number += 1;
         _compute_U_n();
         _compute_A_n();
@@ -81,6 +87,16 @@ void Ball::update_pos(double dt) {
         _compute_T_n();
         _compute_v_x(normal);
         _compute_v_z(normal);
+        min_height = fmin(radius, min_height + 0.2 * (radius - min_height));
+}
+
+void Ball::update_pos(double dt) {
+    //double height = (*ground)(Center.x, Center.z);
+    crt_time += dt;
+    if ((bouncing) && (crt_time >= T)) {
+        double old_t = T;
+        double max_t;
+        _compute_parameters();
         max_t = (T - old_t)/2.0 + old_t;
 
         if (((A + B * max_t - (G_CTE / 2) * max_t * max_t + min_height) < radius)) {

Ball.hpp

@@ -12,6 +12,8 @@
 #pragma once
 #include <cstddef>
 #include "spheroid.hpp"
+#include "FlatGround.hpp"
+#include "PerlinGround.hpp"
 
 #define G_CTE 9.81
 
@@ -19,6 +21,7 @@ class Ball {
     private:
         Point Center;
         Spheroid surface;
+        const Ground* ground;
         double radius;
         double init_h;
         double min_height;
@@ -27,6 +30,7 @@ class Ball {
         double A, B, U, T; // Coefficients for the physical model.
         double v_x, v_z;
         void _compute_pos(double dt);
+        void _compute_parameters();
         void _compute_v_x(Point normal);
         void _compute_v_z(Point normal);
         void _compute_A_n();
@@ -34,7 +38,7 @@ class Ball {
         void _compute_U_n();
         void _compute_T_n();
     public:
-        Ball(const Point& _center, double _min_height, double _v_x,
+        Ball(const Point& _center, const Ground* ground, double _min_height, double _v_x,
                 double _v_z, double p, double q);
         void update_pos(double dt);
         Point getCenter() const {return Center;}

FlatGround.cpp

@@ -0,0 +1,10 @@
+#include "FlatGround.hpp"
+
+
+double FlatGround::operator() (double, double) const {
+    return 0. ;
+}
+
+Point FlatGround::get_normal(double, double) const {
+    return Point(0, 1, 0);
+}

FlatGround.hpp

@@ -0,0 +1,8 @@
+#pragma once
+#include "Ground.hpp"
+
+class FlatGround : public Ground {
+    public:
+        double operator() (double, double) const;
+        Point get_normal(double x, double y) const;
+};

Ground.hpp

@@ -0,0 +1,9 @@
+#pragma once
+
+#include "common_structures.hpp"
+
+class Ground {
+    public:
+        virtual double operator() (double, double) const = 0;
+        virtual Point get_normal(double x, double y) const = 0;
+};

GroundFlatMesh.cpp

@@ -0,0 +1,45 @@
+#include "GroundFlatMesh.hpp"
+
+#include <vector>
+
+const int GroundFlatMesh::MIN_I = -40;
+const int GroundFlatMesh::MAX_I = 100;
+const int GroundFlatMesh::MIN_J = -40;
+const int GroundFlatMesh::MAX_J = 100;
+
+GroundFlatMesh::GroundFlatMesh(const Point& center, double decay_speed)
+    : center(center), decay_speed(decay_speed)
+{
+    std::vector<std::vector<size_t> > vertice_at(MAX_I - MIN_I + 1,
+            std::vector<size_t>(MAX_J - MIN_J + 1));
+
+    for(int i=MIN_I; i < MAX_I + 1; ++i) {
+        for(int j=MIN_J; j < MAX_J + 1; ++j) {
+            vertice_at[i - MIN_I][j - MIN_J] =
+                output.add_vertice(
+                        tile_position(i, j),
+                        Point(0., 1., 0.));
+        }
+    }
+
+    for(int i=MIN_I; i < MAX_I; ++i) {
+        for(int j=MIN_J; j < MAX_J; ++j) {
+            output.add_face(Face(
+                        vertice_at[i - MIN_I][j - MIN_J],
+                        vertice_at[i - MIN_I][j + 1 - MIN_J],
+                        vertice_at[i + 1 - MIN_I][j + 1 - MIN_J]));
+            output.add_face(Face(
+                        vertice_at[i - MIN_I][j - MIN_J],
+                        vertice_at[i + 1 - MIN_I][j + 1 - MIN_J],
+                        vertice_at[i + 1- MIN_I][j - MIN_J]));
+        }
+    }
+}
+
+double GroundFlatMesh::ith_dist(int i) const {
+    return ((i < 0) ? -1 : 1) * decay_speed * i * i;
+}
+
+Point GroundFlatMesh::tile_position(int i, int j) const {
+    return Point(ith_dist(i), 0., ith_dist(j));
+}

GroundFlatMesh.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include "Mesh.hpp"
+
+class GroundFlatMesh {
+    public:
+        GroundFlatMesh(const Point& center, double decay_speed);
+
+        Mesh* get_mesh() { return &output; }
+
+    private: //meth
+        double ith_dist(int i) const;
+        Point tile_position(int i, int j) const;
+
+    private:
+        static const int MIN_I, MAX_I, MIN_J, MAX_J;
+        Mesh output;
+        const Point& center;
+        double decay_speed;
+};

Implicit.hpp

@@ -8,9 +8,16 @@ class ImplicitSurface {
         double operator()(const Point& pt) const;
         Point getCenter() const { return center;}
 
+        const Color& get_color() const { return color; }
+        void set_color(const Color& _color) { color = _color; }
+
+        virtual Point location_hint() const = 0;
+
         /// Compute the normal vector to the isosurface at `pt`
         Point normal_at(const Point& pt) const;
     protected:
         Point center;
+        Color color;
+
         ImplicitSurface(Point _center) : center(_center) {}
 };

Makefile

@@ -11,11 +11,15 @@ OBJS=Implicit.o \
 	 Mesh.o \
 	 spheroid.o \
 	 Ball.o \
+	 Ground.o \
+	 FlatGround.o \
+	 PerlinGround.o \
 	 PerlinNoise.o \
 	 util/ObjParser.o \
 	 MarchingCubes.o \
 	 _gen/marching_cubes_data.o \
 	 periodic_updates.o \
+	 GroundFlatMesh.o \
 	 tests/TestImplicitSphere.o \
 	 render/GlutRender.o
 

MarchingCubes.cpp

@@ -1,6 +1,9 @@
 #include "MarchingCubes.hpp"
 
 #include <cassert>
+#include <queue>
+
+#include "GL/gl.h"
 
 #include "GL/gl.h"
 
@@ -13,14 +16,27 @@ MarchingCubes::MarchingCubes(
     resolution(resolution)
 {}
 
-void MarchingCubes::add_hint(const Cuboid& hint) {
+MarchingCubes& MarchingCubes::add_hint(const Point& hint) {
     hints.push_back(hint);
+    return *this;
 }
 
 Mesh MarchingCubes::operator()() {
     Mesh output;
 
-    without_hints(output);
+    output.set_color(surface.get_color());
+
+    perf_counter = 0;
+
+    if(hints.empty())
+        without_hints(output);
+    else
+        with_hints(output);
+
+#ifdef MC_SHOW_PERF
+    fprintf(stderr, "Explored cubes: %ld\n", perf_counter);
+#endif // MC_SHOW_PERF
+
     output.translate(surface.getCenter());
 
     return output;
@@ -34,11 +50,123 @@ void MarchingCubes::without_hints(Mesh& output) {
             for(double z = box.low(2); z < box.high(2) + resolution;
                     z += resolution) {
                 march_at(x, y, z, output);
+                perf_counter++;
             }
         }
     }
 }
 
+void MarchingCubes::with_hints(Mesh& output) {
+    std::set<PointLoc> explored_cubes;
+    std::queue<PointLoc> process;
+
+    for(const Point& hint: hints) {
+        PointLoc coords = coords_of(hint);
+        PointLoc start_point =
+            seek_closest_intersection(coords, explored_cubes);
+        explored_cubes.erase(start_point);
+        process.push(start_point);
+    }
+
+    while(!process.empty()) {
+        PointLoc pos = process.front();
+        process.pop();
+        if(explored_cubes.find(pos) != explored_cubes.end())
+            continue;
+        explored_cubes.insert(pos);
+
+        Point pt = point_at_coords(pos);
+        Intersections inters =
+            mk_intersection_cube(pt.x, pt.y, pt.z, resolution);
+        if(!inters.has_inters())
+            continue;
+
+        march_at(pt.x, pt.y, pt.z, output);
+
+        for(int dx=-1; dx <= 1; ++dx) {
+            for(int dy=-1; dy <= 1; ++dy) {
+                for(int dz=-1; dz <= 1; ++dz) {
+                    if(dx == 0 && dy == 0 && dz == 0)
+                        continue;
+                    PointLoc nPos = pos + PointLoc(dx, dy, dz);
+                    if(coords_in_box(nPos))
+                        process.push(nPos);
+                }
+            }
+        }
+    }
+    perf_counter = explored_cubes.size();
+}
+
+PointLoc MarchingCubes::seek_closest_intersection(
+        const PointLoc& loc, std::set<PointLoc>& visited) const {
+    std::queue<PointLoc> process;
+    process.push(loc);
+
+    while(!process.empty()) {
+        PointLoc pos = process.front();
+        process.pop();
+        if(visited.find(pos) != visited.end())
+            continue;
+        visited.insert(pos);
+
+        Point pt = point_at_coords(pos);
+        Intersections inters =
+            mk_intersection_cube(pt.x, pt.y, pt.z, resolution);
+        if(inters.has_inters())
+            return pos;
+
+        for(int dx=-1; dx <= 1; ++dx) {
+            for(int dy=-1; dy <= 1; ++dy) {
+                for(int dz=-1; dz <= 1; ++dz) {
+                    if(dx == 0 && dy == 0 && dz == 0)
+                        continue;
+                    PointLoc nPos = pos + PointLoc(dx, dy, dz);
+                    if(coords_in_box(nPos))
+                        process.push(nPos);
+                }
+            }
+        }
+    }
+
+    throw SurfaceNotFound();
+}
+
+Point MarchingCubes::align_on_bb(const Point& pt) const {
+    auto align = [resolution = resolution](double v) {
+        return resolution * floor(v / resolution);
+    };
+    Point dist = pt - box.low_pt();
+    dist.x = align(dist.x);
+    dist.y = align(dist.y);
+    dist.z = align(dist.z);
+    return box.low_pt() + dist;
+}
+
+bool MarchingCubes::coords_in_box(const PointLoc& pt) const {
+    Point low = point_at_coords(pt),
+          high = point_at_coords(pt + PointLoc(1, 1, 1));
+    return
+        low.x >= box.low(0) && low.y >= box.low(1) && low.z >= box.low(2)
+        && high.x <= box.high(0) && high.y <= box.high(1)
+        && high.z <= box.high(2);
+}
+
+PointLoc MarchingCubes::coords_of(const Point& pt) const {
+    Point dist = pt - box.low_pt();
+    return PointLoc(
+            dist.x / resolution,
+            dist.y / resolution,
+            dist.z / resolution);
+}
+
+Point MarchingCubes::point_at_coords(const PointLoc& coords) const {
+    return box.low_pt()
+        + Point(coords.x * resolution,
+                coords.y * resolution,
+                coords.z * resolution);
+}
+
 Point MarchingCubes::CubeEdge::at(double pos,
         double bx, double by, double bz, double resolution) const
 {

MarchingCubes.hpp

@@ -18,6 +18,9 @@ class MarchingCubes {
     private:
         typedef std::set<PointLoc> LocSet;
 
+        class SurfaceNotFound : public std::exception {
+        };
+
         class Intersections {
             public:
                 typedef unsigned char intersect_t;
@@ -34,6 +37,10 @@ class MarchingCubes {
                  * cube of side 1 placed at (0, 0, 0).
                  */
 
+                bool has_inters() const {
+                    return inters != 0 && inters != 0xff;
+                }
+
                 void set_corner(bool x, bool y, bool z, bool val) {
                     intersect_t mask = 1 << (shift(x, y, z));
                     if(val)
@@ -63,15 +70,12 @@ class MarchingCubes {
 
         /** Add a starting point hint
          *
-         * A hint is a cuboid that should intersect at least once the surface,
-         * such that the marching cube will find the surface there and will be
-         * able to follow it.
-         * If at least a hint is given, the algorithm will expect that at least
-         * a hint per disjoint surface is given, ie. that it is safe to only
-         * follow the surface starting from the hints, and ignoring the parts
-         * of the grid that are "far" from the hints.
+         * A hint is a location that should be close to the surface. From each
+         * hint, a BFS-like search will be performed, and will stop when a
+         * surface is first encountered. From there, the MC algorithm will
+         * perform, until the whole surface is explored.
          */
-        void add_hint(const Cuboid& hint);
+        MarchingCubes& add_hint(const Point& hint);
 
         Mesh operator()();
 
@@ -120,6 +124,20 @@ class MarchingCubes {
         /// triangle was added.
 
         void without_hints(Mesh& output);
+        void with_hints(Mesh& output);
+
+        PointLoc seek_closest_intersection(
+                const PointLoc& loc,
+                std::set<PointLoc>& visited) const;
+
+        /** Returns a point close to `pt`, which distance to the bounding box
+         * lower corner are integer multiples of `resolution`. */
+        Point align_on_bb(const Point& pt) const;
+
+        bool coords_in_box(const PointLoc& pt) const;
+
+        PointLoc coords_of(const Point& pt) const;
+        Point point_at_coords(const PointLoc& coords) const;
 
         Intersections mk_intersection_cube(double bx, double by, double bz,
                 double side_len) const;
@@ -134,5 +152,7 @@ class MarchingCubes {
         Cuboid box;
         double resolution;
 
-        std::vector<Cuboid> hints;
+        std::vector<Point> hints;
+
+        size_t perf_counter;
 };

Mesh.hpp

@@ -41,6 +41,9 @@ class Mesh {
         /// Gets the various faces
         const std::vector<Face>& get_faces() const;
 
+        const Color& get_color() const { return color; }
+        void set_color(const Color& _color) { color = _color; }
+
     private: //struct
         struct NormalVect {
             NormalVect(const Point& vect, bool manual=false)
@@ -63,4 +66,5 @@ class Mesh {
         std::vector<std::vector<size_t> > faces_with_vert;
 
         Point center;
+        Color color;
 };

PerlinGround.cpp

@@ -0,0 +1,14 @@
+#include "PerlinGround.hpp"
+
+PerlinGround::PerlinGround() : surface(PerlinNoise()) {}
+PerlinGround::PerlinGround(unsigned int seed) : surface(PerlinNoise(seed)) {}
+
+double PerlinGround::operator() (double x, double z) const {
+    return surface.noise(x, z);
+}
+
+
+Point PerlinGround::get_normal(double x, double z) const {
+    const Point pt(x, surface.noise(x, z), z);
+    return surface.normal_at(pt);
+}

PerlinGround.hpp

@@ -0,0 +1,14 @@
+#pragma once
+#include "Ground.hpp"
+#include "PerlinNoise.hpp"
+
+class PerlinGround : public Ground {
+    public:
+        PerlinGround();
+        PerlinGround(unsigned int seed);
+        double operator() (double, double) const;
+        Point get_normal(double x, double y) const;
+        PerlinNoise* get_surface() { return &surface;}
+    private:
+        PerlinNoise surface;
+};

PerlinNoise.cpp

@@ -44,7 +44,7 @@ PerlinNoise::PerlinNoise(unsigned int seed) : ImplicitSurface(Point(0,0,0)) {
 }
 
 double PerlinNoise::operator() (double x, double y, double z) const {
-    return z - noise(x, y, 0);
+    return (y - fBm(x, 0, z, 2, 0.3, 0.9));
 }
 
 double PerlinNoise::noise(double _x, double _y, double _z) const {
@@ -60,31 +60,31 @@ double PerlinNoise::noise(double _x, double _y, double _z) const {
     double v = fade(_y);
     double w = fade(_z);
 
-    int A = p[X] + Y;
-    int AA = p[A] + Z;
-    int AB = p[A + 1] + Z;
-    int B = p[X + 1] + Y;
-    int BA = p[B] + Z;
-    int BB = p[B + 1] + Z;
+    int A = p[X] + Z;
+    int AA = p[A] + Y;
+    int AB = p[A + 1] + Y;
+    int B = p[X + 1] + Z;
+    int BA = p[B] + Y;
+    int BB = p[B + 1] + Y;
 
     double res = lerp(
        w,
        lerp(
            v,
-           lerp(u, grad(p[AA], _x, _y, _z), grad(p[BA], _x-1, _y, _z)),
-           lerp(u, grad(p[AB], _x, _y-1, _z), grad(p[BB], _x-1, _y-1, _z))
+           lerp(u, grad(p[AA], _x, _z, _y), grad(p[BA], _x-1, _z, _y)),
+           lerp(u, grad(p[AB], _x, _z-1, _y), grad(p[BB], _x-1, _z-1, _y))
        ),
        lerp(
            v,
            lerp(
                u,
-               grad(p[AA+1], _x, _y, _z-1),
-               grad(p[BA+1], _x-1, _y, _z-1)
+               grad(p[AA+1], _x, _z, _y-1),
+               grad(p[BA+1], _x-1, _z, _y-1)
            ),
            lerp(
                u,
-               grad(p[AB+1], _x, _y-1, _z-1),
-               grad(p[BB+1], _x-1, _y-1, _z-1)
+               grad(p[AB+1], _x, _z-1, _y-1),
+               grad(p[BB+1], _x-1, _z-1, _y-1)
            )
        )
     );
@@ -92,7 +92,7 @@ double PerlinNoise::noise(double _x, double _y, double _z) const {
 }
 
 double PerlinNoise::noise(double _x, double _y) const {
-    return noise(_x, _y, 0);
+    return noise(_x, 0, _y);
 }
 
 double PerlinNoise::fade(double t) const {
@@ -111,6 +111,21 @@ double PerlinNoise::grad(int hash, double x, double y, double z)  const {
 	return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
 }
 
+double PerlinNoise::fBm(double x, double y, double z, int octaves, double
+        lacunarity, double gain) const {
+    double amplitude = 1.0;
+    double frequency = 1.0;
+    double sum = 0.0;
+    for(int i = 0; i < octaves; ++i)
+    {
+//        sum += amplitude * (1-abs(noise(x * frequency, y * frequency, z * frequency)));
+        sum += amplitude * noise(x * frequency, y * frequency, z * frequency);
+        amplitude *= gain;
+        frequency *= lacunarity;
+    }
+    return sum;
+}
+
 Point PerlinNoise::normal_vector(double x, double y) {
     double n_x_y = noise(x, y, 0);
     double n_dx_y = noise(x + 0.01, y, 0);
@@ -121,3 +136,11 @@ Point PerlinNoise::normal_vector(double x, double y) {
     );
     return product;
 }
+
+Point PerlinNoise::location_hint() const {
+    return location_hint(0, 0);
+}
+
+Point PerlinNoise::location_hint(double x, double z) const {
+    return Point(x, noise(x, z), z);
+}

PerlinNoise.hpp

@@ -1,6 +1,7 @@
 /**
  * Perlin Noise implementation for the ground (header file)
  **/
+#pragma once
 #include "Implicit.hpp"
 #include <cmath>
 #include <random>
@@ -15,6 +16,11 @@ class PerlinNoise : public ImplicitSurface {
         double operator() (double _x, double _y, double _z) const;
         double noise(double x, double y) const;
         Point normal_vector(double x, double y);
+        double fBm(double x, double y, double z, int octaves, double
+                lacunarity, double gain) const;
+
+        virtual Point location_hint() const;
+        Point location_hint(double x, double z) const;
     private:
         double noise(double x, double y, double z) const;
         double fade(double t) const;

README.md

@@ -1,4 +1,40 @@
 # mpri-graphics-project
 
-Un repo pour le projet de graphics du MPRI. On va faire rebondir des trucs. Ça
-va être distrayant. On va vraiment s'amuser. Et valider un cours. Peut-être.
+
+A repository for the MPRI's Computer Graphics and Visualization course project.
+A jelly-like ball bounces around.
+
+[Course page](https://wikimpri.dptinfo.ens-cachan.fr/doku.php?id=cours:c-2-39)
+|
+[Slides](https://tobast.fr/m2/graphics/jelly.pdf)
+
+## Compiling
+
+The basic compilation is achieved by a simple `make`.
+
+You can add additional compilation flags with `ADD_FLAGS="..." make`.
+
+The produced binary is `./bounce.bin`
+
+## Compilation flags
+
+* `-DDEBUG_DISPLAY_WIREFRAME`: display the wireframe of meshes
+* `-DDEBUG_DISPLAY_NORMAL`: display the computed normal vectors for each
+  vertice
+* `-DMC_SHOW_PERF`: display Marching Cubes performance stats
+
+## Run flags
+
+* `-perlin`: replace flat ground by perlin noise
+* `-qwerty`, `-azerty`: change the used keymap
+
+For instance, `./bounce.bin -perlin -qwerty` will use perlin floor and a qwerty
+keymap.
+
+## In-app commands
+
+* space: play/pause
+* <, >: slow down/speed up animation
+* 0: reset to default speed
+* w, a, s, d (z, q, s, d in azerty): move camera
+* q, e (a, e in azerty): rotate camera around y axis

common_structures.cpp

@@ -23,6 +23,14 @@ double Cuboid::high(unsigned dim) const {
     return highBound[dim];
 }
 
+Point Cuboid::low_pt() const {
+    return lowBound;
+}
+
+Point Cuboid::high_pt() const {
+    return highBound;
+}
+
 double Cuboid::length(unsigned dim) const {
     assert(dim < 3);
     return high(dim) - low(dim);

common_structures.hpp

@@ -143,12 +143,29 @@ class Cuboid {
         double high(unsigned dim) const; ///< Higher bound for a dimension
         double length(unsigned dim) const;
 
+        Point low_pt() const;
+        Point high_pt() const;
+
         double volume() const;
         bool is_empty() const;
     private:
         Point lowBound, highBound;
 };
 
+struct Color {
+    Color() : r(0), g(0), b(0) {}
+    Color(double r, double g, double b) : r(r), g(g), b(b) {}
+    double r, g, b;
+    double operator[](int i) const {
+        switch(i % 3) {
+            case 0: return r;
+            case 1: return g;
+            case 2: return b;
+        }
+        return 0; // won't happen
+    }
+};
+
 namespace std {
 	template<> struct hash<Face>
     {

main_ball.cpp

@@ -5,13 +5,16 @@
 #include "Ball.hpp"
 #include <cstdio>
 #include <iostream>
+#include "FlatGround.hpp"
 
 using namespace std;
 
 int main() {
     int i;
     Point center(0,10,0);
-    Ball ball(center, 0, 0.25, 0,  1, 1);
+    FlatGround* flat = new FlatGround();
+
+    Ball ball(center,flat, 0, 0.25, 0,  1, 1);
     for(i=0; i< 10000; i++) {
         ball.update_pos(0.001);
         cout << ball << "\n";

main_bounce.cpp

@@ -2,26 +2,66 @@
  * bouncing implicit-surface defined sphere.
  **/
 
+#include <cstring>
+
 #include "render/GlutRender.hpp"
 #include "Ball.hpp"
+#include "FlatGround.hpp"
 #include "MarchingCubes.hpp"
+#include "GroundFlatMesh.hpp"
 
 #include "periodic_updates.hpp"
 
 int main(int argc, char** argv) {
+    bool perlin = false;
+    for(int pos=1; pos < argc; ++pos) {
+        if(strcmp(argv[pos], "-perlin") == 0)
+            perlin = true;
+        if(strcmp(argv[pos], "-azerty") == 0)
+            set_key_mapping(KeyMappings::azerty());
+        if(strcmp(argv[pos], "-qwerty") == 0)
+            set_key_mapping(KeyMappings::qwerty());
+    }
+
+    // Last minute switch, this code is ugly, please close your eyes until
+    // stated otherwise.
+    PerlinGround perlin_ground;
+    FlatGround flat_ground;
+    GroundFlatMesh ground_mesh(Point(0., 0., 0.), 0.05);
+    Ground* ball_ground = nullptr;
+    if(perlin)
+        ball_ground = &perlin_ground;
+    else
+        ball_ground = &flat_ground;
+    // You can open your eyes, now.
+
     GlutRender& render = GlutRender::get_instance();
     render.init(&argc, argv, 640, 480, "Bouncing stuff");
 
-    Ball ball(Point(0, 5, 0), 0.75, -.5, -.7, 1, 1);
+    Ball ball(Point(0, 5, 0), ball_ground, 0.55, -.5, -.7, 1, 1);
+    ball.get_surface()->set_color(Color(1., 0., 0.));
 
     Cuboid bbox = ball.get_surface()->max_bounding_box();
+    Cuboid bbox_2(Point(-20, -3, -20), Point(20,3,20));
     printf("%.2lf %.2lf %.2lf | %.2lf %.2lf %.2lf\n",
             bbox.low(0), bbox.low(1), bbox.low(2),
             bbox.high(0), bbox.high(1), bbox.high(2));
     render.add_surface(ball.get_surface(), bbox);
 
+    if(perlin) {
+        perlin_ground.get_surface()->set_color(Color(0.13, 0.82, 0.21));
+        render.add_surface(perlin_ground.get_surface(), bbox_2, false, 0.50);
+    }
+    else {
+        ground_mesh.get_mesh()->set_color(Color(0.13, 0.82, 0.21));
+        render.add_mesh(ground_mesh.get_mesh());
+    }
+    //render.follow_implicit_position(ball.get_surface());
+
     render.set_idle_func(periodic_update);
-    init_periodic_static(&ball);
+    render.add_kb_handler(periodic_kb_handler);
+    render.add_kb_up_handler(periodic_kb_up_handler);
+    init_periodic_static(&ball, &render);
     render.run();
 
     return 0;

periodic_updates.cpp

@@ -1,17 +1,66 @@
 #include "periodic_updates.hpp"
 
 #include <ctime>
+#include <cmath>
 #include <chrono>
 #include <GL/glut.h>
 
+struct Movement {
+    Movement()  :
+        fwd(false),
+        bck(false),
+        left(false),
+        right(false),
+        turn_l(false),
+        turn_r(false),
+        sight_angle(3.14159)
+    {}
+
+    bool
+        fwd,
+        bck,
+        left,
+        right,
+        turn_l,
+        turn_r;
+
+    double sight_angle;
+
+    void tick() {
+        sight_angle += (turn_l - turn_r) * .05;
+    }
+
+    Point movement() {
+        Point front_dir = Point(-sin(sight_angle), 0, -cos(sight_angle)),
+              left_dir = Point(-cos(sight_angle), 0, sin(sight_angle));
+        return (fwd - bck) * front_dir
+            + (left - right) * left_dir;
+    }
+    Point sight() {
+        return Point(-sin(sight_angle), 0., -cos(sight_angle));
+    }
+};
+
 static Ball* _ball = nullptr;
 static std::chrono::time_point<std::chrono::steady_clock>
     _last_time, _init_clocks;
 
-void init_periodic_static(Ball* ball) {
+static bool is_paused = false;
+static double speed_factor = 1.;
+static KeyMappings _keymap = KeyMappings::qwerty();
+static Movement _movement;
+static Point _position(0., 2.5, -10.);
+static GlutRender* _render = nullptr;
+
+void init_periodic_static(Ball* ball, GlutRender* render) {
     _last_time = std::chrono::steady_clock::now();
     _init_clocks = std::chrono::steady_clock::now();
     _ball = ball;
+    _render = render;
+}
+
+void set_key_mapping(const KeyMappings& mapping) {
+    _keymap = mapping;
 }
 
 double ellapsed_double(
@@ -24,8 +73,52 @@ double ellapsed_double(
 
 void periodic_update() {
     auto now = std::chrono::steady_clock::now();
-    _ball->update_pos(ellapsed_double(_last_time, now));
+
+    if(!is_paused)
+        _ball->update_pos(speed_factor * ellapsed_double(_last_time, now));
+
+    _movement.tick();
+    _position += _movement.movement();
+    Point look_at = _position + _movement.sight();
+
+    _render->set_camera(_position, look_at);
 
     _last_time = now;
     glutPostRedisplay();
 }
+
+void mvt_keys_update(bool up, unsigned key) {
+    if(key == _keymap.mv_fwd)
+        _movement.fwd = !up;
+    else if(key == _keymap.mv_bck)
+        _movement.bck = !up;
+    else if(key == _keymap.mv_left)
+        _movement.left = !up;
+    else if(key == _keymap.mv_right)
+        _movement.right = !up;
+    else if(key == _keymap.turn_left)
+        _movement.turn_l = !up;
+    else if(key == _keymap.turn_right)
+        _movement.turn_r = !up;
+}
+
+void periodic_kb_handler(unsigned char key, int, int) {
+    if(key == ' ')
+        is_paused = !is_paused;
+    else if(key == '<') {
+        speed_factor -= .1;
+        if(speed_factor <= 0.05)
+            speed_factor = .1;
+    }
+    else if(key == '>') {
+        speed_factor += .1;
+    }
+    else if(key == '0') {
+        speed_factor = 1.;
+    }
+    mvt_keys_update(false, key);
+}
+
+void periodic_kb_up_handler(unsigned char key, int, int) {
+    mvt_keys_update(true, key);
+}

periodic_updates.hpp

@@ -1,6 +1,45 @@
 #pragma once
 
 #include "Ball.hpp"
+#include "render/GlutRender.hpp"
 
-void init_periodic_static(Ball* ball);
+struct KeyMappings {
+    static KeyMappings qwerty() {
+        return KeyMappings('w', 's', 'a', 'd', 'q', 'e');
+    }
+    static KeyMappings azerty() {
+        return KeyMappings('z', 's', 'q', 'd', 'a', 'e');
+    }
+
+    KeyMappings(
+            unsigned char mv_fwd,
+            unsigned char mv_bck,
+            unsigned char mv_left,
+            unsigned char mv_right,
+            unsigned char turn_left,
+            unsigned char turn_right)
+        :
+            mv_fwd(mv_fwd),
+            mv_bck(mv_bck),
+            mv_left(mv_left),
+            mv_right(mv_right),
+            turn_left(turn_left),
+            turn_right(turn_right)
+    {}
+
+
+    unsigned char
+        mv_fwd,
+        mv_bck,
+        mv_left,
+        mv_right,
+        turn_left,
+        turn_right;
+};
+
+void init_periodic_static(Ball* ball, GlutRender* render);
+void set_key_mapping(const KeyMappings& mapping);
 void periodic_update();
+
+void periodic_kb_handler(unsigned char key, int, int);
+void periodic_kb_up_handler(unsigned char key, int, int);

render/GlutRender.cpp

@@ -12,7 +12,18 @@ GlutRender& GlutRender::get_instance() {
     return instance;
 }
 
-GlutRender::GlutRender() {
+Mesh GlutRender::SurfaceDetails::render() const {
+    return MarchingCubes(*surface, box,
+            resolution)
+        .add_hint(surface->location_hint())
+        ();
+}
+
+GlutRender::GlutRender()
+    : followed_implicit(nullptr),
+    camera_position(0., 2.5, -10.),
+    camera_sight(0, 2.5, 0)
+{
     std::default_random_engine rand_engine(time(NULL));
     std::uniform_real_distribution<double> distribution;
     rand_color = std::bind(distribution, rand_engine);
@@ -29,12 +40,14 @@ void GlutRender::init(int* argc, char** argv,
     // ==== Callbacks ====
     glutDisplayFunc(display_handle);
     glutReshapeFunc(reshape_handle);
+    glutKeyboardFunc(kb_evt_handle);
+    glutKeyboardUpFunc(kb_evt_up_handle);
 
     // ==== Lighting ====
-    GLfloat light0_pos[] = {10., 15., 10.};
+    GLfloat light0_pos[] = {30., 35., 20., 0.};
     GLfloat light0_ambient[] = {0., 0., 0., 1.};
-    GLfloat light0_diffuse[] = {1., 1., .65, 1.};
-    GLfloat light0_specular[] = {5., 5., .33, 1.};
+    GLfloat light0_diffuse[] = {1., 1., .85, 1.};
+    GLfloat light0_specular[] = {5., 5., .43, 1.};
     glLightfv(GL_LIGHT0, GL_POSITION, light0_pos);
     glLightfv(GL_LIGHT0, GL_AMBIENT, light0_ambient);
     glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse);
@@ -47,15 +60,16 @@ void GlutRender::init(int* argc, char** argv,
     glEnable(GL_LIGHT0);
     glShadeModel(GL_SMOOTH); // Enable smooth shading
 
+    glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
     GLfloat material_specular[] = {1., 1., 1., 1.};
     GLfloat material_emission[] = {0., 0., 0., 1.};
-    glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
     glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, material_specular);
     glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, material_emission);
     glEnable(GL_COLOR_MATERIAL);
 
     // ==== Misc ====
-    glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // Background color
+    //glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // Background color
+    glClearColor(0.15f, 0.08f, 0.5f, 1.0f); // Background color
     glClearDepth(1.0f);                   // Set background depth to farthest
     glEnable(GL_DEPTH_TEST);
     glDepthFunc(GL_LEQUAL);
@@ -73,12 +87,19 @@ void GlutRender::remove_mesh(Mesh* mesh) {
     meshes.erase(mesh);
 }
 
-void GlutRender::add_surface(ImplicitSurface* surf, const Cuboid& box) {
-    surfaces.insert(SurfaceDetails(surf, box));
+void GlutRender::add_surface(ImplicitSurface* surf, const Cuboid& box,
+        bool always_render, double resolution)
+{
+    surfaces.push_back(SurfaceDetails(surf, box, resolution, always_render));
 }
 
 void GlutRender::remove_surface(ImplicitSurface* surf) {
-    surfaces.erase(SurfaceDetails(surf, Cuboid::empty()));
+    for(auto it=surfaces.begin(); it != surfaces.end(); ++it) {
+        if(*it == SurfaceDetails(surf, Cuboid::empty(), 0.1, true)) {
+            surfaces.erase(it);
+            break;
+        }
+    }
 }
 
 void GlutRender::run() {
@@ -89,6 +110,23 @@ void GlutRender::set_idle_func(void (*func)(void)) {
     glutIdleFunc(func);
 }
 
+void GlutRender::add_kb_handler(GlutRender::kb_handler_t handler) {
+    kb_handlers.push_back(handler);
+}
+
+void GlutRender::add_kb_up_handler(GlutRender::kb_handler_t handler) {
+    kb_up_handlers.push_back(handler);
+}
+
+void GlutRender::follow_implicit_position(const ImplicitSurface* surf) {
+    followed_implicit = surf;
+}
+
+void GlutRender::set_camera(const Point& location, const Point& sight) {
+    camera_position = location;
+    camera_sight = sight;
+}
+
 void GlutRender::reshape(int wid, int hei) {
     if (hei == 0)
         hei = 1;
@@ -149,8 +187,15 @@ void GlutRender::display_mesh(Mesh& mesh) const {
 #endif // DEBUG_DISPLAY_WIREFRAME
 
 
+    const Color& color = mesh.get_color();
+    GLfloat material_specular[] = {
+        (float)color[0],
+        (float)color[1],
+        (float)color[2],
+        1.};
     glBegin(GL_TRIANGLES);
-    glColor3f(1., 1., 1.);
+    glColor3f(color[0], color[1], color[2]);
+    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, material_specular);
     for(const Face& face: mesh.get_faces()) {
         Point p0 = face.pt(0, points) + mesh_center,
               p1 = face.pt(1, points) + mesh_center,
@@ -172,26 +217,54 @@ void GlutRender::display() {
 
     // Camera position and orientation
     glLoadIdentity();
-    //glTranslatef(1., 2., -10.);
-    gluLookAt(0, 2.5, -10,
-            0, 2.5, 0,
-            0, 1, 0);
+
+    if(followed_implicit != nullptr) {
+        const Point& look_at = followed_implicit->getCenter();
+        gluLookAt(camera_position.x, camera_position.y, camera_position.z,
+                look_at.x, camera_position.y, look_at.z,
+                0, 1, 0);
+    }
+    else {
+        gluLookAt(camera_position.x, camera_position.y, camera_position.z,
+                camera_sight.x, camera_position.y, camera_sight.z,
+                0, 1, 0);
+    }
 
     for(Mesh* mesh: meshes) {
         display_mesh(*mesh);
     }
 
-    for(const SurfaceDetails& surface: surfaces) {
-        Mesh mesh = MarchingCubes(*surface.surface, surface.box)();
-        display_mesh(mesh);
+    for(SurfaceDetails& surface: surfaces) {
+        if(!surface.always_render) {
+            if(surface.prerender == nullptr)
+                surface.self_render();
+            display_mesh(*surface.prerender);
+        }
+        else {
+            Mesh mesh = surface.render();
+            display_mesh(mesh);
+        }
     }
 
     glutSwapBuffers();
 }
 
+void GlutRender::on_kb_evt(bool up, unsigned char key, int x, int y) {
+    std::vector<kb_handler_t>& handlers = up? kb_up_handlers : kb_handlers;
+    for(auto& handler: handlers) {
+        handler(key, x, y);
+    }
+}
+
 void GlutRender::reshape_handle(int wid, int hei) {
     get_instance().reshape(wid, hei);
 }
 void GlutRender::display_handle() {
     get_instance().display();
 }
+void GlutRender::kb_evt_handle(unsigned char key, int x, int y) {
+    get_instance().on_kb_evt(false, key, x, y);
+}
+void GlutRender::kb_evt_up_handle(unsigned char key, int x, int y) {
+    get_instance().on_kb_evt(true, key, x, y);
+}

render/GlutRender.hpp

@@ -5,9 +5,13 @@
 #include "../Mesh.hpp"
 #include "../Implicit.hpp"
 #include <set>
+#include <vector>
+#include <functional>
 
 class GlutRender {
     public:
+        typedef std::function<void(unsigned char, int, int)> kb_handler_t;
+
         static GlutRender& get_instance();
 
         GlutRender(GlutRender const&) = delete;
@@ -18,22 +22,46 @@ class GlutRender {
         void cleanup();
         void add_mesh(Mesh* mesh);
         void remove_mesh(Mesh* mesh);
-        void add_surface(ImplicitSurface* surf, const Cuboid& box);
+        void add_surface(ImplicitSurface* surf, const Cuboid& box,
+                bool always_render=true, double resolution=0.1);
         void remove_surface(ImplicitSurface* surf);
         void run();
 
         void set_idle_func(void (*func)(void));
 
+        void add_kb_handler(kb_handler_t handler);
+        void add_kb_up_handler(kb_handler_t handler);
+
+        void follow_implicit_position(const ImplicitSurface* surf);
+
+        void set_camera(const Point& location, const Point& sight);
+
     private:
         struct SurfaceDetails {
-            SurfaceDetails(ImplicitSurface* surf, const Cuboid& box):
-                surface(surf), box(box) {}
+            SurfaceDetails(ImplicitSurface* surf, const Cuboid& box,
+                    double resolution,
+                    bool always_render):
+                surface(surf), box(box), resolution(resolution),
+                always_render(always_render), prerender(nullptr)
+            {}
+            ~SurfaceDetails() {
+                if(prerender != nullptr)
+                    delete prerender;
+            }
 
             ImplicitSurface* surface;
             Cuboid box;
+            double resolution;
 
-            bool operator<(const SurfaceDetails& oth) const {
-                return surface < oth.surface;
+            bool always_render;
+            Mesh* prerender;
+
+            bool operator==(const SurfaceDetails& oth) const {
+                return surface == oth.surface;
+            }
+            Mesh render() const;
+            void self_render() {
+                prerender = new Mesh(render());
             }
         };
 
@@ -46,13 +74,24 @@ class GlutRender {
     protected:
         void reshape(int wid, int hei);
         void display();
+        void on_kb_evt(bool up, unsigned char key, int x, int y);
 
         static void reshape_handle(int wid, int hei);
         static void display_handle();
 
+        static void kb_evt_handle(unsigned char key, int x, int y);
+        static void kb_evt_up_handle(unsigned char key, int x, int y);
+
     private: //attr
         std::function<double()> rand_color;
 
         std::set<Mesh*> meshes;
-        std::set<SurfaceDetails> surfaces;
+        std::vector<SurfaceDetails> surfaces;
+
+        std::vector<kb_handler_t> kb_handlers;
+        std::vector<kb_handler_t> kb_up_handlers;
+
+        const ImplicitSurface* followed_implicit;
+
+        Point camera_position, camera_sight;
 };

spheroid.cpp

@@ -6,7 +6,7 @@
 #include <iostream>
 
 Spheroid::Spheroid(const Point& _center, double _min_p, double _p, double _q) :
-    ImplicitSurface(_center), normal_vector(Point(0,0,1)), min_p(_min_p),
+    ImplicitSurface(_center), min_p(_min_p),
     init_p(_p), p(_p), q(_q), stiffness(0)
 {
     _compute_volume();
@@ -22,35 +22,19 @@ void Spheroid::update_radius() {
     q = sqrt((3./4.) * V / PI / p);
 }
 
-void Spheroid::update_height() {
-    q = (3./4.) * V / PI / (p*p);
-}
-
 
 void Spheroid::update_center_pos(Point& _center) {
     center = _center;
 }
 
-void Spheroid::check_horizontal_plan_collision(double height) {
-    if (((center.y - p) <= height) || (p < init_p)) {
-        p = fmin(init_p, center.y-height);
+void Spheroid::check_ground_collision(const Ground* ground) {
+    double height = (*ground)(center.x, center.z);
+    if (((center.y -p) <= height) || (p < init_p)) {
+        p = fmin(init_p, center.y - height);
         update_radius();
     }
 }
 
-void Spheroid::check_vertical_plan_collision(double& abscissa) {
-    if (center.x <= abscissa) {
-        if ((center.x + q) >= abscissa) {
-            q = abscissa - center.x;
-            update_height();
-        }
-    } else {
-        if ((center.x - q) <= abscissa) {
-            q = center.x - abscissa;
-            update_height();
-        }
-    }
-}
 
 Cuboid Spheroid::max_bounding_box() const {
     double max_radius = sqrt((3./4.) * V / PI / min_p);
@@ -60,18 +44,13 @@ Cuboid Spheroid::max_bounding_box() const {
     return Cuboid(_bd1, _bd2);
 }
 
-void Spheroid::check_perlin_collision(PerlinNoise perlin) {
-    double height = perlin.noise(center.x, center.z);
-    if ((center.y - p) <= height) {
-        p = fmin(init_p, center.z- height );
-        normal_vector = perlin.normal_vector(center.x, center.z);
-        update_radius();
-    }
-}
-
 
 double Spheroid::operator() (double _x, double _y, double _z) const {
             return (pow(_x, 2) / pow(q, 2)
                   + pow(_y, 2) / pow(p, 2)
                   + pow(_z, 2) / pow(q, 2) -1);
 }
+
+Point Spheroid::location_hint() const {
+    return Point(0, p, 0);
+}

spheroid.hpp

@@ -8,6 +8,8 @@
 #include "Implicit.hpp"
 #include "common_structures.hpp"
 #include "PerlinNoise.hpp"
+#include "FlatGround.hpp"
+#include "PerlinGround.hpp"
 
 const double PI = 3.141592653589793;
 
@@ -16,13 +18,12 @@ class Spheroid : public ImplicitSurface {
         Spheroid(const Point& _center, double _min_p, double _p, double _q);
         void update_center_pos(Point& _center);
         void update_radius();
-        void update_height();
-        void check_horizontal_plan_collision(double height);
-        void check_vertical_plan_collision(double& abscissa);
+        void check_ground_collision(const Ground* ground);
         Cuboid max_bounding_box() const;
         void check_perlin_collision(PerlinNoise perlin);
-        Point getNormalVector() const { return normal_vector;}
         double operator() (double _x, double _y, double _z) const;
+
+        virtual Point location_hint() const;
     private:
         /**
          * p corresponds to the half-height of the ball,
@@ -30,7 +31,6 @@ class Spheroid : public ImplicitSurface {
          * V is the volume. Extremely useful to have a constant volume in the
          * ball
         **/
-        Point normal_vector;
         double min_p, init_p, p, q;
         size_t stiffness;
         double V;

tests/TestImplicitSphere.cpp

@@ -7,3 +7,7 @@ double TestImplicitSphere::operator()(double x, double y, double z) const {
             + sq(center.z - z))
         + sq(radius);
 }
+
+Point TestImplicitSphere::location_hint() const {
+    return Point(center.x + radius, center.y, center.z);
+}

tests/TestImplicitSphere.hpp

@@ -8,6 +8,8 @@ class TestImplicitSphere: public ImplicitSurface {
             ImplicitSurface(center), radius(r) {}
         virtual double operator()(double x, double y, double z) const;
 
+        virtual Point location_hint() const;
+
     private:
         double radius;
 };