Initial commit
authorFrancois Fleuret <francois@fleuret.org>
Fri, 7 Oct 2016 10:03:24 +0000 (12:03 +0200)
committerFrancois Fleuret <francois@fleuret.org>
Fri, 7 Oct 2016 10:03:24 +0000 (12:03 +0200)
13 files changed:
Makefile [new file with mode: 0644]
canvas.cc [new file with mode: 0644]
canvas.h [new file with mode: 0644]
canvas_cairo.cc [new file with mode: 0644]
canvas_cairo.h [new file with mode: 0644]
flatland.cc [new file with mode: 0644]
misc.cc [new file with mode: 0644]
misc.h [new file with mode: 0644]
polygon.cc [new file with mode: 0644]
polygon.h [new file with mode: 0644]
run.sh [new file with mode: 0755]
universe.cc [new file with mode: 0644]
universe.h [new file with mode: 0644]

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..304f403
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,51 @@
+
+# dyncnn is a deep-learning algorithm for the prediction of
+# interacting object dynamics
+#
+# Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+# Written by Francois Fleuret <francois.fleuret@idiap.ch>
+#
+# This file is part of dyncnn.
+#
+# dyncnn is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License version 3 as
+# published by the Free Software Foundation.
+#
+# dyncnn is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+
+ifeq ($(DEBUG),yes)
+ CXXFLAGS = -fPIC -Wall -g -DDEBUG
+else
+ CXXFLAGS = -fPIC -Wall -g -O3
+endif
+
+LDFLAGS = -lz -ldl
+
+CXXFLAGS += -I/usr/include/cairo -DCAIRO_SUPPORT
+
+LDFLAGS += -lcairo
+
+all: flatland TAGS
+
+TAGS: *.cc *.h
+       etags *.cc *.h
+
+flatland:      flatland.o misc.o \
+       universe.o \
+       polygon.o \
+       canvas.o canvas_cairo.o
+       $(CXX) $(CXXFLAGS) $(LDFLAGS) -lX11 -o $@ $^
+
+Makefile.depend: *.h *.cc Makefile
+       $(CC) $(CXXFLAGS) -M *.cc > Makefile.depend
+
+clean:
+       \rm -f flatland *.o *.so Makefile.depend
+
+-include Makefile.depend
diff --git a/canvas.cc b/canvas.cc
new file mode 100644 (file)
index 0000000..358bd7e
--- /dev/null
+++ b/canvas.cc
@@ -0,0 +1,25 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "canvas.h"
diff --git a/canvas.h b/canvas.h
new file mode 100644 (file)
index 0000000..f54f6d0
--- /dev/null
+++ b/canvas.h
@@ -0,0 +1,37 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef CANVAS_H
+#define CANVAS_H
+
+#include "misc.h"
+
+class Canvas {
+public:
+  virtual void set_line_width(scalar_t w) = 0;
+  virtual void set_drawing_color(scalar_t r, scalar_t g, scalar_t b) = 0;
+  virtual void draw_polygon(int filled, int nb, scalar_t *x, scalar_t *y) = 0;
+};
+
+#endif
diff --git a/canvas_cairo.cc b/canvas_cairo.cc
new file mode 100644 (file)
index 0000000..ec26f60
--- /dev/null
@@ -0,0 +1,91 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "canvas_cairo.h"
+
+CanvasCairo::CanvasCairo(scalar_t scale, int width, int height) {
+  const int actual_width = int(width * scale);
+  const int actual_height = int(height * scale);
+  const int depth = 4;
+
+  _data = new unsigned char [actual_width * actual_height * depth];
+
+  _image = cairo_image_surface_create_for_data(_data,
+                                               CAIRO_FORMAT_RGB24,
+                                               actual_width,
+                                               actual_height,
+                                               actual_width * depth);
+
+  _context_resource = cairo_create(_image);
+
+  cairo_scale(_context_resource, scale, scale);
+
+  cairo_set_source_rgb(_context_resource, 1.0, 1.0, 1.0);
+  cairo_set_line_width (_context_resource, 1.0);
+
+  cairo_rectangle(_context_resource, 0, 0, width, height);
+
+  cairo_fill(_context_resource);
+}
+
+CanvasCairo::~CanvasCairo() {
+  cairo_destroy(_context_resource);
+  cairo_surface_destroy(_image);
+  delete[] _data;
+}
+
+
+void CanvasCairo::set_line_width(scalar_t w) {
+  cairo_set_line_width (_context_resource, w);
+}
+
+void CanvasCairo::set_drawing_color(scalar_t r, scalar_t g, scalar_t b) {
+  cairo_set_source_rgb(_context_resource, r, g, b);
+}
+
+void CanvasCairo::draw_polygon(int filled, int nb, scalar_t *x, scalar_t *y) {
+  cairo_move_to(_context_resource, x[0], y[0]);
+  for(int n = 0; n < nb; n++) {
+    cairo_line_to(_context_resource, x[n], y[n]);
+  }
+  cairo_close_path(_context_resource);
+
+  if(filled) {
+    cairo_stroke_preserve(_context_resource);
+    cairo_fill(_context_resource);
+  } else {
+    cairo_stroke(_context_resource);
+  }
+}
+
+static cairo_status_t write_cairo_to_file(void *closure,
+                                          const unsigned char *data,
+                                          unsigned int length) {
+  fwrite(data, 1, length, (FILE *) closure);
+  return CAIRO_STATUS_SUCCESS;
+}
+
+void CanvasCairo::write_png(FILE *file) {
+  cairo_surface_write_to_png_stream(_image, write_cairo_to_file, file);
+}
diff --git a/canvas_cairo.h b/canvas_cairo.h
new file mode 100644 (file)
index 0000000..3b8f06b
--- /dev/null
@@ -0,0 +1,48 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef CANVAS_CAIRO_H
+#define CANVAS_CAIRO_H
+
+#include "canvas.h"
+
+#include <stdio.h>
+#include <cairo.h>
+
+class CanvasCairo : public Canvas {
+  unsigned char *_data;
+  cairo_surface_t *_image;
+  cairo_t *_context_resource;
+
+public:
+  CanvasCairo(scalar_t scale, int width, int height);
+  ~CanvasCairo();
+
+  virtual void set_line_width(scalar_t w);
+  virtual void set_drawing_color(scalar_t r, scalar_t g, scalar_t b);
+  virtual void draw_polygon(int filled, int nb, scalar_t *x, scalar_t *y);
+  virtual void write_png(FILE *file);
+};
+
+#endif
diff --git a/flatland.cc b/flatland.cc
new file mode 100644 (file)
index 0000000..c27bd97
--- /dev/null
@@ -0,0 +1,283 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+#include <errno.h>
+#include <string.h>
+#include <sys/stat.h>
+
+using namespace std;
+
+#include "misc.h"
+#include "universe.h"
+#include "canvas_cairo.h"
+
+void generate_png(Universe *universe, scalar_t scale, FILE *file) {
+  CanvasCairo canvas(scale, universe->width(), universe->height());
+  canvas.set_line_width(1.0 / scale);
+  universe->draw(&canvas);
+  canvas.write_png(file);
+}
+
+FILE *safe_fopen(const char *name, const char *mode) {
+  FILE *file = fopen(name, mode);
+  if(!file) {
+    cerr << "Cannot open " << name << endl;
+    exit(1);
+  }
+  return file;
+}
+
+void print_help(const char *command) {
+  cerr << command << " <nb sequences to generate> [--dir <dir>] [--seed <seed>]]" << endl;
+  exit(1);
+}
+
+int main(int argc, char **argv) {
+  const scalar_t world_width = 400;
+  const scalar_t world_height = 400;
+  const scalar_t block_size = 80;
+
+  const scalar_t dt = 0.1;
+  const int nb_iterations_per_steps = 20;
+
+  //////////////////////////////////////////////////////////////////////
+
+  // We will generate images { 0, every_nth, 2 * every_nth, ..., nb_frames - 1 }
+
+  // The framerate every_nth may be set to smaller value to generate
+  // nice materials for presentations or papers.
+
+  int every_nth = 4;
+
+  int nb_frames = 5;
+
+  int multi_grasp = 0;
+  int nb_shapes = 1;
+  char data_dir[1024] = "/tmp/";
+
+  //////////////////////////////////////////////////////////////////////
+
+  Universe *universe;
+  Polygon *grabbed_polygon;
+
+  if(argc < 2) {
+    print_help(argv[0]);
+  }
+
+  int nb_sequences = atoi(argv[1]);
+
+  int i = 2;
+  while(i < argc) {
+    if(strcmp(argv[i], "--dir") == 0) {
+      i++;
+      if(i == argc) { print_help(argv[0]);}
+      strncpy(data_dir, argv[i], sizeof(data_dir) / sizeof(char) - 1);
+      i++;
+    }
+
+    else if(strcmp(argv[i], "--seed") == 0) {
+      i++;
+      if(i == argc) { print_help(argv[0]); }
+      srand48(atoi(argv[i]));
+      i++;
+    }
+
+    else if(strcmp(argv[i], "--nb_shapes") == 0) {
+      i++;
+      if(i == argc) { print_help(argv[0]);}
+      nb_shapes = atoi(argv[i]);
+      i++;
+    }
+
+    else if(strcmp(argv[i], "--multi_grasp") == 0) {
+      multi_grasp = 1;
+      i++;
+    }
+
+    else if(strcmp(argv[i], "--every_nth") == 0) {
+      i++;
+      if(i == argc) { print_help(argv[0]);}
+      every_nth = atoi(argv[i]);
+      i++;
+    }
+
+    else if(strcmp(argv[i], "--nb_frames") == 0) {
+      i++;
+      if(i == argc) { print_help(argv[0]);}
+      nb_frames = atoi(argv[i]);
+      i++;
+    }
+
+    else {
+      cerr << "Unknown option " << argv[i] << "." << endl;
+      abort();
+    }
+  }
+
+  if(nb_shapes < 1 || nb_shapes > 10) {
+    cerr << "nb_shapes has to be in {1, ..., 10}" << endl;
+    abort();
+  }
+
+  if(nb_frames < 0 || nb_frames > 50) {
+    cerr << "nb_frames has to be in {0, ..., 50}" << endl;
+    abort();
+  }
+
+  universe = new Universe(nb_shapes, world_width, world_height);
+
+  for(int n = 0; n < nb_sequences; n++) {
+
+    scalar_t grab_start_x = world_width * 0.5;
+    scalar_t grab_start_y = world_height * 0.75;
+
+    if(multi_grasp) {
+      grab_start_x = world_width * (0.1 + 0.8 * drand48());
+      grab_start_y = world_height * (0.1 + 0.8 * drand48());
+    }
+
+    if((n+1)%100 == 0) {
+      cout << "Created "
+           << n+1 << "/" << nb_sequences << " sequences in "
+           << data_dir
+           << "." << endl;
+    }
+
+    do {
+      universe->clear();
+
+      const int nb_attempts_max = 100;
+      int nb_attempts = 0;
+
+      for(int u = 0; u < nb_shapes; u++) {
+        Polygon *pol = 0;
+
+        nb_attempts = 0;
+
+        do {
+          scalar_t x[] = {
+            - block_size * 0.4,
+            + block_size * 0.4,
+            + block_size * 0.4,
+            - block_size * 0.4,
+          };
+
+          scalar_t y[] = {
+            - block_size * 0.6,
+            - block_size * 0.6,
+            + block_size * 0.6,
+            + block_size * 0.6,
+          };
+
+          scalar_t delta = block_size / sqrt(2.0);
+          scalar_t object_center_x = delta + (world_width - 2 * delta) * drand48();
+          scalar_t object_center_y = delta + (world_height - 2 * delta) * drand48();
+          scalar_t red, green, blue;
+          red = 1.00;
+          green = red;
+          blue = red;
+          delete pol;
+          pol = new Polygon(0.5,
+                            red, green, blue,
+                            x, y, sizeof(x)/sizeof(scalar_t));
+          pol->set_position(object_center_x, object_center_y, M_PI * 2 * drand48());
+          pol->set_speed(0, 0, 0);
+          universe->initialize_polygon(pol);
+          nb_attempts++;
+        } while(nb_attempts < nb_attempts_max && universe->collide(pol));
+
+        if(nb_attempts == nb_attempts_max) {
+          delete pol;
+          u = 0;
+          universe->clear();
+          nb_attempts = 0;
+        } else {
+          universe->add_polygon(pol);
+        }
+      }
+
+      grabbed_polygon = universe->pick_polygon(grab_start_x, grab_start_y);
+    } while(!grabbed_polygon);
+
+    const scalar_t scaling = 0.16;
+
+    CanvasCairo grab_trace(scaling, world_width, world_height);
+
+    {
+      char buffer[1024];
+      sprintf(buffer, "%s/%03d/", data_dir, n/1000);
+      mkdir(buffer, 0777);
+    }
+
+    scalar_t grab_relative_x = grabbed_polygon->relative_x(grab_start_x, grab_start_y);
+    scalar_t grab_relative_y = grabbed_polygon->relative_y(grab_start_x, grab_start_y);
+
+    {
+      int n = 36;
+      scalar_t xp[n], yp[n];
+      for(int k = 0; k < n; k++) {
+        scalar_t radius = 1/scaling;
+        scalar_t alpha = 2 * M_PI * scalar_t(k) / scalar_t(n);
+        xp[k] = grab_start_x + radius * cos(alpha);
+        yp[k] = grab_start_y + radius * sin(alpha);
+      }
+      grab_trace.set_drawing_color(0.0, 0.0, 0.0);
+      grab_trace.set_line_width(2.0);
+      grab_trace.draw_polygon(1, n, xp, yp);
+    }
+
+    for(int s = 0; s < nb_frames; s++) {
+      if(s % every_nth == 0) {
+        char buffer[1024];
+        sprintf(buffer, "%s/%03d/dyn_%06d_world_%03d.png", data_dir, n/1000, n, s);
+        FILE *file = safe_fopen(buffer, "w");
+        generate_png(universe, scaling, file);
+        fclose(file);
+      }
+
+      for(int i = 0; i < nb_iterations_per_steps; i++) {
+        scalar_t xf = grabbed_polygon->absolute_x(grab_relative_x, grab_relative_y);
+        scalar_t yf = grabbed_polygon->absolute_y(grab_relative_x, grab_relative_y);
+        grabbed_polygon->apply_force(dt, xf, yf, 0.0, -1.0);
+        universe->update(dt);
+      }
+    }
+
+    {
+      char buffer[1024];
+      sprintf(buffer, "%s/%03d/dyn_%06d_grab.png", data_dir, n/1000, n);
+      FILE *file = safe_fopen(buffer, "w");
+      grab_trace.write_png(file);
+      fclose(file);
+    }
+  }
+
+  delete universe;
+}
diff --git a/misc.cc b/misc.cc
new file mode 100644 (file)
index 0000000..84c010a
--- /dev/null
+++ b/misc.cc
@@ -0,0 +1,31 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "misc.h"
+
+int compare_couple(const void *a, const void *b) {
+  if(((Couple *) a)->value < ((Couple *) b)->value) return -1;
+  else if(((Couple *) a)->value > ((Couple *) b)->value) return 1;
+  else return 0;
+}
diff --git a/misc.h b/misc.h
new file mode 100644 (file)
index 0000000..0595f77
--- /dev/null
+++ b/misc.h
@@ -0,0 +1,53 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef MISC_H
+#define MISC_H
+
+#include <float.h>
+#include <stdlib.h>
+
+#ifdef DEBUG
+#define ASSERT(x, s) if(!(x)) { std::cerr << "ASSERT FAILED IN " << __FILE__ << ":" << __LINE__ << " [" << (s) << "]\n"; abort(); }
+#else
+#define ASSERT(x, s)
+#endif
+
+// typedef float scalar_t;
+typedef double scalar_t;
+
+inline scalar_t sq(scalar_t x) { return x*x; }
+
+inline scalar_t prod_vect(scalar_t x1, scalar_t y1, scalar_t x2, scalar_t y2) {
+  return x1 * y2 - x2 * y1;
+}
+
+struct Couple {
+  int index;
+  scalar_t value;
+};
+
+int compare_couple(const void *a, const void *b);
+
+#endif
diff --git a/polygon.cc b/polygon.cc
new file mode 100644 (file)
index 0000000..69b125e
--- /dev/null
@@ -0,0 +1,533 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <iostream>
+
+using namespace std;
+
+#include <cmath>
+#include "polygon.h"
+
+static const scalar_t dl = 20.0;
+static const scalar_t repulsion_constant = 0.2;
+static const scalar_t dissipation = 0.5;
+
+Polygon::Polygon(scalar_t mass,
+                 scalar_t red, scalar_t green, scalar_t blue,
+                 scalar_t *x, scalar_t *y,
+                 int nv) : _mass(mass),
+                           _moment_of_inertia(0), _radius(0),
+                           _relative_x(new scalar_t[nv]), _relative_y(new scalar_t[nv]),
+                           _total_nb_dots(0),
+                           _nb_dots(new int[nv]),
+                           _effecting_edge(0),
+                           _length(new scalar_t[nv]),
+                           _triangles(new Triangle[nv-2]),
+                           _initialized(false), _nailed(false),
+                           _nb_vertices(nv),
+                           _x(new scalar_t[nv]), _y(new scalar_t[nv]),
+                           _red(red), _green(green), _blue(blue) {
+  _last_center_x = 0;
+  _last_center_y = 0;
+  _last_theta = 0;
+
+  if(x) for(int i = 0; i < nv; i++) _relative_x[i] = x[i];
+  if(y) for(int i = 0; i < nv; i++) _relative_y[i] = y[i];
+}
+
+Polygon::~Polygon() {
+  delete[] _relative_x;
+  delete[] _relative_y;
+  delete[] _x;
+  delete[] _y;
+  delete[] _nb_dots;
+  delete[] _length;
+  delete[] _triangles;
+  delete[] _effecting_edge;
+}
+
+Polygon *Polygon::clone() {
+  return new Polygon(_mass, _red, _green, _blue, _relative_x, _relative_y, _nb_vertices);
+}
+
+#ifdef XFIG_SUPPORT
+void Polygon::color_xfig(XFigTracer *tracer) {
+  tracer->add_color(int(255 * _red), int(255 * _green), int(255 * _blue));
+}
+
+void Polygon::print_xfig(XFigTracer *tracer) {
+  tracer->draw_polygon(int(255 * _red), int(255 * _green), int(255 * _blue),
+                       _nb_vertices, _x, _y);
+}
+#endif
+
+#ifdef X11_SUPPORT
+void Polygon::draw(SimpleWindow *window) {
+  window->color(_red, _green, _blue);
+  int x[_nb_vertices], y[_nb_vertices];
+  for(int n = 0; n < _nb_vertices; n++) {
+    x[n] = int(_x[n]);
+    y[n] = int(_y[n]);
+  }
+  window->fill_polygon(_nb_vertices, x, y);
+}
+
+void Polygon::draw_contours(SimpleWindow *window) {
+  int x[_nb_vertices], y[_nb_vertices];
+  for(int n = 0; n < _nb_vertices; n++) {
+    x[n] = int(_x[n]);
+    y[n] = int(_y[n]);
+  }
+  window->color(0.0, 0.0, 0.0);
+  // window->color(1.0, 1.0, 1.0);
+  for(int n = 0; n < _nb_vertices; n++) {
+    window->draw_line(x[n], y[n], x[(n+1)%_nb_vertices], y[(n+1)%_nb_vertices]);
+  }
+}
+#endif
+
+void Polygon::draw(Canvas *canvas) {
+  canvas->set_drawing_color(_red, _green, _blue);
+  canvas->draw_polygon(1, _nb_vertices, _x, _y);
+}
+
+void Polygon::draw_contours(Canvas *canvas) {
+  canvas->set_drawing_color(0.0, 0.0, 0.0);
+  canvas->draw_polygon(0, _nb_vertices, _x, _y);
+}
+
+void Polygon::set_vertex(int k, scalar_t x, scalar_t y) {
+  _relative_x[k] = x;
+  _relative_y[k] = y;
+}
+
+void Polygon::set_position(scalar_t center_x, scalar_t center_y, scalar_t theta) {
+  _center_x = center_x;
+  _center_y = center_y;
+  _theta = theta;
+}
+
+void Polygon::set_speed(scalar_t dcenter_x, scalar_t dcenter_y, scalar_t dtheta) {
+  _dcenter_x = dcenter_x;
+  _dcenter_y = dcenter_y;
+  _dtheta = dtheta;
+}
+
+bool Polygon::contain(scalar_t x, scalar_t y) {
+  for(int t = 0; t < _nb_vertices-2; t++) {
+    scalar_t xa = _x[_triangles[t].a], ya = _y[_triangles[t].a];
+    scalar_t xb = _x[_triangles[t].b], yb = _y[_triangles[t].b];
+    scalar_t xc = _x[_triangles[t].c], yc = _y[_triangles[t].c];
+    if(prod_vect(x - xa, y - ya, xb - xa, yb - ya) <= 0 &&
+       prod_vect(x - xb, y - yb, xc - xb, yc - yb) <= 0 &&
+       prod_vect(x - xc, y - yc, xa - xc, ya - yc) <= 0) return true;
+  }
+  return false;
+}
+
+void Polygon::triangularize(int &nt, int nb, int *index) {
+  if(nb == 3) {
+
+    if(nt >= _nb_vertices-2) {
+      cerr << "Error type #1 in triangularization." << endl;
+      exit(1);
+    }
+
+    _triangles[nt].a = index[0];
+    _triangles[nt].b = index[1];
+    _triangles[nt].c = index[2];
+    nt++;
+
+  } else {
+    int best_m = -1, best_n = -1;
+    scalar_t best_split = -1, det, s = -1, t = -1;
+
+    for(int n = 0; n < nb; n++) for(int m = 0; m < n; m++) if(n > m+1 && m+nb > n+1) {
+      bool no_intersection = true;
+      for(int k = 0; no_intersection & (k < nb); k++)
+        if(k != n && k != m && (k+1)%nb != n && (k+1)%nb != m) {
+          intersection(_relative_x[index[n]], _relative_y[index[n]],
+                       _relative_x[index[m]], _relative_y[index[m]],
+                       _relative_x[index[k]], _relative_y[index[k]],
+                       _relative_x[index[(k+1)%nb]], _relative_y[index[(k+1)%nb]], det, s, t);
+          no_intersection = det == 0 || s < 0 || s > 1 || t < 0 || t > 1;
+        }
+
+      if(no_intersection) {
+        scalar_t a1 = 0, a2 = 0;
+        for(int k = 0; k < nb; k++) if(k >= m && k < n)
+          a1 += prod_vect(_relative_x[index[k]] - _relative_x[index[m]],
+                          _relative_y[index[k]] - _relative_y[index[m]],
+                          _relative_x[index[k+1]] - _relative_x[index[m]],
+                          _relative_y[index[k+1]] - _relative_y[index[m]]);
+        else
+          a2 += prod_vect(_relative_x[index[k]] - _relative_x[index[m]],
+                          _relative_y[index[k]] - _relative_y[index[m]],
+                          _relative_x[index[(k+1)%nb]] - _relative_x[index[m]],
+                          _relative_y[index[(k+1)%nb]] - _relative_y[index[m]]);
+
+        if((a1 * a2 > 0 && best_split < 0) || (abs(a1 - a2) < best_split)) {
+          best_n = n; best_m = m;
+          best_split = abs(a1 - a2);
+        }
+      }
+    }
+
+    if(best_n >= 0 && best_m >= 0) {
+      int index_neg[nb], index_pos[nb];
+      int neg = 0, pos = 0;
+      for(int k = 0; k < nb; k++) {
+        if(k >= best_m && k <= best_n) index_pos[pos++] = index[k];
+        if(k <= best_m || k >= best_n) index_neg[neg++] = index[k];
+      }
+      if(pos < 3 || neg < 3) {
+        cerr << "Error type #2 in triangularization." << endl;
+        exit(1);
+      }
+      triangularize(nt, pos, index_pos);
+      triangularize(nt, neg, index_neg);
+    } else {
+      cerr << "Error type #3 in triangularization." << endl;
+      exit(1);
+    }
+  }
+}
+
+void Polygon::initialize(int nb_polygons) {
+  scalar_t a;
+
+  _nb_polygons = nb_polygons;
+
+  a = _relative_x[_nb_vertices - 1] * _relative_y[0]
+    - _relative_x[0] * _relative_y[_nb_vertices - 1];
+
+  for(int n = 0; n < _nb_vertices - 1; n++)
+    a += _relative_x[n] * _relative_y[n+1] - _relative_x[n+1] * _relative_y[n];
+  a *= 0.5;
+
+  // Reorder the vertices
+
+  if(a < 0) {
+    a = -a;
+    scalar_t x, y;
+    for(int n = 0; n < _nb_vertices / 2; n++) {
+      x = _relative_x[n];
+      y = _relative_y[n];
+      _relative_x[n] = _relative_x[_nb_vertices - 1 - n];
+      _relative_y[n] = _relative_y[_nb_vertices - 1 - n];
+      _relative_x[_nb_vertices - 1 - n] = x;
+      _relative_y[_nb_vertices - 1 - n] = y;
+    }
+  }
+
+  // Compute the center of mass and moment of inertia
+
+  scalar_t sx, sy, w;
+  w = 0;
+  sx = 0;
+  sy = 0;
+  for(int n = 0; n < _nb_vertices; n++) {
+    int np = (n+1)%_nb_vertices;
+    w =_relative_x[n] * _relative_y[np] - _relative_x[np] * _relative_y[n];
+    sx += (_relative_x[n] + _relative_x[np]) * w;
+    sy += (_relative_y[n] + _relative_y[np]) * w;
+  }
+  sx /= 6 * a;
+  sy /= 6 * a;
+
+  _radius = 0;
+  for(int n = 0; n < _nb_vertices; n++) {
+    _relative_x[n] -= sx;
+    _relative_y[n] -= sy;
+    scalar_t r = sqrt(sq(_relative_x[n]) + sq(_relative_y[n]));
+    if(r > _radius) _radius = r;
+  }
+
+  scalar_t num = 0, den = 0;
+  for(int n = 0; n < _nb_vertices; n++) {
+    int np = (n+1)%_nb_vertices;
+    den += abs(prod_vect(_relative_x[np], _relative_y[np], _relative_x[n], _relative_y[n]));
+    num += abs(prod_vect(_relative_x[np], _relative_y[np], _relative_x[n], _relative_y[n])) *
+      (_relative_x[np] * _relative_x[np] + _relative_y[np] * _relative_y[np] +
+       _relative_x[np] * _relative_x[n] + _relative_y[np] * _relative_y[n] +
+       _relative_x[n] * _relative_x[n] + _relative_y[n] * _relative_y[n]);
+  }
+
+  _moment_of_inertia = num / (6 * den);
+
+  scalar_t vx = cos(_theta), vy = sin(_theta);
+
+  for(int n = 0; n < _nb_vertices; n++) {
+    _x[n] = _center_x + _relative_x[n] * vx + _relative_y[n] * vy;
+    _y[n] = _center_y - _relative_x[n] * vy + _relative_y[n] * vx;
+  }
+
+  scalar_t length;
+  _total_nb_dots = 0;
+  for(int n = 0; n < _nb_vertices; n++) {
+    length = sqrt(sq(_relative_x[n] - _relative_x[(n+1)%_nb_vertices]) +
+                  sq(_relative_y[n] - _relative_y[(n+1)%_nb_vertices]));
+    _length[n] = length;
+    _nb_dots[n] = int(length / dl + 1);
+    _total_nb_dots += _nb_dots[n];
+  }
+
+  delete[] _effecting_edge;
+  _effecting_edge = new int[_nb_polygons * _total_nb_dots];
+  for(int p = 0; p < _nb_polygons * _total_nb_dots; p++) _effecting_edge[p] = -1;
+
+  int nt = 0;
+  int index[_nb_vertices];
+  for(int n = 0; n < _nb_vertices; n++) index[n] = n;
+  triangularize(nt, _nb_vertices, index);
+
+  _initialized = true;
+}
+
+bool Polygon::update(scalar_t dt) {
+  if(!_nailed) {
+    _center_x += _dcenter_x * dt;
+    _center_y += _dcenter_y * dt;
+    _theta += _dtheta * dt;
+  }
+
+  scalar_t d = exp(log(dissipation) * dt);
+  _dcenter_x *= d;
+  _dcenter_y *= d;
+  _dtheta *= d;
+
+  scalar_t vx = cos(_theta), vy = sin(_theta);
+
+  for(int n = 0; n < _nb_vertices; n++) {
+    _x[n] = _center_x + _relative_x[n] * vx + _relative_y[n] * vy;
+    _y[n] = _center_y - _relative_x[n] * vy + _relative_y[n] * vx;
+  }
+
+  if(abs(_center_x - _last_center_x) +
+     abs(_center_y - _last_center_y) +
+     abs(_theta - _last_theta) * _radius > 0.1) {
+    _last_center_x = _center_x;
+    _last_center_y = _center_y;
+    _last_theta = _theta;
+    return true;
+  } else return false;
+}
+
+scalar_t Polygon::relative_x(scalar_t ax, scalar_t ay) {
+  return (ax - _center_x) * cos(_theta) - (ay - _center_y) * sin(_theta);
+}
+
+scalar_t Polygon::relative_y(scalar_t ax, scalar_t ay) {
+  return (ax - _center_x) * sin(_theta) + (ay - _center_y) * cos(_theta);
+}
+
+scalar_t Polygon::absolute_x(scalar_t rx, scalar_t ry) {
+  return _center_x + rx * cos(_theta) + ry * sin(_theta);
+}
+
+scalar_t Polygon::absolute_y(scalar_t rx, scalar_t ry) {
+  return _center_y - rx * sin(_theta) + ry * cos(_theta);
+}
+
+void Polygon::apply_force(scalar_t dt, scalar_t x, scalar_t y, scalar_t fx, scalar_t fy) {
+  _dcenter_x += fx / _mass * dt;
+  _dcenter_y += fy / _mass * dt;
+  _dtheta -= prod_vect(x - _center_x, y - _center_y, fx, fy) / (_mass * _moment_of_inertia) * dt;
+}
+
+void Polygon::apply_border_forces(scalar_t dt, scalar_t xmax, scalar_t ymax) {
+  for(int v = 0; v < _nb_vertices; v++) {
+    int vp = (v+1)%_nb_vertices;
+    for(int d = 0; d < _nb_dots[v]; d++) {
+      scalar_t s = scalar_t(d * dl)/_length[v];
+      scalar_t x = _x[v] * (1 - s) + _x[vp] * s;
+      scalar_t y = _y[v] * (1 - s) + _y[vp] * s;
+      scalar_t vx = 0, vy = 0;
+      if(x < 0) vx = x;
+      else if(x > xmax) vx = x - xmax;
+      if(y < 0) vy = y;
+      else if(y > ymax) vy = y - ymax;
+      apply_force(dt, x, y, - dl * vx * repulsion_constant, - dl * vy * repulsion_constant);
+    }
+  }
+}
+
+void Polygon::apply_collision_forces(scalar_t dt, int n_polygon, Polygon *p) {
+  scalar_t closest_x[_total_nb_dots], closest_y[_total_nb_dots];
+  bool inside[_total_nb_dots];
+  scalar_t distance[_total_nb_dots];
+  int _new_effecting_edge[_total_nb_dots];
+
+  int first_dot = 0;
+
+  for(int v = 0; v < _nb_vertices; v++) {
+    int vp = (v+1)%_nb_vertices;
+    scalar_t x = _x[v], y = _y[v], xp = _x[vp], yp = _y[vp];
+
+    for(int d = 0; d < _nb_dots[v]; d++) {
+      inside[d] = false;
+      distance[d] = FLT_MAX;
+    }
+
+    // First, we tag the dots located inside the polygon p
+
+    for(int t = 0; t < p->_nb_vertices-2; t++) {
+      scalar_t min = 0, max = 1;
+      scalar_t xa = p->_x[p->_triangles[t].a], ya = p->_y[p->_triangles[t].a];
+      scalar_t xb = p->_x[p->_triangles[t].b], yb = p->_y[p->_triangles[t].b];
+      scalar_t xc = p->_x[p->_triangles[t].c], yc = p->_y[p->_triangles[t].c];
+      scalar_t den, num;
+
+      const scalar_t eps = 1e-6;
+
+      den = prod_vect(xp - x, yp - y, xb - xa, yb - ya);
+      num = prod_vect(xa - x, ya - y, xb - xa, yb - ya);
+      if(den > eps) {
+        if(num / den < max) max = num / den;
+      } else if(den < -eps) {
+        if(num / den > min) min = num / den;
+      } else {
+        if(num < 0) { min = 1; max = 0; }
+      }
+
+      den = prod_vect(xp - x, yp - y, xc - xb, yc - yb);
+      num = prod_vect(xb - x, yb - y, xc - xb, yc - yb);
+      if(den > eps) {
+        if(num / den < max) max = num / den;
+      } else if(den < -eps) {
+        if(num / den > min) min = num / den;
+      } else {
+        if(num < 0) { min = 1; max = 0; }
+      }
+
+      den = prod_vect(xp - x, yp - y, xa - xc, ya - yc);
+      num = prod_vect(xc - x, yc - y, xa - xc, ya - yc);
+      if(den > eps) {
+        if(num / den < max) max = num / den;
+      } else if(den < -eps) {
+        if(num / den > min) min = num / den;
+      } else {
+        if(num < 0) { min = 1; max = 0; }
+      }
+
+      for(int d = 0; d < _nb_dots[v]; d++) {
+        scalar_t s = scalar_t(d * dl)/_length[v];
+        if(s >= min && s <= max) inside[d] = true;
+      }
+    }
+
+    // Then, we compute for each dot what is the closest point on
+    // the boundary of p
+
+    for(int m = 0; m < p->_nb_vertices; m++) {
+      int mp = (m+1)%p->_nb_vertices;
+      scalar_t xa = p->_x[m], ya = p->_y[m];
+      scalar_t xb = p->_x[mp], yb = p->_y[mp];
+      scalar_t gamma0 = ((x - xa) * (xb - xa) + (y - ya) * (yb - ya)) / sq(p->_length[m]);
+      scalar_t gamma1 = ((xp - x) * (xb - xa) + (yp - y) * (yb - ya)) / sq(p->_length[m]);
+      scalar_t delta0 = (prod_vect(xb - xa, yb - ya, x - xa, y - ya)) / p->_length[m];
+      scalar_t delta1 = (prod_vect(xb - xa, yb - ya, xp - x, yp - y)) / p->_length[m];
+
+      for(int d = 0; d < _nb_dots[v]; d++) if(inside[d]) {
+        int r = _effecting_edge[(first_dot + d) * _nb_polygons + n_polygon];
+
+        // If there is already a spring, we look only at the
+        // vertices next to the current one
+
+        if(r < 0 || m == r || m == (r+1)%p->_nb_vertices || (m+1)%p->_nb_vertices == r) {
+
+          scalar_t s = scalar_t(d * dl)/_length[v];
+          scalar_t delta = abs(s * delta1 + delta0);
+          if(delta < distance[d]) {
+            scalar_t gamma = s * gamma1 + gamma0;
+            if(gamma < 0) {
+              scalar_t l = sqrt(sq(x * (1 - s) + xp * s - xa) + sq(y * (1 - s) + yp * s - ya));
+              if(l < distance[d]) {
+                distance[d] = l;
+                closest_x[d] = xa;
+                closest_y[d] = ya;
+                _new_effecting_edge[first_dot + d] = m;
+              }
+            } else if(gamma > 1) {
+              scalar_t l = sqrt(sq(x * (1 - s) + xp * s - xb) + sq(y * (1 - s) + yp * s - yb));
+              if(l < distance[d]) {
+                distance[d] = l;
+                closest_x[d] = xb;
+                closest_y[d] = yb;
+                _new_effecting_edge[first_dot + d] = m;
+              }
+            } else {
+              distance[d] = delta;
+              closest_x[d] = xa * (1 - gamma) + xb * gamma;
+              closest_y[d] = ya * (1 - gamma) + yb * gamma;
+              _new_effecting_edge[first_dot + d] = m;
+            }
+          }
+        }
+      } else _new_effecting_edge[first_dot + d] = -1;
+    }
+
+    // Apply forces
+
+    for(int d = 0; d < _nb_dots[v]; d++) if(inside[d]) {
+      scalar_t s = scalar_t(d * dl)/_length[v];
+      scalar_t x = _x[v] * (1 - s) + _x[vp] * s;
+      scalar_t y = _y[v] * (1 - s) + _y[vp] * s;
+      scalar_t vx = x - closest_x[d];
+      scalar_t vy = y - closest_y[d];
+
+      p->apply_force(dt,
+                     closest_x[d], closest_y[d],
+                     dl * vx  * repulsion_constant, dl * vy * repulsion_constant);
+
+      apply_force(dt,
+                  x, y,
+                  - dl * vx * repulsion_constant, - dl * vy * repulsion_constant);
+    }
+
+    first_dot += _nb_dots[v];
+  }
+
+  for(int d = 0; d < _total_nb_dots; d++)
+    _effecting_edge[d * _nb_polygons + n_polygon] = _new_effecting_edge[d];
+
+}
+
+bool Polygon::collide(Polygon *p) {
+  for(int n = 0; n < _nb_vertices; n++) {
+    int np = (n+1)%_nb_vertices;
+    for(int m = 0; m < p->_nb_vertices; m++) {
+      int mp = (m+1)%p->_nb_vertices;
+      scalar_t det, s = -1, t = -1;
+      intersection(_x[n], _y[n], _x[np], _y[np],
+                   p->_x[m], p->_y[m], p->_x[mp], p->_y[mp], det, s, t);
+      if(det != 0 && s>= 0 && s <= 1&& t >= 0 && t <= 1) return true;
+    }
+  }
+
+  for(int n = 0; n < _nb_vertices; n++) if(p->contain(_x[n], _y[n])) return true;
+  for(int n = 0; n < p->_nb_vertices; n++) if(contain(p->_x[n], p->_y[n])) return true;
+
+  return false;
+}
diff --git a/polygon.h b/polygon.h
new file mode 100644 (file)
index 0000000..8e23a70
--- /dev/null
+++ b/polygon.h
@@ -0,0 +1,128 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef POLYGON_H
+#define POLYGON_H
+
+#include "misc.h"
+#include "canvas.h"
+
+#ifdef XFIG_SUPPORT
+#include "xfig_tracer.h"
+#endif
+
+#ifdef X11_SUPPORT
+#include "simple_window.h"
+#endif
+
+class Polygon {
+  struct Triangle {
+    int a, b, c;
+  };
+
+  int _nb_max_polygons;
+  int _nb_polygons; // We need to know the total number of polygons in
+                    // the universe
+
+  scalar_t _mass, _moment_of_inertia, _radius;
+
+  scalar_t *_relative_x, *_relative_y;
+  scalar_t _last_center_x, _last_center_y, _last_theta;
+
+  int _total_nb_dots;
+  int *_nb_dots;
+  int *_effecting_edge;
+
+  scalar_t *_length;
+  Triangle *_triangles;
+
+  inline void intersection(scalar_t xa2, scalar_t ya2,
+                           scalar_t xa1, scalar_t ya1,
+                           scalar_t xb2, scalar_t yb2,
+                           scalar_t xb1, scalar_t yb1,
+                           scalar_t &det, scalar_t &s1, scalar_t &s2) {
+    scalar_t m11, m12, m21, m22, n11, n12, n21, n22, v1, v2;
+    m11 = xa1 - xa2; m12 = xb2 - xb1; m21 = ya1 - ya2; m22 = yb2 - yb1;
+    det = m11 * m22 - m12 * m21;
+    if(det != 0) {
+      n11 = + m22 / det; n12 = - m12 / det; n21 = - m21 / det; n22 = + m11 / det;
+      v1 = xb2 - xa2; v2 = yb2 - ya2;
+      s1 = n11 * v1 + n12 * v2; s2 = n21 * v1 + n22 * v2;
+    }
+  }
+
+public:
+  bool _initialized;
+  bool _nailed;
+  int _nb_vertices;
+  scalar_t *_x, *_y;
+  scalar_t _red, _green, _blue;
+  scalar_t _center_x, _center_y, _theta;
+  scalar_t _dcenter_x, _dcenter_y, _dtheta;
+
+  // The x and y parameters are ignored if null. Useful to initialize
+  // directly.
+
+  Polygon(scalar_t mass,
+          scalar_t red, scalar_t green, scalar_t blue,
+          scalar_t *x, scalar_t *y,
+          int nv);
+
+  ~Polygon();
+
+  Polygon *clone();
+
+#ifdef XFIG_SUPPORT
+  void color_xfig(XFigTracer *tracer);
+  void print_xfig(XFigTracer *tracer);
+#endif
+
+#ifdef X11_SUPPORT
+  void draw(SimpleWindow *window);
+  void draw_contours(SimpleWindow *window);
+#endif
+
+  void draw(Canvas *canvas);
+  void draw_contours(Canvas *canvas);
+
+  void set_vertex(int k, scalar_t x, scalar_t y);
+  void set_position(scalar_t center_x, scalar_t center_y, scalar_t theta);
+  void set_speed(scalar_t dcenter_x, scalar_t dcenter_y, scalar_t dtheta);
+  bool contain(scalar_t x, scalar_t y);
+  void triangularize(int &nt, int nb, int *index);
+  void initialize(int nb_polygons);
+  bool update(scalar_t dt);
+  scalar_t relative_x(scalar_t ax, scalar_t ay);
+  scalar_t relative_y(scalar_t ax, scalar_t ay);
+  scalar_t absolute_x(scalar_t rx, scalar_t ry);
+  scalar_t absolute_y(scalar_t rx, scalar_t ry);
+
+  void apply_force(scalar_t dt, scalar_t x, scalar_t y, scalar_t fx, scalar_t fy);
+  void apply_border_forces(scalar_t dt, scalar_t xmax, scalar_t ymax);
+  void apply_collision_forces(scalar_t dt, int n_polygon, Polygon *p);
+
+  bool collide(Polygon *p);
+};
+
+#endif
diff --git a/run.sh b/run.sh
new file mode 100755 (executable)
index 0000000..d905ef9
--- /dev/null
+++ b/run.sh
@@ -0,0 +1,68 @@
+#!/bin/bash
+
+# dyncnn is a deep-learning algorithm for the prediction of
+# interacting object dynamics
+#
+# Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+# Written by Francois Fleuret <francois.fleuret@idiap.ch>
+#
+# This file is part of dyncnn.
+#
+# dyncnn is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License version 3 as
+# published by the Free Software Foundation.
+#
+# dyncnn is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+
+# This script creates the synthetic data-set for shape collision
+
+set -e
+set -o pipefail
+
+[[ "${TORCH_NB_THREADS}" ]] || echo "You can set \$TORCH_NB_THREADS to the proper value (default 1)."
+[[ "${TORCH_USE_GPU}" ]] || echo "You can set \$TORCH_USE_GPU to 'yes' or 'no' (default 'no')."
+[[ "${DYNCNN_DATA_DIR}" ]] || DYNCNN_DATA_DIR="./data/10p-mg"
+[[ "${DYNCNN_RESULT_DIR}" ]] || DYNCNN_RESULT_DIR="./results"
+
+if [[ -d "${DYNCNN_DATA_DIR}" ]]; then
+    echo "Found ${DYNCNN_DATA_DIR}, checking the number of images in there."
+    if [[ $(find "${DYNCNN_DATA_DIR}" -name "dyn_*.png" | wc -l) == 150000 ]]; then
+        echo "Looks good !"
+    else
+        echo "I do not find the proper number of images. Please remove the dir and re-run this scripts, or fix manually."
+        exit 1
+    fi
+else
+    # Creating the data-base
+    make -j -k
+    mkdir -p "${DYNCNN_DATA_DIR}"
+    ./flatland 50000 \
+               --every_nth 4 --nb_frames 5 \
+               --multi_grasp --nb_shapes 10 \
+               --dir "${DYNCNN_DATA_DIR}"
+fi
+
+# Train the model (takes days)
+
+if [[ ! -f "${DYNCNN_RESULT_DIR}"/epoch_01000_model ]]; then
+    ./dyncnn.lua --heavy --dataDir="${DYNCNN_DATA_DIR}" \
+                 --resultFreq=100 \
+                 --resultDir "${DYNCNN_RESULT_DIR}" \
+                 --nbEpochs 1000
+fi
+
+# Create the images of internal activations
+
+for n in 2 12; do
+    ./dyncnn.lua --heavy --dataDir=./data/10p-mg/ \
+                 --learningStateFile="${DYNCNN_RESULT_DIR}"/epoch_01000_model \
+                 --resultDir="${DYNCNN_RESULT_DIR}" \
+                 --noLog \
+                 --exampleInternals=${n}
+done
diff --git a/universe.cc b/universe.cc
new file mode 100644 (file)
index 0000000..658c32c
--- /dev/null
@@ -0,0 +1,196 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <string.h>
+
+#include "universe.h"
+
+Universe::Universe(int nb_max_polygons,
+                   scalar_t width, scalar_t height) : _width(width),
+                                                      _height(height),
+                                                      _nb_max_polygons(nb_max_polygons),
+                                                      _nb_polygons(0) {
+  _polygons = new Polygon *[_nb_max_polygons];
+  for(int n = 0; n < _nb_max_polygons; n++) _polygons[n] = 0;
+}
+
+Universe::~Universe() {
+  clear();
+  delete[] _polygons;
+}
+
+void Universe::initialize_polygon(Polygon *p) {
+  p->initialize(_nb_max_polygons);
+}
+
+void Universe::clear() {
+  for(int n = 0; n < _nb_polygons; n++) if(_polygons[n]) delete _polygons[n];
+  _nb_polygons = 0;
+}
+
+void Universe::add_polygon(Polygon *p) {
+  if(_nb_polygons < _nb_max_polygons) {
+    if(!p->_initialized) {
+      cerr << "You can not add a non-initialized polygon." << endl;
+      abort();
+    }
+    _polygons[_nb_polygons++] = p;
+  } else {
+    cerr << "Too many polygons!" << endl;
+    exit(1);
+  }
+}
+
+bool Universe::collide(Polygon *p) {
+  for(int n = 0; n < _nb_polygons; n++)
+    if(_polygons[n] && _polygons[n]->collide(p)) return true;
+
+  return false;
+}
+
+void Universe::compute_pseudo_collisions(int nb_axis, int *nb_colliding_axis) {
+  Couple couples[_nb_polygons * 2];
+  int in[_nb_polygons];
+  memset((void *) nb_colliding_axis, 0, _nb_polygons * _nb_polygons * sizeof(int));
+
+  for(int a = 0; a < nb_axis; a++) {
+    scalar_t alpha = M_PI * scalar_t(a) / scalar_t(nb_axis);
+    scalar_t vx = cos(alpha), vy = sin(alpha);
+
+    for(int n = 0; n < _nb_polygons; n++) {
+      scalar_t *x = _polygons[n]->_x, *y = _polygons[n]->_y;
+      scalar_t min = x[0] * vx + y[0] * vy, max = min;
+
+      for(int v = 1; v < _polygons[n]->_nb_vertices; v++) {
+        scalar_t s = x[v] * vx + y[v] * vy;
+        if(s > max) max = s;
+        if(s < min) min = s;
+      }
+
+      couples[2 * n + 0].value = min;
+      couples[2 * n + 0].index = n;
+      couples[2 * n + 1].value = max;
+      couples[2 * n + 1].index = n;
+    }
+
+    qsort((void *) couples, 2 * _nb_polygons, sizeof(Couple), compare_couple);
+
+    int nb_in = 0;
+    memset((void *) in, 0, _nb_polygons * sizeof(int));
+    for(int k = 0; k < 2 * _nb_polygons; k++) {
+      int i = couples[k].index;
+      in[i] = !in[i];
+      if(in[i]) {
+        nb_in++;
+        if(nb_in > 1) {
+          for(int j = 0; j < i; j++)
+            if(j != i && in[j]) nb_colliding_axis[j + i * _nb_polygons]++;
+          for(int j = i+1; j < _nb_polygons; j++)
+            if(j != i && in[j]) nb_colliding_axis[i + j * _nb_polygons]++;
+        }
+      } else nb_in--;
+    }
+  }
+
+  for(int i = 0; i < _nb_polygons; i++) {
+    for(int j = 0; j < i; j++) {
+      if(nb_colliding_axis[j + i * _nb_polygons] > nb_colliding_axis[i + i * _nb_polygons])
+        nb_colliding_axis[i + i * _nb_polygons] = nb_colliding_axis[j + i * _nb_polygons];
+      nb_colliding_axis[i + j * _nb_polygons] = nb_colliding_axis[j + i * _nb_polygons];
+    }
+  }
+}
+
+bool Universe::update(scalar_t dt) {
+  bool result = false;
+  apply_collision_forces(dt);
+  for(int n = 0; n < _nb_polygons; n++) if(_polygons[n]) {
+    _polygons[n]->apply_border_forces(dt, _width, _height);
+    result |= _polygons[n]->update(dt);
+  }
+  return result;
+}
+
+Polygon *Universe::pick_polygon(scalar_t x, scalar_t y) {
+  for(int n = 0; n < _nb_polygons; n++)
+    if(_polygons[n] && _polygons[n]->contain(x, y)) return _polygons[n];
+  return 0;
+}
+
+#ifdef XFIG_SUPPORT
+void Universe::print_xfig(XFigTracer *tracer) {
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->color_xfig(tracer);
+    }
+  }
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->print_xfig(tracer);
+    }
+  }
+}
+#endif
+
+#ifdef X11_SUPPORT
+void Universe::draw(SimpleWindow *window) {
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->draw(window);
+    }
+  }
+
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->draw_contours(window);
+    }
+  }
+}
+#endif
+
+void Universe::draw(Canvas *canvas) {
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->draw(canvas);
+    }
+  }
+
+  for(int n = 0; n < _nb_polygons; n++) {
+    if(_polygons[n]) {
+      _polygons[n]->draw_contours(canvas);
+    }
+  }
+}
+
+void Universe::apply_collision_forces(scalar_t dt) {
+  const int nb_axis = 2;
+  int nb_collision[_nb_polygons * _nb_polygons];
+
+  compute_pseudo_collisions(nb_axis, nb_collision);
+
+  for(int n = 0; n < _nb_polygons; n++) if(_polygons[n])
+    for(int m = 0; m < _nb_polygons; m++)
+      if(m != n && _polygons[m] && nb_collision[n + _nb_polygons * m] == nb_axis)
+        _polygons[n]->apply_collision_forces(dt, m, _polygons[m]);
+}
diff --git a/universe.h b/universe.h
new file mode 100644 (file)
index 0000000..6cb4193
--- /dev/null
@@ -0,0 +1,78 @@
+
+/*
+ *  dyncnn is a deep-learning algorithm for the prediction of
+ *  interacting object dynamics
+ *
+ *  Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
+ *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ *  This file is part of dyncnn.
+ *
+ *  dyncnn is free software: you can redistribute it and/or modify it
+ *  under the terms of the GNU General Public License version 3 as
+ *  published by the Free Software Foundation.
+ *
+ *  dyncnn is distributed in the hope that it will be useful, but
+ *  WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ *  General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with dyncnn.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef UNIVERSE_H
+#define UNIVERSE_H
+
+#include <iostream>
+#include <cmath>
+
+#include "misc.h"
+#include "canvas.h"
+#include "polygon.h"
+
+#ifdef X11_SUPPORT
+#include "simple_window.h"
+#endif
+
+using namespace std;
+
+class Universe {
+  scalar_t _width, _height;
+public:
+  int _nb_max_polygons, _nb_polygons;
+  Polygon **_polygons;
+
+  inline scalar_t width() { return _width; }
+  inline scalar_t height() { return _height; }
+
+  Universe(int nb_max_polygons, scalar_t width, scalar_t height);
+  // The destructor deletes all the added polygons
+  ~Universe();
+
+  void initialize_polygon(Polygon *p);
+  void clear();
+  void add_polygon(Polygon *p);
+  bool collide(Polygon *p);
+
+  // Compute collisions between projections of the polygons on a few
+  // axis to speed up the computation
+  void compute_pseudo_collisions(int nb_axis, int *nb_colliding_axis);
+  void apply_collision_forces(scalar_t dt);
+  bool update(scalar_t dt);
+
+  Polygon *pick_polygon(scalar_t x, scalar_t y);
+
+#ifdef XFIG_SUPPORT
+  void print_xfig(XFigTracer *tracer);
+#endif
+
+#ifdef X11_SUPPORT
+  void draw(SimpleWindow *window);
+#endif
+
+  void draw(Canvas *canvas);
+};
+
+#endif