From 26abc735ca75a78974c7ced091035e8c27ca1c00 Mon Sep 17 00:00:00 2001 From: Francois Fleuret Date: Tue, 28 Feb 2017 16:06:26 +0100 Subject: [PATCH 1/1] Initial commit. --- Makefile | 52 +++++ build.py | 23 +++ canvas.cc | 27 +++ canvas.h | 39 ++++ canvas_cairo.cc | 153 ++++++++++++++ canvas_cairo.h | 56 +++++ flatland.cc | 249 ++++++++++++++++++++++ flatland.h | 14 ++ misc.cc | 31 +++ misc.h | 53 +++++ mylib.c | 63 ++++++ mylib.h | 2 + polygon.cc | 533 ++++++++++++++++++++++++++++++++++++++++++++++++ polygon.h | 128 ++++++++++++ test.py | 19 ++ universe.cc | 165 +++++++++++++++ universe.h | 66 ++++++ 17 files changed, 1673 insertions(+) create mode 100644 Makefile create mode 100755 build.py create mode 100644 canvas.cc create mode 100644 canvas.h create mode 100644 canvas_cairo.cc create mode 100644 canvas_cairo.h create mode 100644 flatland.cc create mode 100644 flatland.h create mode 100644 misc.cc create mode 100644 misc.h create mode 100644 mylib.c create mode 100644 mylib.h create mode 100644 polygon.cc create mode 100644 polygon.h create mode 100755 test.py create mode 100644 universe.cc create mode 100644 universe.h diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..e800949 --- /dev/null +++ b/Makefile @@ -0,0 +1,52 @@ + +# 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 +# +# 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 . + +ifeq ($(DEBUG),yes) + CXXFLAGS = -fPIC -Wall -g -DDEBUG +else + CXXFLAGS = -fPIC -Wall -g -O3 +endif + +CXXFLAGS += -I/usr/include/cairo -DCAIRO_SUPPORT + +LDFLAGS += -lcairo + +all: flatland.so TAGS + +TAGS: *.cc *.h + etags *.cc *.h + +flatland.so: flatland.o misc.o \ + universe.o \ + polygon.o \ + canvas.o canvas_cairo.o + $(CXX) $(CXXFLAGS) $(LDFLAGS) -shared -fPIC -o $@ $^ + +test: test.c flatland.so + $(CC) flatland.so -o test test.c + +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/build.py b/build.py new file mode 100755 index 0000000..4f00b97 --- /dev/null +++ b/build.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python-for-pytorch + +import os +from torch.utils.ffi import create_extension + +this_file = os.path.dirname(__file__) + +print('__file__', __file__) + +ffi = create_extension( + '_ext.mylib', + headers = [ 'mylib.h' ], + sources = [ 'mylib.c' ], + extra_objects = [ '/home/fleuret/sources/python/flatland/flatland.so' ], + libraries = [ ], + library_dirs = [ ], + define_macros = [ ], + relative_to=this_file, + with_cuda = False +) + +if __name__ == '__main__': + ffi.build() diff --git a/canvas.cc b/canvas.cc new file mode 100644 index 0000000..58d4019 --- /dev/null +++ b/canvas.cc @@ -0,0 +1,27 @@ + +/* + * 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 + * + * 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 . + * + */ + +#include "canvas.h" + +Canvas::~Canvas() {} diff --git a/canvas.h b/canvas.h new file mode 100644 index 0000000..55390ba --- /dev/null +++ b/canvas.h @@ -0,0 +1,39 @@ + +/* + * 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 + * + * 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 . + * + */ + +#ifndef CANVAS_H +#define CANVAS_H + +#include "misc.h" + +class Canvas { +public: + virtual ~Canvas(); + virtual void clear() = 0; + 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 index 0000000..aeb7b9f --- /dev/null +++ b/canvas_cairo.cc @@ -0,0 +1,153 @@ + +/* + * 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 + * + * 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 . + * + */ + +#include "canvas_cairo.h" + +#include + +#define MAX(x, y) ((x >= y) ? (x) : (y)) + +CanvasCairo::CanvasCairo(scalar_t scale, int nb_rows, int nb_cols, CanvasCairo **ca) { + _actual_width = 0; + _actual_height = 0; + + for(int i = 0; i < nb_rows; i++) { + int row_height = 0, row_width = 0; + for(int j = 0; j < nb_cols; j++) { + CanvasCairo *this_ca = ca[i * nb_cols + j]; + row_height = MAX(row_height, this_ca->_actual_height); + row_width += this_ca->_actual_width; + } + _actual_width = MAX(_actual_width, row_width); + _actual_height += row_height; + } + + _data = new unsigned char [_actual_width * _actual_height * _depth]; + + int x0, y0 = 0; + for(int i = 0; i < nb_rows; i++) { + x0 = 0; + int row_height = 0; + for(int j = 0; j < nb_cols; j++) { + CanvasCairo *this_ca = ca[i * nb_cols + j]; + for(int y = 0; y < this_ca->_actual_height; y++) { + for(int x = 0; x < this_ca->_actual_width; x++) { + for(int d = 0; d < _depth; d++) { + _data[(x0 + x + _actual_width * (y0 + y))* _depth + d] = + this_ca->_data[(x + this_ca->_actual_width * y)* _depth + d]; + + } + } + } + row_height = MAX(row_height, this_ca->_actual_height); + x0 += this_ca->_actual_width; + } + y0 += row_height; + } + + _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); + + clear(); + // cairo_set_source_rgb(_context_resource, 1.0, 1.0, 1.0); + // cairo_set_line_width (_context_resource, 1.0); +} + +CanvasCairo::CanvasCairo(scalar_t scale, int width, int height) { + _actual_width = int(width * scale); + _actual_height = int(height * scale); + _scale = scale; + + _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_line_width (_context_resource, 1.0); + + clear(); + // cairo_set_source_rgb(_context_resource, 1.0, 1.0, 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::clear() { + cairo_set_source_rgb(_context_resource, 1.0, 1.0, 1.0); + cairo_rectangle(_context_resource, 0, 0, _actual_width / _scale, _actual_height / _scale); + cairo_fill(_context_resource); +} + +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 index 0000000..764d049 --- /dev/null +++ b/canvas_cairo.h @@ -0,0 +1,56 @@ + +/* + * 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 + * + * 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 . + * + */ + +#ifndef CANVAS_CAIRO_H +#define CANVAS_CAIRO_H + +#include "canvas.h" + +#include +#include + +class CanvasCairo : public Canvas { + scalar_t _scale; + cairo_surface_t *_image; + cairo_t *_context_resource; + +public: + const static int _depth = 4; + int _actual_width, _actual_height; + unsigned char *_data; + + CanvasCairo(scalar_t scale, int nb_rows, int nb_cols, CanvasCairo **ca); + CanvasCairo(scalar_t scale, int width, int height); + CanvasCairo(int nb_rows, int nb_cols, CanvasCairo **x); + + virtual ~CanvasCairo(); + + virtual void clear(); + 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 index 0000000..92a5ebb --- /dev/null +++ b/flatland.cc @@ -0,0 +1,249 @@ + +/* + * 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 + * + * 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 . + * + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +#include "misc.h" +#include "universe.h" +#include "canvas_cairo.h" + +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 << " [--dir ] [--seed ]]" << endl; + exit(1); +} + +////////////////////////////////////////////////////////////////////// + +void draw_universe_on_canvas(CanvasCairo *canvas, scalar_t scaling, + Universe *universe) { + canvas->set_line_width(1.0 / scaling); + universe->draw(canvas); +} + +void draw_grabbing_point_on_canvas(CanvasCairo *canvas, scalar_t scaling, + scalar_t xg, scalar_t yg, + scalar_t r, scalar_t g, scalar_t b) { + scalar_t radius = 1/scaling; + int n = 36; + scalar_t xp[n], yp[n]; + for(int k = 0; k < n; k++) { + scalar_t alpha = 2 * M_PI * scalar_t(k) / scalar_t(n); + xp[k] = xg + radius * cos(alpha); + yp[k] = yg + radius * sin(alpha); + } + canvas->set_drawing_color(r, g, b); + canvas->set_line_width(2.0); + canvas->draw_polygon(1, n, xp, yp); +} + +////////////////////////////////////////////////////////////////////// + +extern "C" void fl_generate_sequences(int nb_sequences, + int nb_images_per_sequence, + int width, int height, + unsigned char *output) { + + const scalar_t world_width = width * 8; + const scalar_t world_height = height * 8; + const scalar_t scaling = 0.125; + + const scalar_t dt = 0.1; + const int nb_iterations_per_steps = 5; + + ////////////////////////////////////////////////////////////////////// + + // We will generate images { 0, every_nth, 2 * every_nth, ..., k * every_nth < nb_simulated_frames } + + // The framerate every_nth may be set to smaller value to generate + // nice materials for presentations or papers. + + int every_nth = 16; + int nb_simulated_frames = 1 + (nb_images_per_sequence - 1) * every_nth; + int random_grasp = 1; + int random_shape_size = 0; + int nb_shapes = 10; + // char data_dir[1024] = "/tmp/"; + // int show_grabbing_point = 0; + int skip = -1; + + for(int n = 0; n < nb_sequences; n++) { + + Universe *universe; + Polygon *grabbed_polygon; + + universe = new Universe(nb_shapes, world_width, world_height); + + const int nb_saved_frames = (nb_simulated_frames + every_nth - 1) / every_nth; + if(nb_saved_frames != nb_images_per_sequence) { + cerr << "It makes no sense." << endl; + abort(); + } + + CanvasCairo *canvases[nb_saved_frames * 2]; + + for(int s = 0; s < 2 * nb_saved_frames; s++) { + canvases[s] = new CanvasCairo(scaling, universe->width(), universe->height()); + } + + scalar_t grab_start_x, grab_start_y; + + if(random_grasp) { + grab_start_x = world_width * (0.1 + 0.8 * drand48()); + grab_start_y = world_height * (0.1 + 0.8 * drand48()); + } else { + grab_start_x = world_width * 0.5; + grab_start_y = world_height * 0.75; + } + + 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; + + scalar_t shape_size; + + if(random_shape_size) { + shape_size = 40 + 80 * drand48(); + } else { + shape_size = 80; + } + + do { + scalar_t x[] = { - shape_size * 0.4, + shape_size * 0.4, + + shape_size * 0.4, - shape_size * 0.4 }; + + scalar_t y[] = { - shape_size * 0.6, - shape_size * 0.6, + + shape_size * 0.6, + shape_size * 0.6 }; + + scalar_t delta = shape_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(); + + delete pol; + pol = new Polygon(0.5, 1.0, 1.0, 1.0, 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 = -1; + universe->clear(); + nb_attempts = 0; + } else { + universe->add_polygon(pol); + } + } + + grabbed_polygon = universe->pick_polygon(grab_start_x, grab_start_y); + } while(!grabbed_polygon); + + if(skip < 0 || n >= skip) { + + 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); + + for(int s = 0; s < nb_simulated_frames; s++) { + if(s % every_nth == 0) { + int t = s / every_nth; + // 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); + + // canvases[2 * t + 0]->clear(); + // draw_grabbing_point_on_canvas(canvases[2 * t + 0], scaling, + // xf, yf, 0.0, 0.0, 0.0); + // canvases[2 * t + 1]->clear(); + // draw_universe_on_canvas(canvases[2 * t + 1], scaling, universe); + + canvases[t]->clear(); + draw_universe_on_canvas(canvases[t], scaling, universe); + + // if(show_grabbing_point) { + // draw_grabbing_point_on_canvas(canvases[2 * t + 1], scaling, + // xf, yf, 1.0, 0.0, 0.0); + // } + } + + if(s < nb_simulated_frames - 1) { + // Run the simulation + 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); + } + } + } + + for(int t = 0; t < nb_images_per_sequence; t++) { + memcpy(output + (n * nb_images_per_sequence + t) * 3 * width * height, + canvases[t]->_data, + 3 * width * height); + } + + // char buffer[1024]; + // sprintf(buffer, "%s/%03d/dyn_%06d.png", data_dir, n / 1000, n); + // FILE *file = safe_fopen(buffer, "w"); + // main_canvas.write_png(file); + // fclose(file); + } + + for(int t = 0; t < 2 * nb_saved_frames; t++) { + delete canvases[t]; + } + + delete universe; + } +} diff --git a/flatland.h b/flatland.h new file mode 100644 index 0000000..9a9247a --- /dev/null +++ b/flatland.h @@ -0,0 +1,14 @@ + +#ifdef __cplusplus +extern "C" { +#endif + +void fl_generate_sequences(int nb_sequences, + int nb_images_per_sequence, + int width, int height, + unsigned char *output); + +#ifdef __cplusplus +} +#endif + diff --git a/misc.cc b/misc.cc new file mode 100644 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 + * + * 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 . + * + */ + +#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 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 + * + * 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 . + * + */ + +#ifndef MISC_H +#define MISC_H + +#include +#include + +#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/mylib.c b/mylib.c new file mode 100644 index 0000000..61a461d --- /dev/null +++ b/mylib.c @@ -0,0 +1,63 @@ + +#include + +/* + + Example of FFI extension I started from: + + https://github.com/pytorch/extension-ffi.git + + There is this tutorial + + https://github.com/pytorch/tutorials/blob/master/Creating%20Extensions%20using%20FFI.md + + And TH's Tensor definition are here in my install: + + anaconda3/lib/python3.5/site-packages/torch/lib/include/TH/generic/THTensor.h + + */ + +#include "flatland.h" + +int generate_sequence(THByteTensor *output) { + long nb_sequences = 1; + long nb_images_per_sequence = 5; + long depth = 3; + long width = 64; + long height = 64; + long s; + unsigned char *a, *b; + int c, k, i, j, st0, st1, st2, st3, st4; + + THByteTensor_resize5d(output, nb_sequences, nb_images_per_sequence, depth, height, width); + + st0 = THByteTensor_stride(output, 0); + st1 = THByteTensor_stride(output, 1); + st2 = THByteTensor_stride(output, 2); + st3 = THByteTensor_stride(output, 3); + st4 = THByteTensor_stride(output, 4); + + a = + THByteTensor_storage(output)->data + THByteTensor_storageOffset(output); + + for(s = 0; s < nb_sequences; s++) { + unsigned char result[nb_images_per_sequence * depth * width * height]; + unsigned char *r = result; + fl_generate_sequences(1, nb_images_per_sequence, width, height, result); + for(k = 0; k < nb_images_per_sequence; k++) { + for(c = 0; c < depth; c++) { + for(i = 0; i < height; i++) { + b = a + + s * st0 + k * st1 + c * st2 + i * st3; + for(j = 0; j < width; j++) { + *b = (unsigned char) (*r); + r++; + b += st4; + } + } + } + } + } + + return 1; +} diff --git a/mylib.h b/mylib.h new file mode 100644 index 0000000..9c9208b --- /dev/null +++ b/mylib.h @@ -0,0 +1,2 @@ + +int generate_sequence(THByteTensor *output); diff --git a/polygon.cc b/polygon.cc new file mode 100644 index 0000000..69b125e --- /dev/null +++ b/polygon.cc @@ -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 + * + * 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 . + * + */ + +#include + +using namespace std; + +#include +#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 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 + * + * 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 . + * + */ + +#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/test.py b/test.py new file mode 100755 index 0000000..d6d2df1 --- /dev/null +++ b/test.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python-for-pytorch + +import torch +import torchvision +from torchvision import datasets + +from _ext import mylib + +x = torch.ByteTensor(4, 5).fill_(0) + +print(x.size()) + +mylib.generate_sequence(x) + +print(x.size()) + +x = x.float().sub_(128).div_(128) + +torchvision.utils.save_image(x[0], 'example.png') diff --git a/universe.cc b/universe.cc new file mode 100644 index 0000000..2b1383d --- /dev/null +++ b/universe.cc @@ -0,0 +1,165 @@ + +/* + * 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 + * + * 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 . + * + */ + +#include + +#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; +} + +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 index 0000000..11dacc7 --- /dev/null +++ b/universe.h @@ -0,0 +1,66 @@ + +/* + * 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 + * + * 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 . + * + */ + +#ifndef UNIVERSE_H +#define UNIVERSE_H + +#include +#include + +#include "misc.h" +#include "canvas.h" +#include "polygon.h" + +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); + + void draw(Canvas *canvas); +}; + +#endif -- 2.39.5