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

@@ -0,0 +1,118 @@
+#include "Ball.hpp"
+#include <iostream>
+#include <cmath>
+
+bool bouncing = true;
+double stop_bouncing;
+
+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),
+    bounce_number(0.0),
+    crt_time(0),
+    A(_center.y),
+    B(0),
+    U(sqrt(2 * G_CTE * _center.y)),
+    T(sqrt(2.0 * _center.y / G_CTE)),
+    v_x(_v_x),
+    v_z(_v_z)
+{
+}
+
+void Ball::_compute_pos(double dt) {
+    double height = (*ground)(Center.x, Center.z);
+    if (bouncing) {
+        Center.y = fmax(
+            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;
+        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_ground_collision(ground);
+}
+
+void Ball::_compute_T_n() {
+    double update = (2. * U / G_CTE);
+    T += update;
+}
+void Ball::_compute_B_n() {
+    B = G_CTE * T + U;
+}
+
+void Ball::_compute_A_n() {
+    A = - G_CTE * T * T / 2 - (U * T);
+}
+
+void Ball::_compute_U_n() {
+    U *= 0.8;
+}
+
+void Ball::_compute_v_x(Point normal) {
+    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 = 0;
+    if (normal.x + normal.z != 0) {
+        factor = normal.z / (4*(normal.x + normal.z));
+    }
+    v_z *= (0.5 + factor);
+}
+
+void Ball::_compute_parameters() {
+        Point normal = (*ground).get_normal(Center.x, Center.z);
+        bounce_number += 1;
+        _compute_U_n();
+        _compute_A_n();
+        _compute_B_n();
+        _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)) {
+            stop_bouncing = crt_time;
+            bouncing = false;
+        }
+    }
+    _compute_pos(dt);
+}
+
+std::ostream& operator<< (std::ostream &out, Ball const& data) {
+    Point center = data.getCenter();
+    out << "t:" << data.access_crt_time() << ":";
+    out << "T:" << data.accessT() << ":";
+    out << center.x << ':';
+    out << center.z << ':';
+    out << center.y << ':';
+    return out;
+}

Ball.hpp

@@ -0,0 +1,50 @@
+/**
+ * Header file defining a ball.
+ * A ball is the kind of real object, with a physical model, and the implicit
+ * surface.
+ * The physical model is strongly inspired from
+ * http://www.physics-online.info/book1/chapt1/sect1a/problem1a-12/Problem1A-12.htm,
+ * apart for the x-axis drift.
+ * For this one, the assumption done is that the speed is divided by two at
+ * each bounce, which is a reasonable assumption since it is the same for the
+ * vertical bounces.
+ **/
+#pragma once
+#include <cstddef>
+#include "spheroid.hpp"
+#include "FlatGround.hpp"
+#include "PerlinGround.hpp"
+
+#define G_CTE 9.81
+
+class Ball {
+    private:
+        Point Center;
+        Spheroid surface;
+        const Ground* ground;
+        double radius;
+        double init_h;
+        double min_height;
+        size_t bounce_number;
+        double crt_time;
+        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();
+        void _compute_B_n();
+        void _compute_U_n();
+        void _compute_T_n();
+    public:
+        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;}
+        double accessT() const { return T;}
+        double access_crt_time() const { return crt_time;}
+        Spheroid* get_surface() { return &surface; }
+};
+
+std::ostream& operator << (std::ostream &out, Ball const& data);

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.cpp

@@ -3,3 +3,18 @@
 double ImplicitSurface::operator()(const Point& pt) const {
     return operator()(pt.x, pt.y, pt.z);
 }
+
+Point ImplicitSurface::normal_at(const Point& pt) const {
+    // this is simply the gradient.
+    static const double d_dist = 1e-8;
+
+    Point normal(
+            (operator()(pt.x + d_dist, pt.y, pt.z) -
+             operator()(pt.x - d_dist, pt.y, pt.z)) / (2. * d_dist),
+            (operator()(pt.x, pt.y + d_dist, pt.z) -
+             operator()(pt.x, pt.y - d_dist, pt.z)) / (2. * d_dist),
+            (operator()(pt.x, pt.y, pt.z + d_dist) -
+             operator()(pt.x, pt.y, pt.z - d_dist)) / (2. * d_dist));
+    normal.normalize();
+    return normal;
+}

Implicit.hpp

@@ -6,4 +6,18 @@ class ImplicitSurface {
     public:
         virtual double operator() (double x, double y, double z) const = 0;
         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

@@ -1,17 +1,25 @@
 CXX=g++
-CXXFLAGS=-Wall -Wextra -O2 -std=c++14
+CXXFLAGS=$(ADD_FLAGS) -Wall -Wextra -O2 -std=c++14
 CXXLIBS=-lGL -lGLU -lglut
 
 # In `TARGET`, list the names of the `main_[stuff].cpp` you'd like to compile
 # into a `[stuff].bin`.
-TARGETS=glut
+TARGETS=glut ball bounce
 
 OBJS=Implicit.o \
 	 common_structures.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
 
@@ -20,18 +28,21 @@ all: $(TARGETS:=.bin)
 %.bin: main_%.o $(OBJS)
 	$(CXX) $(CXXFLAGS) $(CXXLIBS) $(OBJS) $< -o $@
 
-_gen/marching_cubes_data.cpp: tools/gen_marching_cubes_conf.py _gen
+_gen/marching_cubes_data.cpp: tools/gen_marching_cubes_conf.py
+	mkdir -p _gen
 	python3 $< > $@
 
 %.o: %.cpp
 	$(CXX) $(CXXFLAGS) -c $< -o $@
 
-_gen:
-	mkdir -p _gen
-
 ############################################################
 
+.PRECIOUS: %.o
+
 .PHONY: clean
 
 clean:
 	rm -rf $(OBJS) $(TARGETS:=.bin) _gen
+
+clean_render:
+	rm -rf render/*.o $(TARGETS:=.bin)

MarchingCubes.cpp

@@ -1,6 +1,11 @@
 #include "MarchingCubes.hpp"
 
 #include <cassert>
+#include <queue>
+
+#include "GL/gl.h"
+
+#include "GL/gl.h"
 
 MarchingCubes::MarchingCubes(
         const ImplicitSurface& surface,
@@ -11,14 +16,28 @@ 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;
 }
@@ -31,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
 {
@@ -85,8 +216,10 @@ bool MarchingCubes::march_at(double x, double y, double z, Mesh& output) {
         };
 
         size_t vert_ids[3];
-        for(int i=0; i < 3; ++i)
-            vert_ids[i] = output.add_vertice(verts[i]);
+        for(int i=0; i < 3; ++i) {
+            Point normal = surface.normal_at(verts[i]);
+            vert_ids[i] = output.add_vertice(verts[i], normal);
+        }
 
         for(int i=0; i < 3; ++i) {
             output.add_face(

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)
@@ -59,19 +66,16 @@ class MarchingCubes {
                 const Cuboid& box=Cuboid(
                     Point(-20, -20, -20),
                     Point(20, 20, 20)),
-                double resolution=.25);
+                double resolution=.1);
 
         /** 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.cpp

@@ -6,11 +6,22 @@ Mesh::Mesh()
 }
 
 size_t Mesh::add_vertice(const Point& pt) {
-    vertices.push_back(pt);
-    return vertices.size() - 1;
+    return inner_add_vertice(pt, NormalVect());
+}
+
+size_t Mesh::add_vertice(const Point& pt, const Point& normal) {
+    return inner_add_vertice(pt, NormalVect(normal, true));
 }
 
 void Mesh::add_face(const Face& face) {
+    for(int f_id=0; f_id < 3; ++f_id) {
+        if((unsigned)face[f_id] >= normals.size())
+            throw std::out_of_range("Vertice out of range for face");
+        if(!normals[face[f_id]].manual)
+            normals[face[f_id]].dirty = true;
+
+        faces_with_vert[face[f_id]].push_back(faces.size());
+    }
     faces.push_back(face);
 }
 void Mesh::add_face(size_t f1, size_t f2, size_t f3) {
@@ -50,3 +61,37 @@ const std::vector<Point>& Mesh::get_vertices() const {
 const std::vector<Face>& Mesh::get_faces() const {
     return faces;
 }
+
+const Point& Mesh::get_normal(size_t vert_id) {
+    if(vert_id >= normals.size())
+        throw std::out_of_range("Normal out of range");
+
+    NormalVect& norm = normals[vert_id];
+    if(norm.dirty) {
+        compute_normal(vert_id);
+        norm.dirty = false;
+    }
+    return norm.vect;
+}
+
+size_t Mesh::inner_add_vertice(const Point& pt, const Mesh::NormalVect& normal)
+{
+    vertices.push_back(pt);
+    normals.push_back(normal);
+    faces_with_vert.push_back(std::vector<size_t>());
+    return vertices.size() - 1;
+}
+
+void Mesh::compute_normal(size_t vert) {
+    Point normal(0., 0., 0.);
+    for(size_t f_id: faces_with_vert[vert]) {
+        const Face& face = faces[f_id];
+        Point e1 = vertices[face[1]] - vertices[face[0]],
+              e2 = vertices[face[2]] - vertices[face[0]];
+        normal += e1.cross_prod(e2);
+    }
+
+    normal.normalize();
+    normals[vert].dirty = false;
+    normals[vert].vect = normal;
+}

Mesh.hpp

@@ -15,6 +15,7 @@ class Mesh {
 
         /// Adds a fresh vertice at `pt`, and returns its ID for further use
         size_t add_vertice(const Point& pt);
+        size_t add_vertice(const Point& pt, const Point& normal);
 
         /// Creates a new face out of the three given point IDs
         void add_face(const Face& face);
@@ -31,15 +32,39 @@ class Mesh {
          */
         void normalize_center(bool keep_position=false);
 
+        /// Get the normal vector for vertice `vert_id`
+        const Point& get_normal(size_t vert_id);
+
         /// Gets the various vertices
         const std::vector<Point>& get_vertices() const;
 
         /// 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)
+                : vect(vect), dirty(false), manual(manual) {}
+            NormalVect()
+                : vect(Point(0., 0., 0.)), dirty(true), manual(false) {}
+            Point vect;
+            bool dirty;
+            bool manual;
+        };
+
+    private: //meth
+        size_t inner_add_vertice(const Point& pt, const NormalVect& normal);
+        void compute_normal(size_t vert);
+
     private:
         std::vector<Point> vertices;
+        std::vector<NormalVect> normals;
         std::vector<Face> faces;
+        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

@@ -0,0 +1,146 @@
+#include "PerlinNoise.hpp"
+
+Point vector_product(Point u, Point v) {
+    Point u_v(0,0,0);
+    u_v.x = u.y* v.z - u.z * v.y;
+    u_v.y = u.z* v.x - u.x * v.z;
+    u_v.z = u.x* v.y - u.y * v.x;
+    return u_v;
+}
+
+PerlinNoise::PerlinNoise() : ImplicitSurface(Point(0,0,0)) {
+    // permutation vector (source : http://mrl.nyu.edu/~perlin/noise/)
+    p = {
+		151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,
+		8,99,37,240,21,10,23,190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,
+		35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,
+		134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,
+		55,46,245,40,244,102,143,54, 65,25,63,161,1,216,80,73,209,76,132,187,208, 89,
+		18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,
+		250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,
+		189,28,42,223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167,
+		43,172,9,129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,
+		97,228,251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,
+		107,49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
+        138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180 };
+    // Duplication of the permutation vector.
+    p.insert(p.end(), p.begin(), p.end());
+}
+
+PerlinNoise::PerlinNoise(unsigned int seed) : ImplicitSurface(Point(0,0,0)) {
+	p.resize(256);
+
+	// Fill p with values from 0 to 255
+	std::iota(p.begin(), p.end(), 0);
+
+	// Initialize a random engine with seed
+	std::default_random_engine engine(seed);
+
+	// Shuffle  using the above random engine
+	std::shuffle(p.begin(), p.end(), engine);
+
+	// Duplicate the permutation vector
+	p.insert(p.end(), p.begin(), p.end());
+}
+
+double PerlinNoise::operator() (double x, double y, double z) const {
+    return (y - fBm(x, 0, z, 2, 0.3, 0.9));
+}
+
+double PerlinNoise::noise(double _x, double _y, double _z) const {
+    int X = (int) floor(_x) & 255;
+    int Y = (int) floor(_y) & 255;
+    int Z = (int) floor(_z) & 255;
+
+    _x -= floor(_x);
+    _y -= floor(_y);
+    _z -= floor(_z);
+
+    double u = fade(_x);
+    double v = fade(_y);
+    double w = fade(_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, _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, _z, _y-1),
+               grad(p[BA+1], _x-1, _z, _y-1)
+           ),
+           lerp(
+               u,
+               grad(p[AB+1], _x, _z-1, _y-1),
+               grad(p[BB+1], _x-1, _z-1, _y-1)
+           )
+       )
+    );
+	return (res + 1.0)/2.0;
+}
+
+double PerlinNoise::noise(double _x, double _y) const {
+    return noise(_x, 0, _y);
+}
+
+double PerlinNoise::fade(double t) const {
+    return t * t * t * (t * (t * 6 - 15) + 10);
+}
+
+double PerlinNoise::lerp(double t, double a, double b) const {
+    return a + t * (b-a);
+}
+
+double PerlinNoise::grad(int hash, double x, double y, double z)  const {
+	int h = hash & 15;
+	// Convert lower 4 bits of hash into 12 gradient directions
+	double u = h < 8 ? x : y,
+		   v = h < 4 ? y : h == 12 || h == 14 ? x : z;
+	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);
+    double n_x_dy = noise(x, y + 0.01, 0);
+    Point product = vector_product(
+       Point(0.01, 0, n_dx_y - n_x_y),
+       Point(0, 0.01, (n_x_dy - n_x_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

@@ -0,0 +1,29 @@
+/**
+ * Perlin Noise implementation for the ground (header file)
+ **/
+#pragma once
+#include "Implicit.hpp"
+#include <cmath>
+#include <random>
+#include <vector>
+#include <algorithm>
+
+class PerlinNoise : public ImplicitSurface {
+    std::vector<int> p;
+    public:
+        PerlinNoise();
+        PerlinNoise(unsigned int seed);
+        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;
+        double lerp(double t, double a, double b) const;
+        double grad(int hash, double x, double y, double z) 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

@@ -13,6 +13,17 @@ struct Point {
     Point(double x, double y, double z) : x(x), y(y), z(z) {}
     double x, y, z;
 
+    double norm() const {
+        auto sq = [](double x) { return x * x; };
+        return sqrt(sq(x) + sq(y) + sq(z));
+    }
+    void normalize() {
+        double _norm = norm();
+        x /= _norm;
+        y /= _norm;
+        z /= _norm;
+    }
+
     double operator[](unsigned i) const {
         assert(i < 3);
         switch(i) {
@@ -46,6 +57,10 @@ struct Point {
         return (*this += Point(-pt.x, -pt.y, -pt.z));
     }
 
+    Point operator-(const Point& pt) const {
+        return (*this) + (-1 * pt);
+    }
+
     friend Point operator*(double scalar, const Point& pt) {
         return Point(
                 scalar * pt.x,
@@ -53,6 +68,13 @@ struct Point {
                 scalar * pt.z);
     }
 
+    Point cross_prod(const Point& pt) const {
+        return Point(
+                (*this)[1] * pt[2] - (*this)[2] * pt[1],
+                (*this)[2] * pt[0] - (*this)[0] * pt[2],
+                (*this)[0] * pt[1] - (*this)[1] * pt[0]);
+    }
+
     bool operator<(const Point& pt) const {
         /// Lexicographic order on (x, y, z)
         if(x == pt.x) {
@@ -121,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

@@ -0,0 +1,22 @@
+/**
+ * Test file for the movement of a ball.
+ **/
+
+#include "Ball.hpp"
+#include <cstdio>
+#include <iostream>
+#include "FlatGround.hpp"
+
+using namespace std;
+
+int main() {
+    int i;
+    Point center(0,10,0);
+    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

@@ -0,0 +1,68 @@
+/** An entry-point file using render/GlutRender as a renderer, displaying the
+ * 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), 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);
+    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

@@ -0,0 +1,124 @@
+#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;
+
+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(
+        std::chrono::time_point<std::chrono::steady_clock> beg,
+        std::chrono::time_point<std::chrono::steady_clock> end)
+{
+    std::chrono::duration<double> ellapsed_db = end - beg;
+    return ellapsed_db.count();
+}
+
+void periodic_update() {
+    auto now = std::chrono::steady_clock::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

@@ -0,0 +1,45 @@
+#pragma once
+
+#include "Ball.hpp"
+#include "render/GlutRender.hpp"
+
+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

@@ -1,4 +1,5 @@
 #include "GlutRender.hpp"
+#include "../MarchingCubes.hpp"
 
 #include <GL/gl.h>
 #include <GL/glut.h>
@@ -11,26 +12,68 @@ 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);
+}
 
 void GlutRender::init(int* argc, char** argv,
         int wid, int hei, const char* win_name)
 {
-   glutInit(argc, argv);
-   glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH);
-   glutInitWindowSize(wid, hei);
-   glutCreateWindow(win_name);
+    glutInit(argc, argv);
+    glutInitDisplayMode(GLUT_DOUBLE | GLUT_DEPTH);
+    glutInitWindowSize(wid, hei);
+    glutCreateWindow(win_name);
 
-   glutDisplayFunc(display_handle);
-   glutReshapeFunc(reshape_handle);
+    // ==== Callbacks ====
+    glutDisplayFunc(display_handle);
+    glutReshapeFunc(reshape_handle);
+    glutKeyboardFunc(kb_evt_handle);
+    glutKeyboardUpFunc(kb_evt_up_handle);
 
-   glClearColor(0.0f, 0.0f, 0.0f, 1.0f); // Background color
-   glClearDepth(1.0f);                   // Set background depth to farthest
-   glEnable(GL_DEPTH_TEST);
-   glDepthFunc(GL_LEQUAL);
-   glShadeModel(GL_SMOOTH); // Enable smooth shading
-   glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
+    // ==== Lighting ====
+    GLfloat light0_pos[] = {30., 35., 20., 0.};
+    GLfloat light0_ambient[] = {0., 0., 0., 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);
+    glLightfv(GL_LIGHT0, GL_SPECULAR, light0_specular);
 
+    GLfloat light_ambient[] = {.2, .2, .2, 1.};
+    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, light_ambient);
+
+    glEnable(GL_LIGHTING);
+    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.};
+    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.15f, 0.08f, 0.5f, 1.0f); // Background color
+    glClearDepth(1.0f);                   // Set background depth to farthest
+    glEnable(GL_DEPTH_TEST);
+    glDepthFunc(GL_LEQUAL);
+    glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
 }
 
 void GlutRender::cleanup() {
@@ -44,10 +87,46 @@ void GlutRender::remove_mesh(Mesh* mesh) {
     meshes.erase(mesh);
 }
 
+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) {
+    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() {
     glutMainLoop();
 }
 
+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;
@@ -61,42 +140,131 @@ void GlutRender::reshape(int wid, int hei) {
     gluPerspective(45.0f, aspect, 0.1f, 100.0f);
 }
 
+void GlutRender::display_mesh(Mesh& mesh) const {
+    const Point& mesh_center = mesh.get_center();
+
+    const std::vector<Point>& points = mesh.get_vertices();
+
+#ifdef DEBUG_DISPLAY_NORMAL
+
+    glBegin(GL_LINES);
+    for(size_t vert_id=0; vert_id < points.size(); ++vert_id) {
+        const Point& pt = points[vert_id];
+        Point normal = mesh.get_normal(vert_id);
+        Point locVert = pt + mesh.get_center();
+        Point outwards = locVert + normal;
+        Point inwards = locVert + ((-0.2) * normal);
+        glColor3f(1., 0., 0.);
+        glVertex3f(locVert[0], locVert[1], locVert[2]);
+        glVertex3f(outwards[0], outwards[1], outwards[2]);
+        glColor3f(0., 0., 1.);
+        glVertex3f(locVert[0], locVert[1], locVert[2]);
+        glVertex3f(inwards[0], inwards[1], inwards[2]);
+    }
+    glEnd();
+#endif // DEBUG_DISPLAY_NORMAL
+
+#ifdef DEBUG_DISPLAY_WIREFRAME
+    glBegin(GL_LINES);
+    for(const Face& face: mesh.get_faces()) {
+        Point n0 = mesh.get_normal(face[0]),
+              n1 = mesh.get_normal(face[1]),
+              n2 = mesh.get_normal(face[2]);
+        Point p0 = face.pt(0, points) + mesh_center + 0.01 * n0,
+              p1 = face.pt(1, points) + mesh_center + 0.01 * n1,
+              p2 = face.pt(2, points) + mesh_center + 0.01 * n2;
+        glColor3f(0., 1., 0.);
+        glVertex3f(p0[0], p0[1], p0[2]);
+        glVertex3f(p1[0], p1[1], p1[2]);
+
+        glVertex3f(p1[0], p1[1], p1[2]);
+        glVertex3f(p2[0], p2[1], p2[2]);
+
+        glVertex3f(p2[0], p2[1], p2[2]);
+        glVertex3f(p0[0], p0[1], p0[2]);
+    }
+    glEnd();
+#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(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,
+              p2 = face.pt(2, points) + mesh_center;
+        Point n0 = mesh.get_normal(face[0]),
+              n1 = mesh.get_normal(face[1]),
+              n2 = mesh.get_normal(face[2]);
+
+        glNormal3f(n0[0], n0[1], n0[2]); glVertex3f(p0[0], p0[1], p0[2]);
+        glNormal3f(n1[0], n1[1], n1[2]); glVertex3f(p1[0], p1[1], p1[2]);
+        glNormal3f(n2[0], n2[1], n2[2]); glVertex3f(p2[0], p2[1], p2[2]);
+    }
+    glEnd();
+}
+
 void GlutRender::display() {
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
     glMatrixMode(GL_MODELVIEW);
 
     // Camera position and orientation
     glLoadIdentity();
-    glTranslatef(0., 0., -6.);
 
-    std::default_random_engine rand_engine(time(NULL));
-    std::uniform_real_distribution<double> distribution;
-    auto rand_color = std::bind(distribution, rand_engine);
+    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) {
-        const Point& mesh_center = mesh->get_center();
+        display_mesh(*mesh);
+    }
 
-        const std::vector<Point>& points = mesh->get_vertices();
-
-        glBegin(GL_TRIANGLES);
-        for(const Face& face: mesh->get_faces()) {
-            Point p0 = face.pt(0, points) + mesh_center,
-                  p1 = face.pt(1, points) + mesh_center,
-                  p2 = face.pt(2, points) + mesh_center;
-            glColor3f(rand_color(), rand_color(), rand_color());
-            glVertex3f(p0[0], p0[1], p0[2]);
-            glVertex3f(p1[0], p1[1], p1[2]);
-            glVertex3f(p2[0], p2[1], p2[2]);
+    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);
         }
-        glEnd();
     }
 
     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

@@ -3,10 +3,15 @@
 #pragma once
 
 #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;
@@ -17,17 +22,76 @@ class GlutRender {
         void cleanup();
         void add_mesh(Mesh* mesh);
         void remove_mesh(Mesh* mesh);
+        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,
+                    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 always_render;
+            Mesh* prerender;
+
+            bool operator==(const SurfaceDetails& oth) const {
+                return surface == oth.surface;
+            }
+            Mesh render() const;
+            void self_render() {
+                prerender = new Mesh(render());
+            }
+        };
+
+
     private: //meth
         GlutRender();
 
+        void display_mesh(Mesh& mesh) const;
+
     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();
-    private:
+
+        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::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

@@ -0,0 +1,56 @@
+/**
+ * Implementation of a spheroid
+ **/
+
+#include "spheroid.hpp"
+#include <iostream>
+
+Spheroid::Spheroid(const Point& _center, double _min_p, double _p, double _q) :
+    ImplicitSurface(_center), min_p(_min_p),
+    init_p(_p), p(_p), q(_q), stiffness(0)
+{
+    _compute_volume();
+}
+
+void Spheroid::_compute_volume() {
+    V = (4./3.) * PI * p * q * q;
+}
+
+
+
+void Spheroid::update_radius() {
+    q = sqrt((3./4.) * V / PI / p);
+}
+
+
+void Spheroid::update_center_pos(Point& _center) {
+    center = _center;
+}
+
+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();
+    }
+}
+
+
+Cuboid Spheroid::max_bounding_box() const {
+    double max_radius = sqrt((3./4.) * V / PI / min_p);
+    double max_height = init_p;
+    Point _bd1(max_radius, max_height, max_radius);
+    Point _bd2(-max_radius, -max_height, -max_radius);
+    return Cuboid(_bd1, _bd2);
+}
+
+
+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

@@ -0,0 +1,38 @@
+/**
+ * Defines a spheroid, which is a basic interpretaion of the ball when it is
+ * bouncing.
+ **/
+
+#include <cstddef>
+#include <cmath>
+#include "Implicit.hpp"
+#include "common_structures.hpp"
+#include "PerlinNoise.hpp"
+#include "FlatGround.hpp"
+#include "PerlinGround.hpp"
+
+const double PI = 3.141592653589793;
+
+class Spheroid : public ImplicitSurface {
+    public:
+        Spheroid(const Point& _center, double _min_p, double _p, double _q);
+        void update_center_pos(Point& _center);
+        void update_radius();
+        void check_ground_collision(const Ground* ground);
+        Cuboid max_bounding_box() const;
+        void check_perlin_collision(PerlinNoise perlin);
+        double operator() (double _x, double _y, double _z) const;
+
+        virtual Point location_hint() const;
+    private:
+        /**
+         * p corresponds to the half-height of the ball,
+         * q to the radius of the ball,
+         * V is the volume. Extremely useful to have a constant volume in the
+         * ball
+        **/
+        double min_p, init_p, p, q;
+        size_t stiffness;
+        double V;
+        void _compute_volume();
+};

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

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