Version 0.2
authorAndreas Eversberg <jolly@eversberg.eu>
Wed, 30 Apr 2014 09:03:03 +0000 (11:03 +0200)
committerAndreas Eversberg <jolly@eversberg.eu>
Sun, 15 Mar 2015 14:33:50 +0000 (15:33 +0100)
- various fixes, cleanups and improvements

18 files changed:
gui/colorize.c
gui/flow.c
gui/image.c
gui/image.h
gui/main.c
gui/main.h
gui/timeline.c
lib/Makefile.am
lib/colorize.c [new file with mode: 0644]
lib/colorize.h
lib/darray.c
lib/darray.h
lib/multigrid.c
lib/multigrid.h [new file with mode: 0644]
src/colorize.c
src/opticalflow.c
src/opticalflow.h
version.h

index 2e1edf2..61d396d 100644 (file)
 #include "../src/mark.h"
 #include "../src/ppm.h"
 #include "../src/yuv.h"
+#include "../src/dir_seperator.h"
 #include "../lib/darray.h"
 #include "../lib/colorize.h"
 #include "palette.h"
+#include "timeline.h"
 
 #define min(x,y) ((x < y) ? x : y)
 
@@ -20,6 +22,8 @@ struct colorize_priv {
        GtkWidget *window;
        GdkPixbuf *pixbuf;
        int width, height;
+       char filename[256];
+       char folder[256];
 };
 
 static void colorize_destroy(GtkWidget *widget, gpointer priv)
@@ -63,6 +67,8 @@ static void save_event(gpointer *priv)
                GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
                NULL);
 
+       gtk_file_chooser_set_current_folder(GTK_FILE_CHOOSER(dialog), cp->folder);
+       gtk_file_chooser_set_current_name(GTK_FILE_CHOOSER(dialog), cp->filename);
        if (gtk_dialog_run(GTK_DIALOG(dialog)) == GTK_RESPONSE_ACCEPT) {
                filename = strdup(gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(dialog)));
        }
@@ -100,13 +106,14 @@ static void close_event(gpointer *priv)
 
 void colorize_image(void)
 {
-       darray_t *gI = NULL, *cI = NULL, *markIm = NULL, *ntscIm = NULL, *resultIm = NULL;
+       darray_t *gI = NULL, *cI = NULL, *markIm = NULL, *ntscIm = NULL;
        double *ptr, *ptr2, *ptr3;
        int dims[4];
        int width = img_width, height = img_height;
        int features, change_bc;
        int i, j;
        int c;
+       int rc;
        GdkPixbuf *pixbuf;
        GtkWidget *image;
        guchar *data;
@@ -119,6 +126,7 @@ void colorize_image(void)
        struct colorize_priv *cp;
        int rs;
        int zoom_field;
+       char *p, *q;
 
        if (!img_pixbuf)
                return;
@@ -240,16 +248,16 @@ void colorize_image(void)
        cI = NULL;
 
        /* render u and v change */
-       resultIm = colorize(ntscIm, markIm, NULL, NULL, NULL, NULL, 5, 1, 0);
-       if (!resultIm) {
+       rc = colorize(ntscIm, markIm, NULL, NULL, 5, 1, 0);
+       if (rc < 0) {
                printerror("No memory! Use smaller image or add more memory.");
                goto error;
        }
 
-       /* if we have a change, we apply brightness+contrast from resultIm */
+       /* if we have a change, we apply brightness+contrast from ntscIm */
        if (change_bc) {
-               ptr = darrayGetPr(resultIm);
-               ptr2 = darrayGetPr(resultIm);
+               ptr = darrayGetPr(ntscIm);
+               ptr2 = darrayGetPr(ntscIm);
                for (i = 0; i < img_height; i++) {
                        for (j = 0; j < img_width; j++) {
                                /* apply contrast */
@@ -271,7 +279,7 @@ void colorize_image(void)
 
        /* apply YUV to pixbuffer and display */
        pixbuf = gdk_pixbuf_new(GDK_COLORSPACE_RGB, FALSE, 8, img_width, img_height * zoom_field);
-       ptr = darrayGetPr(resultIm);
+       ptr = darrayGetPr(ntscIm);
        yuv2rgb(ptr, ptr, width, height);
        data = gdk_pixbuf_get_pixels(pixbuf);
        rs = gdk_pixbuf_get_rowstride(pixbuf);
@@ -302,14 +310,15 @@ void colorize_image(void)
        image = gtk_image_new_from_pixbuf(pixbuf);
        gtk_widget_show(image);
 
-       darrayDestroy(resultIm);
        darrayDestroy(markIm);
+       markIm = NULL;
        darrayDestroy(ntscIm);
+       ntscIm = NULL;
 
        darrayDone();
 
        /* private structure */
-       cp = malloc(sizeof(struct colorize_priv));
+       cp = calloc(sizeof(struct colorize_priv), 1);
 
        /* create window */
        colorize_window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
@@ -321,6 +330,12 @@ void colorize_image(void)
        cp->pixbuf = pixbuf;
        cp->width = img_width*zoom_field;
        cp->height = img_height*zoom_field;
+       p = frame_list[timeline_selected].filename;
+       while((q = strchr(p, DIR_SEPERATOR)))
+               p = q + 1;
+       strcpy(cp->folder, frame_list[timeline_selected].filename);
+       cp->folder[p - frame_list[timeline_selected].filename] = '\0';
+       sprintf(cp->filename, "colorized_%s", p);
 
        /* create menu*/
        menu_bar = gtk_menu_bar_new();
@@ -329,7 +344,7 @@ void colorize_image(void)
        root_menu = gtk_menu_item_new_with_mnemonic("_File");
        gtk_widget_show(root_menu);
 
-       menu_item = gtk_image_menu_item_new_from_stock(GTK_STOCK_SAVE, NULL);
+       menu_item = gtk_image_menu_item_new_from_stock(GTK_STOCK_SAVE_AS, NULL);
        gtk_menu_shell_append(GTK_MENU_SHELL(menu), menu_item);
        g_signal_connect_swapped(menu_item, "activate", G_CALLBACK(save_event), cp);
        gtk_widget_show(menu_item);
@@ -353,12 +368,9 @@ void colorize_image(void)
        gtk_container_add(GTK_CONTAINER(colorize_window), vbox);
        gtk_widget_show(colorize_window);
 
-       return;
-
 error:
        darrayDestroy(gI);
        darrayDestroy(cI);
-       darrayDestroy(resultIm);
        darrayDestroy(markIm);
        darrayDestroy(ntscIm);
 
index e1f6740..2d70331 100644 (file)
 #include "timeline.h"
 #include "image.h"
 
-static GtkWidget *flow_window = NULL;
+static GtkWidget *flow_dialog = NULL;
 static GtkWidget *enable_button;
-static GtkObject *steps_adj;
-static GtkObject *win_adj;
-static GtkObject *max_adj;
+static GtkObject *window_adj;
+static GtkObject *vector_adj;
+static GtkObject *uv_adj;
 
 static void flow_destroy(GtkWidget *widget, gpointer priv)
 {
-       gtk_widget_destroy(flow_window);
-       flow_window = NULL;
+       gtk_widget_destroy(flow_dialog);
+       flow_dialog = NULL;
 }
 
 static gboolean flow_key_press(GtkWidget *widget, GdkEventKey *event, gpointer priv)
@@ -40,8 +40,8 @@ static gboolean flow_key_press(GtkWidget *widget, GdkEventKey *event, gpointer p
 
 static void close_event(gpointer *priv)
 {
-       gtk_widget_destroy(flow_window);
-       flow_window = NULL;
+       gtk_widget_destroy(flow_dialog);
+       flow_dialog = NULL;
 }
 
 static void reset_event(gpointer *priv)
@@ -49,9 +49,9 @@ static void reset_event(gpointer *priv)
        anything_modified = 1;
 
        flow_default();
-       gtk_adjustment_set_value(GTK_ADJUSTMENT(steps_adj), flow_steps);
-       gtk_adjustment_set_value(GTK_ADJUSTMENT(win_adj), flow_win);
-       gtk_adjustment_set_value(GTK_ADJUSTMENT(max_adj), flow_max);
+       gtk_adjustment_set_value(GTK_ADJUSTMENT(window_adj), flow_window);
+       gtk_adjustment_set_value(GTK_ADJUSTMENT(vector_adj), flow_view_vector);
+       gtk_adjustment_set_value(GTK_ADJUSTMENT(uv_adj), flow_view_uv);
        if (flowview)
                timeline_select_and_save(timeline_selected, timeline_selected);
 //             draw_image(0, 0, -1, -1);
@@ -62,9 +62,9 @@ static void apply_event(gpointer *priv)
        anything_modified = 1;
 
        flow_enable = GTK_TOGGLE_BUTTON(enable_button)->active;
-       flow_steps = GTK_ADJUSTMENT(steps_adj)->value;
-       flow_win = GTK_ADJUSTMENT(win_adj)->value;
-       flow_max = GTK_ADJUSTMENT(max_adj)->value;
+       flow_window = GTK_ADJUSTMENT(window_adj)->value;
+       flow_view_vector = GTK_ADJUSTMENT(vector_adj)->value;
+       flow_view_uv = GTK_ADJUSTMENT(uv_adj)->value;
        if (flowview)
                timeline_select_and_save(timeline_selected, timeline_selected);
 //             draw_image(0, 0, -1, -1);
@@ -105,19 +105,19 @@ void flow_event(gpointer priv)
        GtkWidget *menu;
        GtkWidget *menu_item;
        GtkWidget *enable_label;
-       GtkWidget *steps_label;
-       GtkWidget *steps_scroll;
-       GtkWidget *win_label;
-       GtkWidget *win_scroll;
-       GtkWidget *max_label;
-       GtkWidget *max_scroll;
+       GtkWidget *window_label;
+       GtkWidget *window_scroll;
+       GtkWidget *vector_label;
+       GtkWidget *vector_scroll;
+       GtkWidget *uv_label;
+       GtkWidget *uv_scroll;
        GtkWidget *button_apply;
        GtkWidget *button_reset;
        GtkWidget *box;
 
-       if (flow_window) {
-               gtk_widget_destroy(flow_window);
-               flow_window = NULL;
+       if (flow_dialog) {
+               gtk_widget_destroy(flow_dialog);
+               flow_dialog = NULL;
                return;
        }
 
@@ -125,11 +125,11 @@ void flow_event(gpointer priv)
                return;
 
        /* create window */
-       flow_window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-       gtk_window_set_default_size(GTK_WINDOW(flow_window), 300, 0);
-       gtk_window_set_title(GTK_WINDOW(flow_window), "Optical Flow");
-       g_signal_connect(flow_window, "destroy", G_CALLBACK(flow_destroy), flow_window);
-       g_signal_connect(flow_window, "key_press_event", G_CALLBACK(flow_key_press), flow_window);
+       flow_dialog = gtk_window_new(GTK_WINDOW_TOPLEVEL);
+       gtk_window_set_default_size(GTK_WINDOW(flow_dialog), 300, 0);
+       gtk_window_set_title(GTK_WINDOW(flow_dialog), "Optical Flow");
+       g_signal_connect(flow_dialog, "destroy", G_CALLBACK(flow_destroy), flow_dialog);
+       g_signal_connect(flow_dialog, "key_press_event", G_CALLBACK(flow_key_press), flow_dialog);
 
        /* create menu*/
        menu_bar = gtk_menu_bar_new();
@@ -140,7 +140,7 @@ void flow_event(gpointer priv)
 
        menu_item = gtk_image_menu_item_new_from_stock(GTK_STOCK_CLOSE, NULL);
        gtk_menu_shell_append(GTK_MENU_SHELL(menu), menu_item);
-       g_signal_connect_swapped(menu_item, "activate", G_CALLBACK(close_event), flow_window);
+       g_signal_connect_swapped(menu_item, "activate", G_CALLBACK(close_event), flow_dialog);
        gtk_widget_show(menu_item);
        gtk_menu_item_set_submenu(GTK_MENU_ITEM(root_menu), menu);
        gtk_menu_shell_append(GTK_MENU_SHELL(menu_bar), root_menu);
@@ -153,23 +153,23 @@ void flow_event(gpointer priv)
 //     g_signal_connect(enable_button, "toggled", G_CALLBACK(enable_event), NULL);
        gtk_widget_show(enable_button);
 
-       steps_label = gtk_label_new("Steps:"); 
-       gtk_widget_show(GTK_WIDGET(steps_label));
-       steps_adj = gtk_adjustment_new(flow_steps, 1, 31, 1, 1, 1);
-       steps_scroll = gtk_hscale_new(GTK_ADJUSTMENT(steps_adj));
-       gtk_widget_show(steps_scroll);
+       window_label = gtk_label_new("Search Window:"); 
+       gtk_widget_show(GTK_WIDGET(window_label));
+       window_adj = gtk_adjustment_new(flow_window, 3, 61, 1, 1, 1);
+       window_scroll = gtk_hscale_new(GTK_ADJUSTMENT(window_adj));
+       gtk_widget_show(window_scroll);
 
-       win_label = gtk_label_new("Search Window:"); 
-       gtk_widget_show(GTK_WIDGET(win_label));
-       win_adj = gtk_adjustment_new(flow_win, 3, 61, 1, 1, 1);
-       win_scroll = gtk_hscale_new(GTK_ADJUSTMENT(win_adj));
-       gtk_widget_show(win_scroll);
+       vector_label = gtk_label_new("Display vectors (1=disable):"); 
+       gtk_widget_show(GTK_WIDGET(vector_label));
+       vector_adj = gtk_adjustment_new(flow_view_vector, 1, 51, 1, 1, 1);
+       vector_scroll = gtk_hscale_new(GTK_ADJUSTMENT(vector_adj));
+       gtk_widget_show(vector_scroll);
 
-       max_label = gtk_label_new("Maximum flow vector length:"); 
-       gtk_widget_show(GTK_WIDGET(max_label));
-       max_adj = gtk_adjustment_new(flow_max, 0, 61, 1, 1, 1);
-       max_scroll = gtk_hscale_new(GTK_ADJUSTMENT(max_adj));
-       gtk_widget_show(max_scroll);
+       uv_label = gtk_label_new("Display next as UV (0=disable):"); 
+       gtk_widget_show(GTK_WIDGET(uv_label));
+       uv_adj = gtk_adjustment_new(flow_view_uv, 0, 101, 1, 1, 1);
+       uv_scroll = gtk_hscale_new(GTK_ADJUSTMENT(uv_adj));
+       gtk_widget_show(uv_scroll);
 
        /* buttons */
        button_reset = gtk_button_new();
@@ -196,17 +196,17 @@ void flow_event(gpointer priv)
        gtk_box_pack_start(GTK_BOX(vbox), menu_bar, FALSE, FALSE, 2);
        gtk_box_pack_start(GTK_BOX(vbox), enable_label, FALSE, FALSE, 2);
        gtk_box_pack_start(GTK_BOX(vbox), enable_button, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), steps_label, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), steps_scroll, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), win_label, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), win_scroll, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), max_label, FALSE, FALSE, 2);
-       gtk_box_pack_start(GTK_BOX(vbox), max_scroll, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), window_label, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), window_scroll, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), vector_label, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), vector_scroll, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), uv_label, FALSE, FALSE, 2);
+       gtk_box_pack_start(GTK_BOX(vbox), uv_scroll, FALSE, FALSE, 2);
        gtk_box_pack_start(GTK_BOX(vbox), hbox, FALSE, FALSE, 20);
        gtk_widget_show(vbox);
 
-       gtk_container_add(GTK_CONTAINER(flow_window), vbox);
-       gtk_widget_show(flow_window);
+       gtk_container_add(GTK_CONTAINER(flow_dialog), vbox);
+       gtk_widget_show(flow_dialog);
 }
 
 
index 70f5cb7..1590279 100644 (file)
@@ -78,11 +78,61 @@ void create_image(const char *filename)
 
 no_mark:
        /* create pixbuf */
-       create_or_reset_pixbuf();
+       create_or_reset_pixbuf(1);
 }
 
-void create_or_reset_pixbuf(void)
+void create_or_reset_pixbuf(int resize)
 {
+       int w, h, dw = 0, dh = 0;
+       GdkScreen *screen = gtk_window_get_screen(GTK_WINDOW(main_window));
+       GdkRectangle max;
+
+       /* calculate the maximum windows size, so that the border will not exceed the screen size */
+       gdk_window_get_frame_extents(gtk_widget_get_window(main_window), &max);
+       max.width -= main_window->allocation.width;
+       max.height -= main_window->allocation.height;
+       if (max.width < 0)
+               max.width = 0;
+       if (max.height < 0)
+               max.height = 0;
+       max.width = gdk_screen_get_width(screen)-max.width;
+       max.height = gdk_screen_get_height(screen)-max.height;
+
+       if (resize && screen) {
+try_smaller:
+               w = img_scroll->allocation.width-20;
+               h = img_scroll->allocation.height-20;
+               if (w < img_width*img_scale_x/16) {
+                       w = main_window->allocation.width - w + img_width*img_scale_x/16;
+                       if (w > max.width) {
+                               if (img_scale_x > 4) {
+                                       int aspect = img_scale_y / img_scale_x;
+                                       img_scale_x -= 4;
+                                       img_scale_y = img_scale_x * aspect;
+                                       goto try_smaller;
+                               }
+                               w = max.width;
+                       }
+                       dw = w - main_window->allocation.width;
+               }
+               if (h < img_height*img_scale_y/16) {
+                       h = main_window->allocation.height - h + img_height*img_scale_y/16;
+                       if (h > max.height) {
+                               if (img_scale_x > 4) {
+                                       int aspect = img_scale_y / img_scale_x;
+                                       img_scale_x -= 4;
+                                       img_scale_y = img_scale_x * aspect;
+                                       goto try_smaller;
+                               }
+                               h = max.height;
+                       }
+                       dh = h - main_window->allocation.height;
+               }
+//             printf("%d %d\n", gdk_screen_get_width(screen), gdk_screen_get_height(screen));
+               if (dw || dh)
+                       gtk_window_resize(GTK_WINDOW(main_window), main_window->allocation.width+dw, main_window->allocation.height+dh);
+       }
+
        if (img_pixbuf) {
                g_object_unref(img_pixbuf);
                img_pixbuf = NULL;
@@ -90,7 +140,6 @@ void create_or_reset_pixbuf(void)
        img_pixbuf = gdk_pixbuf_new(GDK_COLORSPACE_RGB, FALSE, 8, img_width*img_scale_x/16, img_height*img_scale_y/16);
 
        gtk_drawing_area_size(GTK_DRAWING_AREA(img_drawing_area), img_width*img_scale_x/16, img_height*img_scale_y/16);
-//     gtk_widget_set_size_request(drawing_area, img_width, img_height);
 
        /* label */
        timeline_show_name(timeline_selected);
index 2dca9b7..59b4b4a 100644 (file)
@@ -7,7 +7,7 @@ extern int img_scale_x, img_scale_y;
 extern int anything_modified;
 
 void create_image(const char *filename);
-void create_or_reset_pixbuf(void);
+void create_or_reset_pixbuf(int resize);
 void destroy_image(void);
 void draw_image(int x, int y, int w, int h);
 void paint_brush(int x, int y, int size, int paint);
index 00b76b4..71cff5f 100644 (file)
@@ -21,6 +21,7 @@
 
 #define CC_APP_TITLE "Colorize GTK"
 GtkWidget *main_window;
+GtkWidget *img_scroll;
 GtkWidget *img_drawing_area = NULL;
 GtkWidget *timeline_drawing_area;
 GtkWidget *timeline_scroll;
@@ -439,6 +440,19 @@ static void open_event(gpointer priv)
        GtkWidget *dialog;
        char *filename = NULL;
        
+       if (anything_modified) {
+               int ret;
+               GtkWidget *dialog = gtk_message_dialog_new (GTK_WINDOW(main_window),
+                       GTK_DIALOG_DESTROY_WITH_PARENT,
+                       GTK_MESSAGE_ERROR,
+                       GTK_BUTTONS_OK_CANCEL,
+                       "Image has been changed, really open new image?");
+               ret = gtk_dialog_run(GTK_DIALOG (dialog));
+               gtk_widget_destroy(dialog);
+               if (ret != GTK_RESPONSE_OK)
+                       return;
+       }
+
        dialog = gtk_file_chooser_dialog_new("Select first grey image",
                GTK_WINDOW(main_window),
                GTK_FILE_CHOOSER_ACTION_OPEN,
@@ -446,6 +460,8 @@ static void open_event(gpointer priv)
                GTK_STOCK_OPEN, GTK_RESPONSE_ACCEPT,
                NULL);
 
+       if (frame_list)
+               gtk_file_chooser_set_filename(GTK_FILE_CHOOSER(dialog), frame_list[timeline_selected].filename);
        if (gtk_dialog_run(GTK_DIALOG(dialog)) == GTK_RESPONSE_ACCEPT) {
                filename = strdup(gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(dialog)));
        }
@@ -453,7 +469,7 @@ static void open_event(gpointer priv)
        gtk_widget_destroy(dialog);
 
        if (filename) {
-               timeline_select_and_save(timeline_selected, timeline_selected);
+// do not save         timeline_select_and_save(timeline_selected, timeline_selected);
                
                create_timeline(filename);
                if (frame_list)
@@ -562,7 +578,7 @@ static void zoom_in_event(gpointer priv)
                int aspect = img_scale_y / img_scale_x;
                img_scale_x += 4;
                img_scale_y = img_scale_x * aspect;
-               create_or_reset_pixbuf();
+               create_or_reset_pixbuf(0);
        }
 }
 
@@ -572,7 +588,7 @@ static void zoom_out_event(gpointer priv)
                int aspect = img_scale_y / img_scale_x;
                img_scale_x -= 4;
                img_scale_y = img_scale_x * aspect;
-               create_or_reset_pixbuf();
+               create_or_reset_pixbuf(0);
        }
 }
 
@@ -580,7 +596,7 @@ static void zoom_100_event(gpointer priv)
 {
        img_scale_x = 16;
        img_scale_y = 16;
-       create_or_reset_pixbuf();
+       create_or_reset_pixbuf(0);
 }
 
 static void zoom_field_event(gpointer priv)
@@ -589,7 +605,7 @@ static void zoom_field_event(gpointer priv)
                img_scale_y = img_scale_x*2;
        else
                img_scale_y = img_scale_x;
-       create_or_reset_pixbuf();
+       create_or_reset_pixbuf(1);
 }
 
 static void show_highlighted_event(gpointer priv)
@@ -998,7 +1014,8 @@ static void set_button_toggel_by_label(enum toggle_label label, gboolean active)
 int main(int argc, char *argv[])
 {
        GtkWidget *vbox, *tool_box;
-       GtkWidget *scroll;
+       GtkWidget *paned;
+       GtkWidget *pal_scroll;
        GtkWidget *menu_bar;
        GtkWidget *separator;
        GtkTreeSelection *selection;
@@ -1015,7 +1032,7 @@ int main(int argc, char *argv[])
 
        main_window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
        gtk_widget_set_size_request(GTK_WIDGET(main_window), 400, 200);
-       gtk_window_set_default_size(GTK_WINDOW(main_window), 800, 600);
+       gtk_window_set_default_size(GTK_WINDOW(main_window), 800, 500);
        gtk_window_set_title(GTK_WINDOW(main_window), CC_APP_TITLE " Version "
 #include "../version.h"
        );
@@ -1088,7 +1105,7 @@ int main(int argc, char *argv[])
        }
 
        /* add paned view to vbox (middle part of vbox) */
-       GtkWidget *paned = gtk_hpaned_new ();
+       paned = gtk_hpaned_new ();
        gtk_widget_show(paned);
        gtk_box_pack_start(GTK_BOX(vbox), paned, TRUE, TRUE, 2);
        gtk_paned_set_position(GTK_PANED (paned), 250);
@@ -1101,10 +1118,10 @@ int main(int argc, char *argv[])
        gtk_widget_set_events(palette_treeview, GDK_ENTER_NOTIFY_MASK | GDK_LEAVE_NOTIFY_MASK);
        gtk_widget_show(palette_treeview);
 
-       scroll = gtk_scrolled_window_new(NULL, NULL);
-       gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(scroll), GTK_POLICY_AUTOMATIC, GTK_POLICY_ALWAYS);
-       gtk_widget_show(scroll);
-       gtk_container_add(GTK_CONTAINER(scroll), palette_treeview);
+       pal_scroll = gtk_scrolled_window_new(NULL, NULL);
+       gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(pal_scroll), GTK_POLICY_AUTOMATIC, GTK_POLICY_ALWAYS);
+       gtk_widget_show(pal_scroll);
+       gtk_container_add(GTK_CONTAINER(pal_scroll), palette_treeview);
 
        palette_renderer = gtk_cell_renderer_text_new();
        g_object_set(G_OBJECT(palette_renderer), "editable", FALSE, NULL);
@@ -1125,7 +1142,7 @@ int main(int argc, char *argv[])
 
        selection = gtk_tree_view_get_selection(GTK_TREE_VIEW(palette_treeview));
        g_signal_connect(selection, "changed", G_CALLBACK(palette_change), NULL);
-       gtk_paned_pack1 (GTK_PANED (paned), scroll, TRUE, TRUE);
+       gtk_paned_pack1 (GTK_PANED (paned), pal_scroll, FALSE, TRUE);
 
        /* add drawing area (right of hbox) */
        img_drawing_area = gtk_drawing_area_new ();
@@ -1142,11 +1159,11 @@ int main(int argc, char *argv[])
                | GDK_POINTER_MOTION_HINT_MASK);
        gtk_widget_show(img_drawing_area);
 
-       scroll = gtk_scrolled_window_new(NULL, NULL);
-       gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(scroll), GTK_POLICY_AUTOMATIC, GTK_POLICY_AUTOMATIC);
-       gtk_widget_show(scroll);
-       gtk_scrolled_window_add_with_viewport(GTK_SCROLLED_WINDOW(scroll), img_drawing_area);
-       gtk_paned_pack2 (GTK_PANED (paned), scroll, TRUE, TRUE);
+       img_scroll = gtk_scrolled_window_new(NULL, NULL);
+       gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(img_scroll), GTK_POLICY_AUTOMATIC, GTK_POLICY_AUTOMATIC);
+       gtk_widget_show(img_scroll);
+       gtk_scrolled_window_add_with_viewport(GTK_SCROLLED_WINDOW(img_scroll), img_drawing_area);
+       gtk_paned_pack2 (GTK_PANED (paned), img_scroll, TRUE, TRUE);
 
 #if 0
        /* add horizontal seperation to vbox */
index 8aaa1b3..58ff8c0 100644 (file)
@@ -2,6 +2,7 @@ void printerror(const char *fmt, ...);
 void main_destroy(void);
 extern int highlight, preview, rendered, flowview;
 extern GtkWidget *main_window;
+extern GtkWidget *img_scroll;
 extern GtkWidget *img_drawing_area;
 extern GtkWidget *timeline_drawing_area;
 extern GtkWidget *palette_treeview;
index 081113e..5a25b86 100644 (file)
@@ -297,9 +297,9 @@ int timeline_select_and_save(int old_frame, int new_frame)
                        (new_frame < timeline_frames - 1) ? frame_list[new_frame + 1].filename : NULL,
                        img_grey_buffer,
                        img_width, img_height,
-                       flow_steps,
-                       flow_win,
-                       flow_max);
+                       flow_window,
+                       flow_view_vector,
+                       flow_view_uv);
 #endif
 
        if (timeline_frames > 1) {
index 187c279..f45e7b0 100644 (file)
@@ -7,5 +7,6 @@ lib_LTLIBRARIES = libcolorize.la
 
 libcolorize_la_SOURCES = \
        darray.c \
-       multigrid.c
+       multigrid.c \
+       colorize.c
 
diff --git a/lib/colorize.c b/lib/colorize.c
new file mode 100644 (file)
index 0000000..6c8caba
--- /dev/null
@@ -0,0 +1,484 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include "darray.h"
+#include "multigrid.h"
+#include "colorize.h"
+
+/* colorize image array using image_mark
+ *
+ * image:
+ *     first dimansion = width
+ *     second dimension = height
+ *     thrid dimension = components (first component is luminance, other components are features, like u, v)
+ *     fourth dimension = time
+ *
+ * image_mark:
+ *     first dimansion = width
+ *     second dimension = height
+ *     fourth dimension = time
+ */
+int colorize_prepare(struct colorize *col)
+{
+       int dimensions[4];
+       int ww, hh, kk, last_kk, d, temp, z;
+       darray_t *luminance = NULL;
+       int rc = -1;
+
+       ASSERT(darrayGetNumberOfDimensions(col->image) == 4, "image array must have 4 dimensions");
+       col->w = darrayGetDimensions(col->image)[0];
+       col->h = darrayGetDimensions(col->image)[1];
+       col->comp = darrayGetDimensions(col->image)[2];
+       col->k = darrayGetDimensions(col->image)[3];
+       ASSERT(col->comp > 0, "image array requires at least one component (Y-component)");
+
+       ASSERT(darrayGetNumberOfDimensions(col->image_mark) == 3, "image_mark array must have 3 dimensions");
+       ASSERT(col->w == darrayGetDimensions(col->image_mark)[0], "image_mark and image must both have same x-size");
+       ASSERT(col->h == darrayGetDimensions(col->image_mark)[1], "image_mark and image must both have same y-size");
+       ASSERT(col->k == darrayGetDimensions(col->image_mark)[2], "image_mark and image must both have same z-size");
+
+       /* Calculate maximum depth
+        * Scale down dimensions until at least 2x2 remains
+        */
+       temp = min(col->w, col->h);
+       for (col->max_depth = 0; temp > 2; col->max_depth++)
+               temp=(temp+1)>>1;
+       ASSERT(col->max_depth >= 2, "image too small");
+
+       col->nb_list = calloc(col->max_depth, sizeof(struct neighbors *));
+       if (!col->nb_list)
+               goto error;
+       col->marks = calloc(col->max_depth, sizeof(darray_t *));
+       if (!col->marks)
+               goto error;
+       col->values = calloc(col->max_depth, sizeof(darray_t *));
+       if (!col->values)
+               goto error;
+       col->flows = calloc(col->max_depth, sizeof(darray_t *));
+       if (!col->flows)
+               goto error;
+       col->flows_i = calloc(col->max_depth, sizeof(darray_t *));
+       if (!col->flows_i)
+               goto error;
+       if (col->scalexyz) {
+               /* the flows_xy are scaled in xy in the current level */
+               col->flows_xy = calloc(col->max_depth, sizeof(darray_t *));
+               if (!col->flows_xy)
+                       goto error;
+       }
+       dimensions[0] = col->w;
+       dimensions[1] = col->h;
+       dimensions[2] = col->k;
+       luminance = darrayCreate(3, dimensions);
+       if (!luminance)
+               goto error;
+       /* generate sparse multigrid */
+       last_kk = kk = col->k;
+       printf("-> prepare sparse grid");
+       for (d = 0, ww=col->w, hh=col->h; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
+               /* create temporary arrays
+                * copy fist array to level 0 or shrink previous level
+                */
+               printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
+               dimensions[0] = ww;
+               dimensions[1] = hh;
+               dimensions[2] = last_kk;
+               col->marks[d] = darrayCreate(3, dimensions);
+               if (!col->marks[d])
+                       goto error;
+               dimensions[0] = ww;
+               dimensions[1] = hh;
+               dimensions[2] = last_kk-1;
+               dimensions[3] = 2;
+               col->flows[d] = darrayCreate(4, dimensions);
+               if (!col->flows[d])
+                       goto error;
+               col->flows_i[d] = darrayCreate(4, dimensions);
+               if (!col->flows_i[d])
+                       goto error;
+               if (d < col->scalexyz) {
+                       dimensions[0] = (ww+1)>>1;
+                       dimensions[1] = (hh+1)>>1;
+                       dimensions[2] = kk-1;
+                       dimensions[3] = 2;
+                       col->flows_xy[d] = darrayCreate(4, dimensions);
+                       if (!col->flows_xy[d])
+                               goto error;
+               }
+               if (d == 0) {
+                       for (z = 0; z < kk; z++) {
+                               darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
+                               darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
+                               memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
+                       }
+                       if (col->flow && col->flow_i) {
+                               darrayCheckPr(col->flows[d], darrayGetPr(col->flows[d]), col->w*col->h*(kk-1)*2);
+                               darrayCheckPr(col->flow, darrayGetPr(col->flow), col->w*col->h*(kk-1)*2);
+                               memcpy(darrayGetPr(col->flows[d]), darrayGetPr(col->flow), col->w*col->h*(kk-1)*2*sizeof(double));
+                               darrayCheckPr(col->flows_i[d], darrayGetPr(col->flows_i[d]), col->w*col->h*(kk-1)*2);
+                               darrayCheckPr(col->flow_i, darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2);
+                               memcpy(darrayGetPr(col->flows_i[d]), darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2*sizeof(double));
+                       }
+                       darrayCheckPr(col->marks[d], darrayGetPr(col->marks[d]), col->w*col->h*kk);
+                       darrayCheckPr(col->image_mark, darrayGetPr(col->image_mark), col->w*col->h*kk);
+                       memcpy(darrayGetPr(col->marks[d]), darrayGetPr(col->image_mark), col->w*col->h*kk*sizeof(double));
+                       /* invert, so 0 = marked */
+                       quantize(col->marks[d], 0.5, 0, 1);
+               } else {
+                       rc = shrink_xy(luminance, NULL);
+                       if (rc < 0)
+                               goto error;
+                       rc = shrink_xy(col->marks[d-1], col->marks[d]);
+                       if (rc < 0)
+                               goto error;
+                       rc = shrink_xy(col->flows[d-1], col->flows[d]);
+                       if (rc < 0)
+                               goto error;
+                       rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
+                       if (rc < 0)
+                               goto error;
+                       /* shrink length of flow vector on smaller grid.
+                        * This is not performed when comparing with
+                        * the original algorithm. (it is a bug there)
+                        */
+                       multiply(col->flows[d], 0.5);
+                       multiply(col->flows_i[d], 0.5);
+                       if (d-1 < col->scalexyz) { /* process if previous d shall be shrinked in z direction */
+                               /* shrink time */
+                               rc = shrink_z(luminance, NULL, col->flows[d]);
+                               if (rc < 0)
+                                       goto error;
+                               rc = shrink_z(col->marks[d], NULL, col->flows[d]);
+                               if (rc < 0)
+                                       goto error;
+                               rc = shrink_flow_z(col->flows[d], NULL);
+                               if (rc < 0)
+                                       goto error;
+                               rc = shrink_flow_z(col->flows_i[d], NULL);
+                               if (rc < 0)
+                                       goto error;
+                       }
+                       /* If any of the restricted pixles is 0 (marked),
+                        * the result must be 0 (marked) too.
+                        */
+                       quantize(col->marks[d], 0.9999, 1, 0);
+               }
+               if (d < col->scalexyz) {
+                       /* create x and y shrinked flow vectors for restriction and prolonging process */
+                       rc = shrink_xy(col->flows[d], col->flows_xy[d]);
+                       if (rc < 0)
+                               goto error;
+                       multiply(col->flows_xy[d], 0.5);
+               }
+
+               /* generate values array */
+               dimensions[0] = ww;
+               dimensions[1] = hh;
+               dimensions[2] = kk;
+               col->values[d] = darrayCreate(3, dimensions);
+               if (!col->values[d])
+                       goto error;
+
+               /* generate array of neighbors */
+#ifndef TEST_NO_RENDER
+               col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
+               if (!col->nb_list[d])
+                       goto error;
+               /* weighten */
+               weighten(luminance, col->nb_list[d]);
+#endif
+
+               last_kk = kk;
+               if (d < col->scalexyz)
+                       kk = (kk+1)>>1;
+       }
+       darrayDestroy(luminance);
+       luminance = NULL;
+       printf("\n");
+
+       dimensions[0] = col->w;
+       dimensions[1] = col->h;
+       dimensions[2] = col->k;
+       col->init = darrayCreate(3, dimensions);
+       if (!col->init)
+               goto error;
+
+#ifdef TEST_PYRAMID
+       dimensions[0] = col->w;
+       dimensions[1] = col->h;
+       dimensions[2] = col->k;
+       luminance = darrayCreate(3, dimensions);
+       for (z = 0; z < col->k; z++) {
+               darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
+               darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
+               memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
+               darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+               darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
+               memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
+#ifdef TEST_FLOW
+               draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
+#endif
+#ifdef TEST_MARKED
+               paste(col->marks[0], col->values[0], 0, 0, 0);
+#endif
+
+       }
+       for (d = 1; d < col->max_depth; d++) {
+               rc = shrink_xy(luminance, NULL);
+               if (rc < 0)
+                       goto error;
+               if (d-1 < col->scalexyz)
+                       shrink_z(luminance, NULL, col->flows_xy[d-1]);
+               paste(luminance, col->values[0], 0, 0, 0);
+#ifdef TEST_FLOW
+               draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
+#endif
+#ifdef TEST_MARKED
+               paste(col->marks[d], col->values[0], 0, 0, 0);
+#endif
+       }
+       darrayDestroy(luminance);
+       luminance = NULL;
+       /* copy test pyramid image as Y-component */
+       for (z = 0; z < col->k; z++) {
+               darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
+               darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+               memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+       }
+#endif
+
+       rc = 0;
+
+error:
+       darrayDestroy(luminance);
+
+       return rc;
+}
+
+int colorize_solve(struct colorize *col, int quick)
+{
+       int z, d, c, ww, hh, kk;
+       int rc = -1;
+
+       if (!quick)
+               printf("-> apply component to grid and solve");
+       for (c = 1; c < col->comp; c++) {
+               int pixles, j;
+               double *value_ptr, *init_ptr, *mark_ptr;
+
+               printf(" #%d", c); fflush(stdout);
+               /* apply component
+                * - copy component into value array (level 0), otherwise 0
+                * - use maked colors for init array, otherwise use 0
+                */
+               for (z = 0; z < col->k; z++) {
+                       darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+                       darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
+                       memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
+               }
+               for (d = 1, ww=(col->w+1)>>1, hh=(col->h+1)>>1, kk=col->k; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
+                       if (d-1 < col->scalexyz) /* scale kk, if last array was scaled */
+                               kk=(kk+1)>>1;   
+                       value_ptr = darrayGetPr(col->values[d]);
+                       pixles = ww*hh*kk;
+                       for (j = 0; j < pixles; j++) {
+                               darrayCheckPr(col->values[d], value_ptr, 1);
+                               *value_ptr++ = 0.0;
+                       }
+               }
+               value_ptr = darrayGetPr(col->values[0]);
+               init_ptr = darrayGetPr(col->init);
+               mark_ptr = darrayGetPr(col->image_mark);
+               pixles = col->w*col->h*col->k;
+               for (j = 0; j < pixles; j++) {
+                       darrayCheckPr(col->image_mark, mark_ptr, 1);
+                       if (*mark_ptr++ > 0.5) {
+                               darrayCheckPr(col->init, init_ptr, 1);
+                               darrayCheckPr(col->values[0], value_ptr, 1);
+                               *init_ptr++ = *value_ptr;
+                       } else {
+                               darrayCheckPr(col->init, init_ptr, 1);
+                               *init_ptr++ = 0.0;
+                       }
+                       value_ptr++;
+               }
+
+#ifndef TEST_NO_RENDER
+               /* colorize component */
+               if (quick) {
+                       rc = solve_mg(col->max_depth, col->inner_iter, 1, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+                       if (rc < 0)
+                               goto error;
+               } else {
+                       int i;
+                       for (i = 0; i < col->outer_iter; i++) {
+                               rc = solve_mg(col->max_depth, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+                               if (rc < 0)
+                                       goto error;
+                               rc = solve_mg((col->max_depth+1) / 2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+                               if (rc < 0)
+                                       goto error;
+                               rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+                               if (rc < 0)
+                                       goto error;
+                               rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+                               if (rc < 0)
+                                       goto error;
+                       }
+               }
+#endif
+
+#ifdef TEST_PYRAMID
+               for (d = 1; d < col->max_depth; d++)
+                       paste(col->values[d], col->values[0], 0, 0, 0);
+#endif
+
+               /* apply colorized component */
+               for (z = 0; z < col->k; z++) {
+                       darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
+                       darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+                       memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+               }
+       }
+       if (!quick)
+               printf("\n");
+
+#ifdef TEST_SHRINK_XY
+       for (c = 0; c < col->comp; c++) {
+               darray_t *temp;
+               int dimensions[3] = {col->w, col->h, col->k};
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+               dimensions[0] = (col->w+1)>>1;
+               dimensions[1] = (col->h+1)>>1;
+               temp = darrayCreate(3, dimensions);
+               ASSERT(temp, "no memory");
+               rc = shrink_xy(col->values[0], temp);
+               ASSERT(rc >= 0, "no memory");
+               paste(temp, col->values[0], 0, 0, 0);
+               darrayDestroy(temp);
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+       }
+#endif
+
+#ifdef TEST_SHRINK_Z
+       for (c = 0; c < col->comp; c++) {
+               darray_t *temp;
+               int dimensions[3] = {col->w, col->h, col->k};
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+               dimensions[2] = (col->k+1)>>1;
+               temp = darrayCreate(3, dimensions);
+               ASSERT(temp, "no memory");
+               rc = shrink_z(col->values[0], temp, col->flows[0]);
+               ASSERT(rc >= 0, "no memory");
+               paste(temp, col->values[0], 0, 0, 0);
+               darrayDestroy(temp);
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+       }
+#endif
+
+#ifdef TEST_EXPAND_XY
+       for (c = 0; c < col->comp; c++) {
+               darray_t *temp;
+               int dimensions[3] = {col->w, col->h, col->k};
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+               dimensions[0] = (col->w+1)>>1;
+               dimensions[1] = (col->h+1)>>1;
+               temp = darrayCreate(3, dimensions);
+               ASSERT(temp, "no memory");
+               rc = shrink_xy(col->values[0], temp);
+               ASSERT(rc >= 0, "no memory");
+               expand_xy(temp, col->values[0]);
+               darrayDestroy(temp);
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+       }
+#endif
+
+#ifdef TEST_EXPAND_Z
+       for (c = 0; c < col->comp; c++) {
+               darray_t *temp;
+               int dimensions[3] = {col->w, col->h, col->k};
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+               dimensions[2] = (col->k+1)>>1;
+               temp = darrayCreate(3, dimensions);
+               ASSERT(temp, "no memory");
+               rc = shrink_z(col->values[0], temp, col->flows[0]);
+               ASSERT(rc >= 0, "no memory");
+               rc = expand_z(temp, col->values[0], col->flows[0]);
+               ASSERT(rc >= 0, "no memory");
+               darrayDestroy(temp);
+               for (z = 0; z < col->k; z++)
+                       memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+       }
+#endif
+
+       rc = 0;
+
+error:
+       return rc;
+}
+
+void colorize_free(struct colorize *col)
+{
+       int d;
+
+       /* free memory */
+       for (d = 0; d < col->max_depth; d++) {
+               if (col->nb_list)
+                       free(col->nb_list[d]);
+               if (col->marks)
+                       darrayDestroy(col->marks[d]);
+               if (col->values)
+                       darrayDestroy(col->values[d]);
+               if (col->flows_xy)
+                       darrayDestroy(col->flows_xy[d]);
+               if (col->flows)
+                       darrayDestroy(col->flows[d]);
+               if (col->flows_i)
+                       darrayDestroy(col->flows_i[d]);
+       }
+       free(col->nb_list);
+       free(col->marks);
+       free(col->values);
+       free(col->flows_xy);
+       free(col->flows);
+       free(col->flows_i);
+       darrayDestroy(col->init);
+}
+
+int colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow, const darray_t *flow_i, int inner_iter, int outer_iter, int scalexyz)
+{
+       struct colorize col;
+       int rc = -1;
+
+       memset(&col, 0, sizeof(col));
+       col.image = image;
+       col.image_mark = image_mark;
+       col.flow = flow;
+       col.flow_i = flow_i;
+       col.inner_iter = inner_iter;
+       col.outer_iter = outer_iter;
+       col.scalexyz = scalexyz;
+
+       rc = colorize_prepare(&col);
+       if (rc < 0)
+               goto error;
+
+       rc = colorize_solve(&col, 0);
+       if (rc < 0)
+               goto error;
+
+       rc = 0;
+
+error:
+       colorize_free(&col);
+
+       return rc;
+}
index 228e3d6..295bf36 100644 (file)
@@ -1,4 +1,30 @@
 #ifndef ___COLORIZE_H___
 #define ___COLORIZE_H___
-darray_t *colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy, int inner_iter, int outer_iter, int scalexyz);
+
+struct colorize {
+       /* input stuff */
+       const darray_t *image;
+       const darray_t *image_mark;
+       const darray_t *flow;
+       const darray_t *flow_i;
+       int inner_iter;
+       int outer_iter;
+       int scalexyz;
+
+       /* internal stuff */
+       int w, h, k, comp, max_depth;
+       unsigned char **nb_list;
+       darray_t **values;
+       darray_t **marks;
+       darray_t *init;
+       darray_t **flows;
+       darray_t **flows_i;
+       darray_t **flows_xy;
+};
+
+int colorize_prepare(struct colorize *col);
+int colorize_solve(struct colorize *col, int quick);
+void colorize_free(struct colorize *col);
+
+int colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow, const darray_t *flow_i, int inner_iter, int outer_iter, int scalexyz);
 #endif
index adf205c..aed5911 100644 (file)
@@ -241,3 +241,17 @@ void darrayCheckList(void)
        }
 }
 
+#ifdef DARRAY_CHECK
+void darrayCheckPr(const darray_t *array, const double *ptr, int num)
+{
+       if (ptr < array->data) {
+               printf("darrayCheckPr FALED: pointer access before array data\n");
+               abort();
+       }
+       if (((unsigned char *)(array->data) + array->size) < ((unsigned char *)ptr + num*sizeof(double))) {
+               printf("darrayCheckPr FALED: pointer access past array data\n");
+               abort();
+       }
+}
+#endif
+
index f328824..c6d4789 100644 (file)
@@ -1,6 +1,12 @@
 #ifndef _ARRAY_H
 #define _ARRAY_H
 
+/* enable check function for access outside array data
+ * this requres darrayCheckPr to be called before each access
+ * when disabled, the darrayCheckPr call is not compiled at all
+ */
+//#define DARRAY_CHECK
+
 // n-dimensional arrays of double type elements
 //
 
@@ -35,4 +41,10 @@ int compare_darray_file(const darray_t *first, const char *file, double threshol
 void darrayDone(void);
 void darrayCheckList(void);
 
+#ifdef DARRAY_CHECK
+void darrayCheckPr(const darray_t *array, const double *ptr, int num);
+#else
+#define darrayCheckPr(array, ptr, num)
+#endif
+
 #endif /* _ARRAY_H */
index 20f3312..1687d07 100644 (file)
@@ -3,7 +3,11 @@
 #include <stdlib.h>
 #include <math.h>
 #include "darray.h"
-#include "colorize.h"
+#include "multigrid.h"
+
+/* if defined, residual is calculated and loop is aborted, if no improvement
+ * is reached. this does not seem to work, so it should be undefined */
+//#define RESIDUAL
 
 /* just process, but do not generate neighbor list and do not solve multigrid */
 //#define TEST_NO_RENDER
 /* test expanding of X, Y and Z */
 //#define TEST_EXPAND_Z
 
-/* print error message b, if a is false */
-#define ASSERT(a, b) \
-{ if (!(a)) { fprintf(stderr, "%s %s() Error in line %d: %s\n", __FILE__, __func__, __LINE__, (b)); abort(); } }
-
-void ASSERT_FLOW(int w, int h, int k, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+void ASSERT_FLOW(int w, int h, int k, const darray_t *flow, const darray_t *flow_i)
 {
-       if (flow_x) {
-               ASSERT(darrayGetNumberOfDimensions(flow_x) == 3, "flow_x array must have 3 dimensions");
-               ASSERT(w == darrayGetDimensions(flow_x)[0], "x-sizes of flow_x vectors are different");
-               ASSERT(h == darrayGetDimensions(flow_x)[1], "y-sizes of flow_x vectors are different");
-               ASSERT(k-1 == darrayGetDimensions(flow_x)[2], "z-sizes of flow_x vectors are not k-1");
-       }
-       if (flow_y) {
-               ASSERT(darrayGetNumberOfDimensions(flow_y) == 3, "flow_y array must have 3 dimensions");
-               ASSERT(w == darrayGetDimensions(flow_y)[0], "x-sizes of flow_y vectors are different");
-               ASSERT(h == darrayGetDimensions(flow_y)[1], "y-sizes of flow_y vectors are different");
-               ASSERT(k-1 == darrayGetDimensions(flow_y)[2], "z-sizes of flow_y vectors are not k-1");
+       if (flow) {
+               ASSERT(darrayGetNumberOfDimensions(flow) == 4, "flow array must have 4 dimensions");
+               ASSERT(w == darrayGetDimensions(flow)[0], "x-sizes of flow vectors are different");
+               ASSERT(h == darrayGetDimensions(flow)[1], "y-sizes of flow vectors are different");
+               ASSERT(k-1 == darrayGetDimensions(flow)[2], "z-sizes of flow vectors are not k-1");
+               ASSERT(2 == darrayGetDimensions(flow)[3], "number of flow cubes are not 2");
        }
-       if (flow_ix) {
-               ASSERT(darrayGetNumberOfDimensions(flow_ix) == 3, "flow_ix array must have 3 dimensions");
-               ASSERT(w == darrayGetDimensions(flow_ix)[0], "x-sizes of flow_ix vectors are different");
-               ASSERT(h == darrayGetDimensions(flow_ix)[1], "y-sizes of flow_ix vectors are different");
-               ASSERT(k-1 == darrayGetDimensions(flow_ix)[2], "z-sizes of flow_ix vectors are not k-1");
-       }
-       if (flow_iy) {
-               ASSERT(darrayGetNumberOfDimensions(flow_iy) == 3, "flow_iy array must have 3 dimensions");
-               ASSERT(w == darrayGetDimensions(flow_iy)[0], "x-sizes of flow_iy vectors are different");
-               ASSERT(h == darrayGetDimensions(flow_iy)[1], "y-sizes of flow_iy vectors are different");
-               ASSERT(k-1 == darrayGetDimensions(flow_iy)[2], "z-sizes of flow_iy vectors are not k-1");
+       if (flow_i) {
+               ASSERT(darrayGetNumberOfDimensions(flow_i) == 4, "flow_i array must have 4 dimensions");
+               ASSERT(w == darrayGetDimensions(flow_i)[0], "x-sizes of flow_i vectors are different");
+               ASSERT(h == darrayGetDimensions(flow_i)[1], "y-sizes of flow_i vectors are different");
+               ASSERT(k-1 == darrayGetDimensions(flow_i)[2], "z-sizes of flow_i vectors are not k-1");
+               ASSERT(2 == darrayGetDimensions(flow_i)[3], "number of flow_i cubes are not 2");
        }
 }
 
-#define min(a, b) ((a < b) ? a : b)
-#define max(a, b) ((a > b) ? a : b)
-
 #define NEIGHBORS_2D 9
 #define NEIGHBORS_3D 27
 
@@ -142,7 +129,7 @@ static struct neigh_off {
  *
  * src plane 4 ------->  dst plane 2
  */
-static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const darray_t *flow_y)
+int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow)
 {
        int x, y, z, xx, yy, w, h, ks, kd;
        int i, pixles;
@@ -170,7 +157,7 @@ static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const
        ASSERT(w == darrayGetDimensions(dst)[0], "x-size of source and destination must be equal");
        ASSERT(h == darrayGetDimensions(dst)[1], "y-size of source and destination must be equal");
        ASSERT(kd == ((ks+1)>>1), "z-size of destination must be half of the source (+an extra plane if source z-size is odd)");
-       ASSERT_FLOW(w, h, ks, flow_x, flow_y, NULL, NULL);
+       ASSERT_FLOW(w, h, ks, flow, NULL);
 
        sum = darrayCreate(3, darrayGetDimensions(dst));
        if (!sum)
@@ -179,16 +166,22 @@ static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const
        /* clear dst array with zeros, fill sum array with ones */
        pixles = w*h*kd;
        dst_ptr = darrayGetPr(dst);
-       for (i = 0; i < pixles; i++)
+       for (i = 0; i < pixles; i++) {
+               darrayCheckPr(dst, dst_ptr, 1);
                *dst_ptr++ = 0.0;
+       }
        sum_ptr = darrayGetPr(sum);
-       for (i = 0; i < pixles; i++)
+       for (i = 0; i < pixles; i++) {
+               darrayCheckPr(sum, sum_ptr, 1);
                *sum_ptr++ = 1.0;
+       }
 
        /* first copy all even planes */
        for (z = 0; z < ks; z++) {
                if ((z&1))
                        continue;
+               darrayCheckPr(src, darrayGetPr(src) + z*w*h, w*h);
+               darrayCheckPr(dst, darrayGetPr(dst) + (z/2)*w*h, w*h);
                memcpy(darrayGetPr(dst) + (z/2)*w*h, darrayGetPr(src) + z*w*h, w*h*sizeof(double));
        }
 
@@ -200,14 +193,19 @@ static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const
                src_ptr = darrayGetPr(src) + z*w*h;
                dst_ptr = darrayGetPr(dst) + (z/2)*w*h;
                sum_ptr = darrayGetPr(sum) + (z/2)*w*h;
-               flow_x_ptr = darrayGetPr(flow_x) + (z-1)*w*h;
-               flow_y_ptr = darrayGetPr(flow_y) + (z-1)*w*h;
+               flow_x_ptr = darrayGetPr(flow) + (z-1)*w*h;
+               flow_y_ptr = flow_x_ptr + (ks-1)*w*h;
                for (y = 0; y < h; y++) {
                        for (x = 0; x < w; x++) {
+                               darrayCheckPr(flow, flow_x_ptr, 1);
                                xx = x + *flow_x_ptr++;
+                               darrayCheckPr(flow, flow_y_ptr, 1);
                                yy = y + *flow_y_ptr++;
                                if (xx>=0 && yy>=0 && xx<w && yy<h) {
+                                       darrayCheckPr(dst, dst_ptr, 1);
+                                       darrayCheckPr(src, src_ptr + xx + yy*w, 1);
                                        *dst_ptr += src_ptr[xx + yy*w];
+                                       darrayCheckPr(sum, sum_ptr, 1);
                                        (*sum_ptr) += 1.0;
                                }
                                dst_ptr++;
@@ -220,8 +218,11 @@ static int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow_x, const
        pixles = w*h*kd;
        dst_ptr = darrayGetPr(dst);
        sum_ptr = darrayGetPr(sum);
-       for (i = 0; i < pixles; i++)
+       for (i = 0; i < pixles; i++) {
+               darrayCheckPr(dst, dst_ptr, 1);
+               darrayCheckPr(sum, sum_ptr, 1);
                *dst_ptr++ /= *sum_ptr++;
+       }
 
        darrayDestroy(sum);
        if (dst_move)
@@ -318,7 +319,7 @@ error:
  * A = 0.707
  * B = 0.293
  */
-static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x, const darray_t *flow_y)
+static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow)
 {
        int x, y, z, xx, yy, xxx, yyy, w, h, ks, kd, k_;
        double _x, _y;
@@ -337,7 +338,7 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
        ASSERT(w == darrayGetDimensions(dst)[0], "x-size of source and destination must be equal");
        ASSERT(h == darrayGetDimensions(dst)[1], "y-size of source and destination must be equal");
        ASSERT(ks == ((kd+1)>>1), "z-size of source must be half of the destination (+an extra plane if destination z-size is odd)");
-       ASSERT_FLOW(w, h, kd, flow_x, flow_y, NULL, NULL);
+       ASSERT_FLOW(w, h, kd, flow, NULL);
 
        sum = darrayCreate(3, darrayGetDimensions(dst));
        if (!sum)
@@ -353,6 +354,8 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
        k_ = kd&(~1); /* even number of planes */
        if (k_) {
                /* copy first plane with weight 1.0 */
+               darrayCheckPr(dst, darrayGetPr(dst) + z*w*h, w*h);
+               darrayCheckPr(src, darrayGetPr(src) + (z/2)*w*h, w*h);
                memcpy(darrayGetPr(dst) + z*w*h, darrayGetPr(src) + (z/2)*w*h, w*h*sizeof(double));
                /* fill sum array with ones */
                pixles = w*h;
@@ -369,37 +372,53 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
                        sum2_ptr = sum1_ptr + w*h;
                        src1_ptr = darrayGetPr(src) + (z/2)*w*h;
                        src2_ptr = src1_ptr + w*h;
-                       flow_x1_ptr = darrayGetPr(flow_x) + z*w*h;
-                       flow_y1_ptr = darrayGetPr(flow_y) + z*w*h;
+                       flow_x1_ptr = darrayGetPr(flow) + z*w*h;
+                       flow_y1_ptr = flow_x1_ptr + (kd-1)*w*h;
                        flow_x2_ptr = flow_x1_ptr + w*h;
                        flow_y2_ptr = flow_y1_ptr + w*h;
                        for (y = 0; y < h; y++) {
                                for (x = 0; x < w; x++) {
                                        /* even plane (weight A) */
+                                       darrayCheckPr(dst, dst2_ptr + x + y*w, 1);
+                                       darrayCheckPr(src, src2_ptr + x + y*w, 1);
                                        dst2_ptr[x + y*w] += src2_ptr[x + y*w] * A;
+                                       darrayCheckPr(sum, sum2_ptr + x + y*w, 1);
                                        sum2_ptr[x + y*w] += A;
                                        /* get position via firt plane */
+                                       darrayCheckPr(flow, flow_x1_ptr, 1);
                                        _x = *flow_x1_ptr++;
+                                       darrayCheckPr(flow, flow_y1_ptr, 1);
                                        _y = *flow_y1_ptr++;
                                        xx = x + round(_x);
                                        yy = y + round(_y);
                                        if (xx<0 || yy<0 || xx>=w || yy>=h)
                                                continue;
                                        /* odd plane (weight A)*/
+                                       darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+                                       darrayCheckPr(src, src1_ptr + x + y*w, 1);
                                        dst1_ptr[xx + yy*w] += src1_ptr[x + y*w] * A;
+                                       darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
                                        sum1_ptr[xx + yy*w] += A;
                                        /* follow position via second plane */
+                                       darrayCheckPr(flow, flow_x2_ptr + xx + yy*w, 1);
                                        _x += flow_x2_ptr[xx + yy*w];
+                                       darrayCheckPr(flow, flow_y2_ptr + xx + yy*w, 1);
                                        _y += flow_y2_ptr[xx + yy*w];
                                        xxx = x + round(_x);
                                        yyy = y + round(_y);
                                        if (xxx<0 || yyy<0 || xxx>=w || yyy>=h)
                                                continue;
                                        /* even plane (weight B) */
+                                       darrayCheckPr(dst, dst2_ptr + xxx + yyy*w, 1);
+                                       darrayCheckPr(src, src1_ptr + x + y*w, 1);
                                        dst2_ptr[xxx + yyy*w] += src1_ptr[x + y*w] * B;
+                                       darrayCheckPr(sum, sum2_ptr + xxx + yyy*w, 1);
                                        sum2_ptr[xxx + yyy*w] += B;
                                        /* odd plane (weigth B) */
+                                       darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+                                       darrayCheckPr(src, src2_ptr + xxx + yyy*w, 1);
                                        dst1_ptr[xx + yy*w] += src2_ptr[xxx + yyy*w] * B;
+                                       darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
                                        sum1_ptr[xx + yy*w] += B;
                                }
                        }
@@ -409,19 +428,24 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
                dst1_ptr = darrayGetPr(dst) + z*w*h;
                sum1_ptr = darrayGetPr(sum) + z*w*h;
                src1_ptr = darrayGetPr(src) + (z/2)*w*h;
-               flow_x1_ptr = darrayGetPr(flow_x) + z*w*h;
-               flow_y1_ptr = darrayGetPr(flow_y) + z*w*h;
+               flow_x1_ptr = darrayGetPr(flow) + z*w*h;
+               flow_y1_ptr = flow_x1_ptr + (kd-1)*w*h;
                for (y = 0; y < h; y++) {
                        for (x = 0; x < w; x++) {
                                /* get position via firt plane */
+                               darrayCheckPr(flow, flow_x1_ptr, 1);
                                _x = *flow_x1_ptr++;
+                               darrayCheckPr(flow, flow_y1_ptr, 1);
                                _y = *flow_y1_ptr++;
                                xx = x + round(_x);
                                yy = y + round(_y);
                                if (xx<0 || yy<0 || xx>=w || yy>=h)
                                        continue;
                                /* border (odd) plane (weight 1.0)*/
+                               darrayCheckPr(dst, dst1_ptr + xx + yy*w, 1);
+                               darrayCheckPr(src, src1_ptr + x + y*w, 1);
                                dst1_ptr[xx + yy*w] += src1_ptr[x + y*w];
+                               darrayCheckPr(sum, sum1_ptr + xx + yy*w, 1);
                                sum1_ptr[xx + yy*w] += 1.0;
                        }
                }
@@ -430,6 +454,8 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
 
        /* copy extra plane, if destiantion size is odd */
        if (z < kd) {
+               darrayCheckPr(dst, darrayGetPr(dst) + z*w*h, w*h);
+               darrayCheckPr(src, darrayGetPr(src) + (z/2)*w*h, w*h);
                memcpy(darrayGetPr(dst) + z*w*h, darrayGetPr(src) + (z/2)*w*h, w*h*sizeof(double));
                /* fill sum array with ones */
                pixles = w*h;
@@ -444,6 +470,8 @@ static int expand_z(darray_t *src, const darray_t *dst, const darray_t *flow_x,
        sum1_ptr = darrayGetPr(sum);
        for (i = 0; i < pixles; i++) {
                /* only if there is a sum */
+               darrayCheckPr(sum, sum1_ptr, 1);
+               darrayCheckPr(dst, dst1_ptr, 1);
                if (*sum1_ptr > 0.1)
                        *dst1_ptr++ /= *sum1_ptr;
                else
@@ -469,91 +497,81 @@ error:
  *
  * If the number of flow planes is odd, the extra plane is removed.
  */
-static int shrink_flow_z(darray_t *src_x, darray_t *src_y, darray_t *dst_x, darray_t *dst_y)
+int shrink_flow_z(darray_t *src, darray_t *dst)
 {
        int w, h, ks, kd, x, y, z, xx, yy;
        double _x, _y;
        double *src_x1_ptr, *src_x2_ptr, *src_y1_ptr, *src_y2_ptr, *dst_x_ptr, *dst_y_ptr;
-       int dst_x_move = 0, dst_y_move = 0;
+       int dst_move = 0;
 
-       ASSERT(darrayGetNumberOfDimensions(src_x) == 3, "src_x array must have 3 dimensions");
-       ASSERT(darrayGetNumberOfDimensions(src_y) == 3, "src_y array must have 3 dimensions");
-       w = darrayGetDimensions(src_x)[0];
-       h = darrayGetDimensions(src_x)[1];
-       ks = darrayGetDimensions(src_x)[2];
+       ASSERT(darrayGetNumberOfDimensions(src) == 4, "src array must have 4 dimensions");
+       w = darrayGetDimensions(src)[0];
+       h = darrayGetDimensions(src)[1];
+       ks = darrayGetDimensions(src)[2];
+       ASSERT(2 == darrayGetDimensions(src)[3], "number of cubes must be 2");
        kd = ks>>1;
-       ASSERT(w == darrayGetDimensions(src_y)[0], "x-size must be equal");
-       ASSERT(h == darrayGetDimensions(src_y)[1], "y-size must be equal");
-       ASSERT(ks == darrayGetDimensions(src_y)[2], "z-size must be equal");
 
-       if (dst_x == NULL) {
-               int dimensions[3] = {w, h, kd};
+       if (dst == NULL) {
+               int dimensions[4] = {w, h, kd, 2};
 
-               dst_x = darrayCreate(3, dimensions);
-               if (!dst_x)
-                       goto error;
-               dst_x_move = 1;
-       }
-       if (dst_y == NULL) {
-               int dimensions[3] = {w, h, kd};
-
-               dst_y = darrayCreate(3, dimensions);
-               if (!dst_y)
+               dst = darrayCreate(4, dimensions);
+               if (!dst)
                        goto error;
-               dst_y_move = 1;
+               dst_move = 1;
        }
 
-       ASSERT(darrayGetNumberOfDimensions(dst_x) == 3, "flow_x array must have 3 dimensions");
-       ASSERT(darrayGetNumberOfDimensions(dst_y) == 3, "flow_y array must have 3 dimensions");
-       ASSERT(w == darrayGetDimensions(dst_x)[0], "x-size must be equal");
-       ASSERT(h == darrayGetDimensions(dst_x)[1], "y-size must be equal");
-       ASSERT(kd == darrayGetDimensions(dst_x)[2], "z-size of dest must be half of the planes + optional odd plane");
-       ASSERT(w == darrayGetDimensions(dst_y)[0], "x-size must be equal");
-       ASSERT(h == darrayGetDimensions(dst_y)[1], "y-size must be equal");
-       ASSERT(kd == darrayGetDimensions(dst_y)[2], "z-size of dest must be half of the planes + optional odd plane");
+       ASSERT(darrayGetNumberOfDimensions(dst) == 4, "dst array must have 4 dimensions");
+       ASSERT(w == darrayGetDimensions(dst)[0], "x-size must be equal");
+       ASSERT(h == darrayGetDimensions(dst)[1], "y-size must be equal");
+       ASSERT(kd == darrayGetDimensions(dst)[2], "z-size of dest must be half of the planes + optional odd plane");
+       ASSERT(2 == darrayGetDimensions(dst)[3], "number of cubes must be 2");
 
        for (z = 0; z < kd; z++) {
-               src_x1_ptr = darrayGetPr(src_x) + (z*2)*w*h;
-               src_y1_ptr = darrayGetPr(src_y) + (z*2)*w*h;
-               dst_x_ptr = darrayGetPr(dst_x) + z*w*h;
-               dst_y_ptr = darrayGetPr(dst_y) + z*w*h;
-               src_x2_ptr = darrayGetPr(src_x) + (z*2+1)*w*h;
-               src_y2_ptr = darrayGetPr(src_y) + (z*2+1)*w*h;
+               src_x1_ptr = darrayGetPr(src) + (z*2)*w*h;
+               src_y1_ptr = src_x1_ptr + w*h*ks;
+               src_x2_ptr = src_x1_ptr + w*h;
+               src_y2_ptr = src_y1_ptr + w*h;
+               dst_x_ptr = darrayGetPr(dst) + z*w*h;
+               dst_y_ptr = dst_x_ptr + w*h*kd;
                for (y = 0; y < h; y++) {
                        for (x = 0; x < w; x++) {
                                /* get position via first plane */
+                               darrayCheckPr(src, src_x1_ptr, 1);
                                _x = *src_x1_ptr++;
+                               darrayCheckPr(src, src_y1_ptr, 1);
                                _y = *src_y1_ptr++;
                                xx = x + round(_x);
                                yy = y + round(_y);
                                if (xx<0 || yy<0 || xx>=w || yy>=h) {
                                        /* use an outside value when outside in the first plane */
+                                       darrayCheckPr(dst, dst_x_ptr, 1);
                                        *dst_x_ptr++ = _x;
+                                       darrayCheckPr(dst, dst_y_ptr, 1);
                                        *dst_y_ptr++ = _y;
                                        continue;
                                }
                                /* follow position via second plane */
+                               darrayCheckPr(src, src_x2_ptr + xx + yy*w, 1);
                                _x += src_x2_ptr[xx + yy*w];
+                               darrayCheckPr(src, src_y2_ptr + xx + yy*w, 1);
                                _y += src_y2_ptr[xx + yy*w];
                                /* store displacement result */
+                               darrayCheckPr(dst, dst_x_ptr, 1);
                                *dst_x_ptr++ = _x;
+                               darrayCheckPr(dst, dst_y_ptr, 1);
                                *dst_y_ptr++ = _y;
                        }
                }
        }
 
-       if (dst_x_move)
-               darrayMove(dst_x, src_x);
-       if (dst_y_move)
-               darrayMove(dst_y, src_y);
+       if (dst_move)
+               darrayMove(dst, src);
 
        return 0;
 
 error:
-       if (dst_x_move)
-               darrayDestroy(dst_x);
-       if (dst_y_move)
-               darrayDestroy(dst_y);
+       if (dst_move)
+               darrayDestroy(dst);
 
        return -1;
 }
@@ -567,40 +585,52 @@ error:
  * If a dimension of a source matrix is odd, the end of the dimension is just
  * copied and not scaled.
  *
- * The points in a plane are shrinked from 2x2x2 to 1x1x1 by calculating the mean
+ * The points in a plane are shrinked from 2x2 to 1x1 by calculating the mean
  * of all eight values.
+ *
+ * If the array has more than two dimensions, all planes in all dimensions
+ * are processed.
  */
-static int shrink_xy(darray_t *src, darray_t *dst)
+int shrink_xy(darray_t *src, darray_t *dst)
 {
-       int x, xx, y, z, ws, wd, w_, hs, hd, h_, k;
+       int x, xx, y, z, ws, wd, w_, hs, hd, h_, k, dims, i;
        double *s1, *s2, *d;
        int dst_move = 0;
-       ASSERT(darrayGetNumberOfDimensions(src) == 3, "Source array must have 3 dimensions");
+
+       dims = darrayGetNumberOfDimensions(src);
+       ASSERT(dims >= 2, "Source array must have at least 2 dimensions");
        ws = darrayGetDimensions(src)[0];
        w_ = ws>>1;
        hs = darrayGetDimensions(src)[1];
        h_ = hs>>1;
-       k = darrayGetDimensions(src)[2];
+       k = 1;
+       for (i = 2; i < dims; i++)
+               k *= darrayGetDimensions(src)[i];
 
        /* create dst, if not given */
        if (!dst) {
-               int dimensions[3] = {(ws+1)>>1, (hs+1)>>1, k};
-
-               dst = darrayCreate(3, dimensions);
+               int dimensions[dims];
+               
+               dimensions[0] = (ws+1)>>1;
+               dimensions[1] = (hs+1)>>1;
+               for (i = 2; i < dims; i++)
+                       dimensions[i] = darrayGetDimensions(src)[i];
+
+               dst = darrayCreate(dims, dimensions);
                if (!dst)
                        goto error;
                dst_move = 1;
        }
 
-       ASSERT(darrayGetNumberOfDimensions(dst) == 3, "Destination array must have 3 dimensions");
+       ASSERT(darrayGetNumberOfDimensions(dst) == dims, "Destination array must have equal dimensions than source array");
        wd = darrayGetDimensions(dst)[0];
        hd = darrayGetDimensions(dst)[1];
        ASSERT((ws&1) || ws == wd*2, "x-size of source must be twice the size of destination (this is even source size)");
        ASSERT(!(ws&1) || ws-1 == (wd-1)*2, "x-size of source-1 must be twice the size of destination-1 (this odd source size)");
        ASSERT((hs&1) || hs == hd*2, "y-size of source must be twice the size of destination (this is even source size)");
        ASSERT(!(hs&1) || hs-1 == (hd-1)*2, "y-size of source-1 must be twice the size of destination-1 (this odd source size)");
-       ASSERT(k == darrayGetDimensions(dst)[2], "z-size must be equal");
+       for (i = 2; i < dims; i++)
+               ASSERT(darrayGetDimensions(src)[i] == darrayGetDimensions(dst)[i], "dimension 3 and greater must be equal");
 
        for (z = 0; z < k; z++) {
                for (y = 0; y < h_; y++) {
@@ -608,21 +638,35 @@ static int shrink_xy(darray_t *src, darray_t *dst)
                        s1 = darrayGetPr(src) + y*2*ws + z*ws*hs;
                        s2 = s1 + ws;
                        d = darrayGetPr(dst) + y*wd + z*wd*hd;
-                       for (x = 0, xx = 0; x < w_; x++, xx+=2)
+                       for (x = 0, xx = 0; x < w_; x++, xx+=2) {
+                               darrayCheckPr(dst, d + x, 1);
+                               darrayCheckPr(src, s1 + xx, 2);
+                               darrayCheckPr(src, s2 + xx, 2);
                                d[x] = (s1[xx] + s1[xx+1] + s2[xx] + s2[xx+1]) * 0.25;
+                       }
                        /* mean of 1x2x1 */
-                       if (x < wd)
+                       if (x < wd) {
+                               darrayCheckPr(dst, d + x, 1);
+                               darrayCheckPr(src, s1 + xx, 1);
+                               darrayCheckPr(src, s2 + xx, 1);
                                d[x] = (s1[xx] + s2[xx]) * 0.5;
+                       }
                }
                if (y < hd) {
                        /* mean of 2x1x1 */
                        s1 = darrayGetPr(src) + y*2*ws + z*ws*hs;
                        d = darrayGetPr(dst) + y*wd + z*wd*hd;
-                       for (x = 0, xx = 0; x < w_; x++, xx+=2)
+                       for (x = 0, xx = 0; x < w_; x++, xx+=2) {
+                               darrayCheckPr(dst, d + x, 1);
+                               darrayCheckPr(src, s1 + xx, 2);
                                d[x] = (s1[xx] + s1[xx+1]) * 0.5;
+                       }
                        /* mean of 1x1x1 */
-                       if (x < wd)
+                       if (x < wd) {
+                               darrayCheckPr(dst, d + x, 1);
+                               darrayCheckPr(src, s1 + xx, 1);
                                d[x] = s1[xx];
+                       }
                }
        }
 
@@ -692,14 +736,21 @@ static void expand_xy(darray_t *src, const darray_t *dst)
                        d = darrayGetPr(dst) + (y+h_)*wd + z*wd*hd;
                        *d++ = *s; /* x=0 plane */
                        for (x = 0; x < w_-1; x++) {
+                               darrayCheckPr(dst, d, 1);
+                               darrayCheckPr(src, s, 2);
                                *d++ = A*s[0] + B*s[1];
                                *d++ = A*s[1] + B*s[0];
                                s++;
                        }
+                       darrayCheckPr(dst, d, 1);
+                       darrayCheckPr(src, s, 1);
                        *d++ = *s++; /* x=w-1 plane */
                        /* extra plane (odd fine grid) */
-                       if ((wd&1))
+                       if ((wd&1)) {
+                               darrayCheckPr(dst, d, 1);
+                               darrayCheckPr(src, s, 1);
                                *d = *s;
+                       }
                }
        }
 
@@ -709,14 +760,21 @@ static void expand_xy(darray_t *src, const darray_t *dst)
                /* y=h-1 plane */
                s = darrayGetPr(dst) + h_*wd + z*wd*hd;
                d = darrayGetPr(dst) + z*wd*hd;
-               for (x = 0; x < wd; x++)
+               for (x = 0; x < wd; x++) {
+                       darrayCheckPr(dst, d, 1);
+                       darrayCheckPr(dst, s, 1);
                        *d++ = *s++; /* y=0 plane */
+               }
                for (y = 0; y < h_-1; y++) {
                        s1 = darrayGetPr(dst) + (y+h_)*wd + z*wd*hd;
                        s2 = s1 + wd;
                        d1 = darrayGetPr(dst) + (y*2+1)*wd + z*wd*hd;
                        d2 = d1 + wd;
                        for (x = 0; x < wd; x++) {
+                               darrayCheckPr(dst, d1, 1);
+                               darrayCheckPr(dst, d2, 1);
+                               darrayCheckPr(dst, s1, 1);
+                               darrayCheckPr(dst, s2, 1);
                                *d1++ = A * *s1 + B * *s2;
                                *d2++ = A * *s2 + B * *s1;
                                s1++;
@@ -731,18 +789,18 @@ static void expand_xy(darray_t *src, const darray_t *dst)
  *
  * If an element is greater theshold, the element is set to gt, otherwise le.
  */
-static void quantize(const darray_t *array, double threshold, double gt, double le)
+void quantize(const darray_t *array, double threshold, double gt, double le)
 {
-       int n = 1, i;
+       int n = 1, dims, i;
        double *ptr;
 
-       ASSERT(darrayGetNumberOfDimensions(array) == 3, "array must have 3 dimensions");
-       n *= darrayGetDimensions(array)[0];
-       n *= darrayGetDimensions(array)[1];
-       n *= darrayGetDimensions(array)[2];
+       dims = darrayGetNumberOfDimensions(array);
+       for (i = 0; i < dims; i++)
+               n *= darrayGetDimensions(array)[i];
        ptr = darrayGetPr(array);
 
        for (i = 0; i < n; i++, ptr++) {
+               darrayCheckPr(array, ptr, 1);
                if (*ptr > threshold)
                        *ptr = gt;
                else
@@ -752,26 +810,27 @@ static void quantize(const darray_t *array, double threshold, double gt, double
 
 /* multiply the array content with a value
  */
-static void multiply(const darray_t *array, double factor)
+void multiply(const darray_t *array, double factor)
 {
-       int n = 1, i;
+       int n = 1, dims, i;
        double *ptr;
 
-       ASSERT(darrayGetNumberOfDimensions(array) == 3, "array must have 3 dimensions");
-       n *= darrayGetDimensions(array)[0];
-       n *= darrayGetDimensions(array)[1];
-       n *= darrayGetDimensions(array)[2];
+       dims = darrayGetNumberOfDimensions(array);
+       for (i = 0; i < dims; i++)
+               n *= darrayGetDimensions(array)[i];
        ptr = darrayGetPr(array);
 
-       for (i = 0; i < n; i++, ptr++)
+       for (i = 0; i < n; i++, ptr++) {
+               darrayCheckPr(array, ptr, 1);
                *ptr *= factor;
+       }
 }
 
 /* create an array of neighbors and use flow vectors
  * flow_x and flow_x are the vectors that point to neighbors of the next z plane
  * flow_x and flow_x are the vectors that point to neighbors of the previous z plane
  */
-static unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow, const darray_t *flow_i)
 {
        int x, y, z, xx, yy, zz, i, j, n;
        double _x, _y;
@@ -782,11 +841,11 @@ static unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow_x,
        int nb_num;
        int *i_index;
 
-       ASSERT_FLOW(w, h, k, flow_x, flow_y, flow_ix, flow_iy);
-       flow_x_ptr = darrayGetPr(flow_x);
-       flow_y_ptr = darrayGetPr(flow_y);
-       flow_ix_ptr = darrayGetPr(flow_ix);
-       flow_iy_ptr = darrayGetPr(flow_iy);
+       ASSERT_FLOW(w, h, k, flow, flow_i);
+       flow_x_ptr = darrayGetPr(flow);
+       flow_y_ptr = flow_x_ptr + w*h*(k-1);
+       flow_ix_ptr = darrayGetPr(flow_i);
+       flow_iy_ptr = flow_ix_ptr + w*h*(k-1);
 
        nb_num = (k==1) ? NEIGHBORS_2D : NEIGHBORS_3D;
        nb_size = (k==1) ? sizeof(struct neighbors_2d) : sizeof(struct neighbors_3d);
@@ -816,12 +875,16 @@ static unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow_x,
                                        if (xx >= w || yy >= h || zz >= k)
                                                continue;
                                        if (zz > z) {
+                                               darrayCheckPr(flow, flow_x_ptr + xx + yy*w + (zz-1)*w*h, 1);
                                                _x = flow_x_ptr[xx + yy*w + (zz-1)*w*h];
+                                               darrayCheckPr(flow, flow_y_ptr + xx + yy*w + (zz-1)*w*h, 1);
                                                _y = flow_y_ptr[xx + yy*w + (zz-1)*w*h];
                                                xx = xx + round(_x);
                                                yy = yy + round(_y);
                                        } else if (zz < z) {
+                                               darrayCheckPr(flow_i, flow_ix_ptr + xx + yy*w + zz*w*h, 1);
                                                _x = flow_ix_ptr[xx + yy*w + zz*w*h];
+                                               darrayCheckPr(flow_i, flow_iy_ptr + xx + yy*w + zz*w*h, 1);
                                                _y = flow_iy_ptr[xx + yy*w + zz*w*h];
                                                xx = xx + round(_x);
                                                yy = yy + round(_y);
@@ -853,7 +916,7 @@ static unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow_x,
 }
 
 /* weighten neighbors */
-static void weighten(const darray_t *luminance, unsigned char *neighbors)
+void weighten(const darray_t *luminance, unsigned char *neighbors)
 {
        int x, y, z, i, w, h, k, num;
        double val, nval, mean, variance, sigma, sum, weight, min;
@@ -894,6 +957,7 @@ static void weighten(const darray_t *luminance, unsigned char *neighbors)
                                variance = 0;
                                for (i = 0; i < num; i++) {
                                        /* value of neighbor pixles, including current pixle */
+                                       darrayCheckPr(luminance, lum_ptr + i_index[i], 1);
                                        nval = lum_ptr[i_index[i]];
                                        /* squared diffs between current pixle and each neighbor */
                                        sqdiff[i] = (nval-val) * (nval-val);
@@ -981,25 +1045,36 @@ static void smooth(int level, darray_t *value, const darray_t *mark, const darra
                                        w_index = ((struct neighbors_3d *)neighbors)->weight;
                                }
                                /* if pixle is marked */
+                               darrayCheckPr(mark, mark_line_ptr + x, 1);
                                if (mark_line_ptr[x] < 0.5) {
                                        /* always use initial color on marked plane, otherwise do not smooth */
-                                       if (level == 0)
+                                       if (level == 0) {
+                                               darrayCheckPr(value, value_line_ptr + x, 1);
+                                               darrayCheckPr(init, init_line_ptr + x, 1);
                                                value_line_ptr[x] = init_line_ptr[x];
+                                       }
                                        continue;
                                }
                                /* Calculate value */
                                val = 0;
-                               for (i = 1; i < num; i++) /* sum value without point itself */
+                               for (i = 1; i < num; i++) { /* sum value without point itself */
+                                       darrayCheckPr(value, value_ptr + i_index[i], 1);
                                        val += w_index[i] * value_ptr[i_index[i]];
-                               if (level == 0)
+                               }
+                               if (level == 0) {
+                                       darrayCheckPr(value, value_line_ptr + x, 1);
+                                       darrayCheckPr(init, init_line_ptr + x, 1);
                                        value_line_ptr[x] = (init_line_ptr[x] - val) / w_index[0];
-                               else
+                               } else {
+                                       darrayCheckPr(value, value_line_ptr + x, 1);
                                        value_line_ptr[x] = (-val) / w_index[0];
+                               }
                        }
                }
        }
 }
 
+#ifdef RESIDUAL
 /* Calculate residual
  */
 static double residual(int level, darray_t *value, const darray_t *mark, const darray_t *init, unsigned char *neighbors)
@@ -1038,15 +1113,19 @@ static double residual(int level, darray_t *value, const darray_t *mark, const d
                                        w_index = ((struct neighbors_3d *)neighbors)->weight;
                                }
                                /* if pixle is marked, no residual */
+                               darrayCheckPr(mark, mark_line_ptr + x, 1);
                                if (mark_line_ptr[x] < 0.5)
                                        continue;
                                /* Calculate value */
                                val = 0;
-                               for (i = 0; i < num; i++) /* sum value including point itself */
+                               for (i = 0; i < num; i++) { /* sum value including point itself */
+                                       darrayCheckPr(value, value_ptr + i_index[i], 1);
                                        val += w_index[i] * value_ptr[i_index[i]];
-                               if (level == 0)
+                               }
+                               if (level == 0) {
+                                       darrayCheckPr(init, init_line_ptr + x, 1);
                                        val = init_line_ptr[x] - val;
-                               else
+                               else
                                        val = -val;
                                residual += val * val;
                        }
@@ -1055,6 +1134,7 @@ static double residual(int level, darray_t *value, const darray_t *mark, const d
 
        return residual;
 }
+#endif
 
 /* Restrict
  *
@@ -1064,9 +1144,10 @@ static double residual(int level, darray_t *value, const darray_t *mark, const d
  * different depth (z dimension), a temporary array (temp) is used to
  * perform the intermediate step of scaling x,y,z.
  *
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
  */
-static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flow_x, const darray_t *flow_y)
+static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flow)
 {
        int ks, wd, hd, kd;
        darray_t *temp = NULL;
@@ -1079,7 +1160,7 @@ static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flo
        kd = darrayGetDimensions(nvalue)[2];
        ks = darrayGetDimensions(value)[2];
        if (ks != kd)
-               ASSERT(flow_x && flow_y, "when scaling z, flow vectors must be given");
+               ASSERT(flow, "when scaling z, flow vectors must be given");
 
        if (ks != kd) {
                int dimensions[3] = {wd, hd, ks};
@@ -1092,7 +1173,7 @@ static int restrict(const darray_t *value, darray_t *nvalue, const darray_t *flo
                rc = shrink_xy((darray_t *)value, temp);
                if (rc < 0)
                        goto error;
-               rc = shrink_z(temp, nvalue, flow_x, flow_y);
+               rc = shrink_z(temp, nvalue, flow);
                if (rc < 0)
                        goto error;
        } else {
@@ -1117,7 +1198,8 @@ error:
  * different depth (z dimension), a temporary array (temp3) is used to
  * perform the intermediate step of scaling x,y,z.
  *
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
  *
  * First the value is shrinked to the size of the value of the next level.
  * The error is calculated by substracting the values of the next level
@@ -1125,7 +1207,7 @@ error:
  * substracted from the value. So the changes in the next level's value are
  * prolonged to the value.
  */
-static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow_x, const darray_t *flow_y)
+static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow)
 {
        int ws, hs, ks, wd, hd, kd;
        int i, pixles;
@@ -1142,7 +1224,7 @@ static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow
        hd = darrayGetDimensions(value)[1];
        kd = darrayGetDimensions(value)[2];
        if (ks != kd)
-               ASSERT(flow_x && flow_y, "when scaling z, flow vectors must be given");
+               ASSERT(flow, "when scaling z, flow vectors must be given");
 
        /* generate temporary arrays */
        temp1 = darrayCreate(3, darrayGetDimensions(value));
@@ -1165,7 +1247,7 @@ static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow
                rc = shrink_xy(value, temp3);
                if (rc < 0)
                        goto error;
-               rc = shrink_z(temp3, temp2, flow_x, flow_y);
+               rc = shrink_z(temp3, temp2, flow);
                if (rc < 0)
                        goto error;
        } else {
@@ -1178,12 +1260,15 @@ static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow
        temp2_ptr = darrayGetPr(temp2);
        nvalue_ptr = darrayGetPr(nvalue);
        pixles = ws*hs*ks;
-       for (i = 0; i < pixles; i++)
+       for (i = 0; i < pixles; i++) {
+               darrayCheckPr(temp2, temp2_ptr, 1);
+               darrayCheckPr(nvalue, nvalue_ptr, 1);
                *temp2_ptr++ -= *nvalue_ptr++;
+       }
 
        /* Expand (temp2 -> temp1) */
        if (ks != kd) {
-               rc = expand_z(temp2, temp3, flow_x, flow_y);
+               rc = expand_z(temp2, temp3, flow);
                if (rc < 0)
                        goto error;
                expand_xy(temp3, temp1);
@@ -1194,8 +1279,11 @@ static int prolong(darray_t *value, const darray_t *nvalue, const darray_t *flow
        temp1_ptr = darrayGetPr(temp1);
        value_ptr = darrayGetPr(value);
        pixles = wd*hd*kd;
-       for (i = 0; i < pixles; i++)
+       for (i = 0; i < pixles; i++) {
+               darrayCheckPr(temp1, temp1_ptr, 1);
+               darrayCheckPr(value, value_ptr, 1);
                *value_ptr++ -= *temp1_ptr++;
+       }
 
        rc = 0;
 
@@ -1226,6 +1314,8 @@ void paste(const darray_t *src, darray_t *dst, int xoff, int yoff, int zoff)
        for (z = 0; z < ks; z++) {
                for (y = 0; y < hs; y++) {
                        for (x = 0; x < ws; x++) {
+                               darrayCheckPr(dst, dst_ptr + (x+xoff) + (y+yoff)*wd + (z+zoff)*wd*hd, 1);
+                               darrayCheckPr(src, src_ptr + x + y*ws + z*ws*hs, 1);
                                dst_ptr[(x+xoff) + (y+yoff)*wd + (z+zoff)*wd*hd] = src_ptr[x + y*ws + z*ws*hs];
                        }
                }
@@ -1271,7 +1361,7 @@ inline void draw_line(double *ptr, int w, int h, int x0, int y0, int x1, int y1,
  * The flow vectors are drawn to the upper left corner of the dst, which
  * may be larger in size than the flow arrays.
  */
-void draw_flow(darray_t *dst, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy)
+void draw_flow(darray_t *dst, const darray_t *flow, const darray_t *flow_i)
 {
        int wd, hd, ws, hs, ks, x, y, z;
        double *flow_x_ptr, *flow_y_ptr, *flow_ix_ptr, *flow_iy_ptr;
@@ -1279,13 +1369,13 @@ void draw_flow(darray_t *dst, const darray_t *flow_x, const darray_t *flow_y, co
 
        wd = darrayGetDimensions(dst)[0];
        hd = darrayGetDimensions(dst)[1];
-       ws = darrayGetDimensions(flow_x)[0];
-       hs = darrayGetDimensions(flow_x)[1];
-       ks = darrayGetDimensions(flow_x)[2];
-       flow_x_ptr = darrayGetPr(flow_x);
-       flow_y_ptr = darrayGetPr(flow_y);
-       flow_ix_ptr = darrayGetPr(flow_ix);
-       flow_iy_ptr = darrayGetPr(flow_iy);
+       ws = darrayGetDimensions(flow)[0];
+       hs = darrayGetDimensions(flow)[1];
+       ks = darrayGetDimensions(flow)[2];
+       flow_x_ptr = darrayGetPr(flow);
+       flow_y_ptr = flow_x_ptr + ws*hs*ks;
+       flow_ix_ptr = darrayGetPr(flow_i);
+       flow_iy_ptr = flow_ix_ptr + ws*hs*ks;
 
        for (z = 0; z <= ks; z++) {
                for (y = step/2; y < hs-step/2+1; y+=step) {
@@ -1303,25 +1393,33 @@ void draw_flow(darray_t *dst, const darray_t *flow_x, const darray_t *flow_y, co
 
 /* Solve multigrid
  *
- * Note that the flow vectors must be scaled xy already!
+ * Note that the flow vectors must be scaled xy already (if x and y is shrinked
+ * before z)!
+ *
  */
-static int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows_xy_x, darray_t **flows_xy_y, int scalexyz)
+int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows, int scalexyz)
 {
+#ifdef RESIDUAL
        double curr_residual, prev_residual = 10000000000000.0;
+#endif
        int cycle, level, i;
        int rc = -1;
 
        for (cycle = 0; cycle < ncycle; cycle++) {
                /* go down */
+#ifdef RESIDUAL
                curr_residual = 0;
+#endif
                for (level = 0; level < nlevel-1; level++) {
                        for (i = 0; i < iters; i++)
                                smooth(level, values[level], marks[level], init, nb_list[level]);
+#ifdef RESIDUAL
                        curr_residual += residual(level, values[level], marks[level], init, nb_list[level]);
+#endif
                        if (level < scalexyz)
-                               rc = restrict(values[level], values[level+1], flows_xy_x[level], flows_xy_y[level]);
+                               rc = restrict(values[level], values[level+1], flows[level]);
                        else
-                               rc = restrict(values[level], values[level+1], NULL, NULL);
+                               rc = restrict(values[level], values[level+1], NULL);
                        if (rc < 0)
                                goto error;
                }
@@ -1333,9 +1431,9 @@ static int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray
                /* go up */
                for (level = nlevel-2; level >= 0; level--) {
                        if (level < scalexyz)
-                               rc = prolong(values[level], values[level+1], flows_xy_x[level], flows_xy_y[level]);
+                               rc = prolong(values[level], values[level+1], flows[level]);
                        else
-                               rc = prolong(values[level], values[level+1], NULL, NULL);
+                               rc = prolong(values[level], values[level+1], NULL);
                        if (rc < 0)
                                goto error;
                        for (i = 0; i < iters; i++)
@@ -1343,9 +1441,11 @@ static int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray
                }
 
                /* end if residual raises again */
-//             if (prev_residual < curr_residual)
-//                     break;
+#ifdef RESIDUAL
+               if (prev_residual < curr_residual)
+                       break;
                curr_residual = prev_residual;
+#endif
        }
 
        return 0;
@@ -1354,457 +1454,3 @@ error:
        return rc;
 }
 
-/* colorize image array using image_mark
- *
- * image:
- *     first dimansion = width
- *     second dimension = height
- *     thrid dimension = components (first component is luminance, other components are features, like u, v)
- *     fourth dimension = time
- *
- * image_mark:
- *     first dimansion = width
- *     second dimension = height
- *     fourth dimension = time
- */
-darray_t *colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow_x, const darray_t *flow_y, const darray_t *flow_ix, const darray_t *flow_iy, int inner_iter, int outer_iter, int scalexyz)
-{
-       unsigned char **nb_list = NULL;
-       darray_t **values = NULL, **marks = NULL, *init = NULL, **flows_x = NULL, **flows_y = NULL, **flows_ix = NULL, **flows_iy = NULL, **flows_xy_x = NULL, **flows_xy_y = NULL;
-       darray_t *result = NULL, *luminance = NULL;
-       int dimensions[4];
-       int w, h, k, ww, hh, kk, last_kk, comp, d, max_depth, temp, z, c;
-       int rc = -1;
-
-       ASSERT(darrayGetNumberOfDimensions(image) == 4, "image array must have 4 dimensions");
-       w = darrayGetDimensions(image)[0];
-       h = darrayGetDimensions(image)[1];
-       comp = darrayGetDimensions(image)[2];
-       k = darrayGetDimensions(image)[3];
-       ASSERT(comp > 0, "image array requires at least one component (Y-component)");
-
-       ASSERT(darrayGetNumberOfDimensions(image_mark) == 3, "image_mark array must have 3 dimensions");
-       ASSERT(w == darrayGetDimensions(image_mark)[0], "image_mark and image must both have same x-size");
-       ASSERT(h == darrayGetDimensions(image_mark)[1], "image_mark and image must both have same y-size");
-       ASSERT(k == darrayGetDimensions(image_mark)[2], "image_mark and image must both have same z-size");
-
-       /* Calculate maximum depth
-        * Scale down dimensions until at least 2x2 remains
-        */
-       temp = min(w, h);
-       for (max_depth = 0; temp > 2; max_depth++)
-               temp=(temp+1)>>1;
-       ASSERT(max_depth >= 2, "image too small");
-
-       nb_list = calloc(max_depth, sizeof(struct neighbors *));
-       if (!nb_list)
-               goto error;
-       marks = calloc(max_depth, sizeof(darray_t *));
-       if (!marks)
-               goto error;
-       ASSERT(marks, "no memory");
-       values = calloc(max_depth, sizeof(darray_t *));
-       if (!values)
-               goto error;
-       ASSERT(values, "no memory");
-       flows_x = calloc(max_depth, sizeof(darray_t *));
-       if (!flows_x)
-               goto error;
-       ASSERT(flows_x, "no memory");
-       flows_y = calloc(max_depth, sizeof(darray_t *));
-       if (!flows_y)
-               goto error;
-       ASSERT(flows_y, "no memory");
-       flows_ix = calloc(max_depth, sizeof(darray_t *));
-       if (!flows_ix)
-               goto error;
-       ASSERT(flows_ix, "no memory");
-       flows_iy = calloc(max_depth, sizeof(darray_t *));
-       if (!flows_iy)
-               goto error;
-       ASSERT(flows_iy, "no memory");
-       if (scalexyz) {
-               /* the flows_xy are scaled in xy in the current level */
-               flows_xy_x = calloc(max_depth, sizeof(darray_t *));
-               if (!flows_xy_x)
-                       goto error;
-               flows_xy_y = calloc(max_depth, sizeof(darray_t *));
-               if (!flows_xy_y)
-                       goto error;
-       }
-       dimensions[0] = w;
-       dimensions[1] = h;
-       dimensions[2] = k;
-       luminance = darrayCreate(3, dimensions);
-       if (!luminance)
-               goto error;
-       /* generate sparse multigrid */
-       last_kk = kk = k;
-       printf("-> prepare sparse grid");
-       for (d = 0, ww=w, hh=h; d < max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
-               /* create temporary arrays
-                * copy fist array to level 0 or shrink previous level
-                */
-               printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
-               dimensions[0] = ww;
-               dimensions[1] = hh;
-               dimensions[2] = last_kk;
-               marks[d] = darrayCreate(3, dimensions);
-               if (!marks[d])
-                       goto error;
-               dimensions[0] = ww;
-               dimensions[1] = hh;
-               dimensions[2] = last_kk-1;
-               flows_x[d] = darrayCreate(3, dimensions);
-               if (!flows_x[d])
-                       goto error;
-               flows_y[d] = darrayCreate(3, dimensions);
-               if (!flows_y[d])
-                       goto error;
-               flows_ix[d] = darrayCreate(3, dimensions);
-               if (!flows_ix[d])
-                       goto error;
-               flows_iy[d] = darrayCreate(3, dimensions);
-               if (!flows_iy[d])
-                       goto error;
-               if (d < scalexyz) {
-                       dimensions[0] = (ww+1)>>1;
-                       dimensions[1] = (hh+1)>>1;
-                       dimensions[2] = kk-1;
-                       flows_xy_x[d] = darrayCreate(3, dimensions);
-                       if (!flows_xy_x[d])
-                               goto error;
-                       flows_xy_y[d] = darrayCreate(3, dimensions);
-                       if (!flows_xy_y[d])
-                               goto error;
-               }
-               if (d == 0) {
-                       for (z = 0; z < kk; z++)
-                               memcpy(darrayGetPr(luminance) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-                       if (flow_x && flow_y && flow_ix && flow_iy) {
-                               memcpy(darrayGetPr(flows_x[d]), darrayGetPr(flow_x), w*h*(kk-1)*sizeof(double));
-                               memcpy(darrayGetPr(flows_y[d]), darrayGetPr(flow_y), w*h*(kk-1)*sizeof(double));
-                               memcpy(darrayGetPr(flows_ix[d]), darrayGetPr(flow_ix), w*h*(kk-1)*sizeof(double));
-                               memcpy(darrayGetPr(flows_iy[d]), darrayGetPr(flow_iy), w*h*(kk-1)*sizeof(double));
-                       }
-                       memcpy(darrayGetPr(marks[d]), darrayGetPr(image_mark), w*h*kk*sizeof(double));
-                       /* invert, so 0 = marked */
-                       quantize(marks[d], 0.5, 0, 1);
-               } else {
-                       rc = shrink_xy(luminance, NULL);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(marks[d-1], marks[d]);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(flows_x[d-1], flows_x[d]);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(flows_y[d-1], flows_y[d]);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(flows_ix[d-1], flows_ix[d]);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(flows_iy[d-1], flows_iy[d]);
-                       if (rc < 0)
-                               goto error;
-                       /* shrink length of flow vector on smaller grid.
-                        * This is not performed when comparing with
-                        * the original algorithm. (it is a bug there)
-                        */
-                       multiply(flows_x[d], 0.5);
-                       multiply(flows_y[d], 0.5);
-                       multiply(flows_ix[d], 0.5);
-                       multiply(flows_iy[d], 0.5);
-                       if (d-1 < scalexyz) { /* process if previous d shall be shrinked in z direction */
-                               /* shrink time */
-                               rc = shrink_z(luminance, NULL, flows_x[d], flows_y[d]);
-                               if (rc < 0)
-                                       goto error;
-                               rc = shrink_z(marks[d], NULL, flows_x[d], flows_y[d]);
-                               if (rc < 0)
-                                       goto error;
-                               rc = shrink_flow_z(flows_x[d], flows_y[d], NULL, NULL);
-                               if (rc < 0)
-                                       goto error;
-                               rc = shrink_flow_z(flows_ix[d], flows_iy[d], NULL, NULL);
-                               if (rc < 0)
-                                       goto error;
-                       }
-                       /* If any of the restricted pixles is 0 (marked),
-                        * the result must be 0 (marked) too.
-                        */
-                       quantize(marks[d], 0.9999, 1, 0);
-               }
-               if (d < scalexyz) {
-                       /* create x and y shrinked flow vectors for restriction and prolonging process */
-                       rc = shrink_xy(flows_x[d], flows_xy_x[d]);
-                       if (rc < 0)
-                               goto error;
-                       rc = shrink_xy(flows_y[d], flows_xy_y[d]);
-                       if (rc < 0)
-                               goto error;
-                       multiply(flows_xy_x[d], 0.5);
-                       multiply(flows_xy_y[d], 0.5);
-               }
-
-               /* generate values array */
-               dimensions[0] = ww;
-               dimensions[1] = hh;
-               dimensions[2] = kk;
-               values[d] = darrayCreate(3, dimensions);
-               if (!values[d])
-                       goto error;
-
-               /* generate array of neighbors */
-#ifndef TEST_NO_RENDER
-               nb_list[d] = gen_neighbor(ww, hh, kk, flows_x[d], flows_y[d], flows_ix[d], flows_iy[d]);
-               if (!nb_list[d])
-                       goto error;
-               /* weighten */
-               weighten(luminance, nb_list[d]);
-#endif
-
-               last_kk = kk;
-               if (d < scalexyz)
-                       kk = (kk+1)>>1;
-       }
-       darrayDestroy(luminance);
-       luminance = NULL;
-       printf("\n");
-
-       dimensions[0] = w;
-       dimensions[1] = h;
-       dimensions[2] = k;
-       init = darrayCreate(3, dimensions);
-       if (!init)
-               goto error;
-
-       /* result array for all components */
-       dimensions[0] = w;
-       dimensions[1] = h;
-       dimensions[2] = comp;
-       dimensions[3] = k;
-       result = darrayCreate(4, dimensions);
-       if (!result)
-               goto error;
-
-#ifdef TEST_PYRAMID
-       dimensions[0] = w;
-       dimensions[1] = h;
-       dimensions[2] = k;
-       luminance = darrayCreate(3, dimensions);
-       for (z = 0; z < k; z++) {
-               memcpy(darrayGetPr(luminance) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-               memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-#ifdef TEST_FLOW
-               draw_flow(values[0], flows_x[0], flows_y[0], flows_ix[0], flows_iy[0]);
-#endif
-#ifdef TEST_MARKED
-               paste(marks[0], values[0], 0, 0, 0);
-#endif
-
-       }
-       for (d = 1; d < max_depth; d++) {
-               rc = shrink_xy(luminance, NULL);
-               if (rc < 0)
-                       goto error;
-               if (d-1 < scalexyz)
-                       shrink_z(luminance, NULL, flows_xy_x[d-1], flows_xy_y[d-1]);
-               paste(luminance, values[0], 0, 0, 0);
-#ifdef TEST_FLOW
-               draw_flow(values[0], flows_x[d], flows_y[d], flows_ix[d], flows_iy[d]);
-#endif
-#ifdef TEST_MARKED
-               paste(marks[d], values[0], 0, 0, 0);
-#endif
-       }
-       darrayDestroy(luminance);
-       luminance = NULL;
-       /* copy test pyramid image as Y-component */
-       for (z = 0; z < k; z++)
-               memcpy(darrayGetPr(result) + z*w*h*comp, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-#else
-       /* copy Y-component */
-       for (z = 0; z < k; z++)
-               memcpy(darrayGetPr(result) + z*w*h*comp, darrayGetPr(image) + z*w*h*comp, w*h*sizeof(double));
-#endif
-
-       printf("-> apply component to grid and solve");
-       for (c = 1; c < comp; c++) {
-               int pixles, j;
-               double *value_ptr, *init_ptr, *mark_ptr;
-
-               printf(" #%d", c); fflush(stdout);
-               /* apply component
-                * - copy component into value array (level 0), otherwise 0
-                * - use maked colors for init array, otherwise use 0
-                */
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(image) + c*w*h + z*w*h*comp, w*h*sizeof(double));
-               for (d = 1, ww=(w+1)>>1, hh=(h+1)>>1, kk=k; d < max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
-                       if (d-1 < scalexyz) /* scale kk, if last array was scaled */
-                               kk=(kk+1)>>1;   
-                       value_ptr = darrayGetPr(values[d]);
-                       pixles = ww*hh*kk;
-                       for (j = 0; j < pixles; j++)
-                               *value_ptr++ = 0.0;
-               }
-               value_ptr = darrayGetPr(values[0]);
-               init_ptr = darrayGetPr(init);
-               mark_ptr = darrayGetPr(image_mark);
-               pixles = w*h*k;
-               for (j = 0; j < pixles; j++) {
-                       if (*mark_ptr++ > 0.5)
-                               *init_ptr++ = *value_ptr;
-                       else
-                               *init_ptr++ = 0.0;
-                       value_ptr++;
-               }
-
-#ifndef TEST_NO_RENDER
-               /* colorize component */
-               int i;
-               for (i = 0; i < outer_iter; i++) {
-                       rc = solve_mg(max_depth, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
-                       if (rc < 0)
-                               goto error;
-                       rc = solve_mg((max_depth+1) / 2, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
-                       if (rc < 0)
-                               goto error;
-                       rc = solve_mg(2, inner_iter, 2, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
-                       if (rc < 0)
-                               goto error;
-                       rc = solve_mg(1, inner_iter, 4, values, marks, init, nb_list, flows_xy_x, flows_xy_y, scalexyz);
-                       if (rc < 0)
-                               goto error;
-               }
-#endif
-
-#ifdef TEST_PYRAMID
-               for (d = 1; d < max_depth; d++)
-                       paste(values[d], values[0], 0, 0, 0);
-#endif
-
-               /* apply colorized component */
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-       }
-       printf("\n");
-
-#ifdef TEST_SHRINK_XY
-       for (c = 0; c < comp; c++) {
-               darray_t *temp;
-               int dimensions[3] = {w, h, k};
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
-               dimensions[0] = (w+1)>>1;
-               dimensions[1] = (h+1)>>1;
-               temp = darrayCreate(3, dimensions);
-               ASSERT(temp, "no memory");
-               rc = shrink_xy(values[0], temp);
-               ASSERT(rc >= 0, "no memory");
-               paste(temp, values[0], 0, 0, 0);
-               darrayDestroy(temp);
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-       }
-#endif
-
-#ifdef TEST_SHRINK_Z
-       for (c = 0; c < comp; c++) {
-               darray_t *temp;
-               int dimensions[3] = {w, h, k};
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
-               dimensions[2] = (k+1)>>1;
-               temp = darrayCreate(3, dimensions);
-               ASSERT(temp, "no memory");
-               rc = shrink_z(values[0], temp, flows_x[0], flows_y[0]);
-               ASSERT(rc >= 0, "no memory");
-               paste(temp, values[0], 0, 0, 0);
-               darrayDestroy(temp);
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-       }
-#endif
-
-#ifdef TEST_EXPAND_XY
-       for (c = 0; c < comp; c++) {
-               darray_t *temp;
-               int dimensions[3] = {w, h, k};
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
-               dimensions[0] = (w+1)>>1;
-               dimensions[1] = (h+1)>>1;
-               temp = darrayCreate(3, dimensions);
-               ASSERT(temp, "no memory");
-               rc = shrink_xy(values[0], temp);
-               ASSERT(rc >= 0, "no memory");
-               expand_xy(temp, values[0]);
-               darrayDestroy(temp);
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-       }
-#endif
-
-#ifdef TEST_EXPAND_Z
-       for (c = 0; c < comp; c++) {
-               darray_t *temp;
-               int dimensions[3] = {w, h, k};
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(values[0]) + z*w*h, darrayGetPr(result) + z*w*h*comp + c*w*h, w*h*sizeof(double));
-               dimensions[2] = (k+1)>>1;
-               temp = darrayCreate(3, dimensions);
-               ASSERT(temp, "no memory");
-               rc = shrink_z(values[0], temp, flows_x[0], flows_y[0]);
-               ASSERT(rc >= 0, "no memory");
-               rc = expand_z(temp, values[0], flows_x[0], flows_y[0]);
-               ASSERT(rc >= 0, "no memory");
-               darrayDestroy(temp);
-               for (z = 0; z < k; z++)
-                       memcpy(darrayGetPr(result) + z*w*h*comp + c*w*h, darrayGetPr(values[0]) + z*w*h, w*h*sizeof(double));
-       }
-#endif
-
-       rc = 0;
-
-error:
-       /* free memory */
-       for (d = 0; d < max_depth; d++) {
-               if (nb_list)
-                       free(nb_list[d]);
-               if (marks)
-                       darrayDestroy(marks[d]);
-               if (values)
-                       darrayDestroy(values[d]);
-               if (flows_xy_x)
-                       darrayDestroy(flows_xy_x[d]);
-               if (flows_xy_y)
-                       darrayDestroy(flows_xy_y[d]);
-               if (flows_x)
-                       darrayDestroy(flows_x[d]);
-               if (flows_y)
-                       darrayDestroy(flows_y[d]);
-               if (flows_ix)
-                       darrayDestroy(flows_ix[d]);
-               if (flows_iy)
-                       darrayDestroy(flows_iy[d]);
-       }
-       free(nb_list);
-       free(marks);
-       free(values);
-       free(flows_xy_x);
-       free(flows_xy_y);
-       free(flows_x);
-       free(flows_y);
-       free(flows_ix);
-       free(flows_iy);
-       darrayDestroy(init);
-       darrayDestroy(luminance);
-       if (rc < 0)
-               darrayDestroy(result);
-
-       return (rc < 0) ? NULL : result;
-}
diff --git a/lib/multigrid.h b/lib/multigrid.h
new file mode 100644 (file)
index 0000000..1253e64
--- /dev/null
@@ -0,0 +1,19 @@
+#ifndef ___MULTIGRID_H___
+#define ___MULTIGRID_H___
+/* print error message b, if a is false */
+#define ASSERT(a, b) \
+{ if (!(a)) { fprintf(stderr, "%s %s() Error in line %d: %s\n", __FILE__, __func__, __LINE__, (b)); abort(); } }
+
+#define min(a, b) ((a < b) ? a : b)
+#define max(a, b) ((a > b) ? a : b)
+
+int shrink_z(darray_t *src, darray_t *dst, const darray_t *flow);
+int shrink_flow_z(darray_t *src, darray_t *dst);
+int shrink_xy(darray_t *src, darray_t *dst);
+void quantize(const darray_t *array, double threshold, double gt, double le);
+void multiply(const darray_t *array, double factor);
+unsigned char *gen_neighbor(int w, int h, int k, const darray_t *flow, const darray_t *flow_i);
+void weighten(const darray_t *luminance, unsigned char *neighbors);
+int solve_mg(int nlevel, int iters, int ncycle, darray_t **values, darray_t **marks, const darray_t *init, unsigned char **nb_list, darray_t **flows, int scalexyz);
+#endif
+
index a96e607..d7f14c7 100644 (file)
@@ -406,17 +406,14 @@ static void print_test_help()
 
 int main(int argc, char *argv[])
 {
-       darray_t *gI = NULL, *cI = NULL, *markIm = NULL, *ntscIm = NULL, *resultIm;
-       darray_t *dx = NULL, *dy = NULL, *idx = NULL, *idy = NULL;
-#ifdef WITH_OPENCV
-       int *flow_map_x = NULL, *flow_map_y = NULL;
-#endif
+       darray_t *gI = NULL, *cI = NULL, *markIm = NULL, *ntscIm = NULL;
+       darray_t *flow = NULL, *flow_i = NULL;
        double *ptr, *ptr2, *ptr3;
        double diff, sum;
        int dims[4];
-       int width = 0, height = 0, w, h, k = 1, index = 0, count;
+       int w = 0, h = 0, load_w, load_h, k = 1, index = 0, z;
        unsigned char *img_buffer, *img_buffer_all = NULL, *mark_buffer = NULL;
-       int rc, i, ii, j;
+       int rc, x, y, i, ii;
        time_t start, end;
        const char* filename;
        char first_filename[256];
@@ -461,9 +458,9 @@ next_sequence:
        time(&start);
 
        change_bc = 0;
-       for (count = 0; count < k; count++) {
+       for (z = 0; z < k; z++) {
                if (sequence) {
-                       filename = sequence[count].filename;
+                       filename = sequence[z].filename;
                        /* first_filename is set by load_sequence */
                } else {
                        filename = argv[1];
@@ -482,23 +479,23 @@ next_sequence:
                }
 #endif
                // load original image and convert their RGB components to double RGB array
-               rc = load_img(-1, &img_buffer, &w, &h, filename, index + count);
+               rc = load_img(-1, &img_buffer, &load_w, &load_h, filename, index + z);
                if (rc) {
                        fprintf(stderr, "Failed to load grey image\n");
                        return 0;
                }
-               scale_img(img_buffer, w, h, scale);
-               if (count == 0) {
-                       width = w / scale;
-                       height = h / scale;
+               scale_img(img_buffer, load_w, load_h, scale);
+               if (z == 0) {
+                       w = load_w / scale;
+                       h = load_h / scale;
                }
-               if (w/scale != width || h/scale != height) {
+               if (load_w/scale != w || load_h/scale != h) {
                        fprintf(stderr, "Error: All input images must have equal dimenstions.\n");
                        return 0;
                }
                // now we know the dimensions, so we can create input arrays
                if (!gI) {
-                       dims[0] = width; dims[1] = height; dims[2] = 3; dims[3] = k;
+                       dims[0] = w; dims[1] = h; dims[2] = 3; dims[3] = k;
                        gI = darrayCreate(4, dims);
                }
                if (!gI) {
@@ -506,24 +503,24 @@ next_sequence:
                        exit (0);
                }
                if (!cI) {
-                       dims[0] = width; dims[1] = height; dims[2] = 3; dims[3] = k;
+                       dims[0] = w; dims[1] = h; dims[2] = 3; dims[3] = k;
                        cI = darrayCreate(4, dims);
                }
                if (!cI) {
                        printf("failed to create marked image array\n");
                        exit (0);
                }
-               ptr = darrayGetPr(gI) + width*height*3*count;
-               img2array(img_buffer, width, height, ptr, width, height);
+               ptr = darrayGetPr(gI) + w*h*3*z;
+               img2array(img_buffer, w, h, ptr, w, h);
 #ifdef WITH_OPENCV
                if (k > 1 && flow_enable) {
                        if (!img_buffer_all)
-                               img_buffer_all = malloc(width*height*3*k);
+                               img_buffer_all = malloc(w*h*3*k);
                        if (!img_buffer_all) {
                                printf("failed to create grey image array\n");
                                exit (0);
                        }
-                       memcpy(img_buffer_all + width*height*3*count, img_buffer, width*height*3);
+                       memcpy(img_buffer_all + w*h*3*z, img_buffer, w*h*3);
                }
 #endif
                free(img_buffer);
@@ -532,10 +529,10 @@ next_sequence:
                        char name[256];
                        unsigned char c;
                        // load marked mask and convert their RGB components to double YUV array
-                       memcpy(darrayGetPr(cI) + width*height*3*count, darrayGetPr(gI) + width*height*3*count, width*height*3 * sizeof(double));
+                       memcpy(darrayGetPr(cI) + w*h*3*z, darrayGetPr(gI) + w*h*3*z, w*h*3 * sizeof(double));
                        /* add extra memory for unscaled data to prevent buffer overflow */
                        if (!mark_buffer)
-                               mark_buffer = (unsigned char *)malloc(width*height*k + w*h);
+                               mark_buffer = (unsigned char *)malloc(w*h*k + load_w*load_h);
                        if (!mark_buffer) {
                                printf("no memory for mark buffer\n");
                                exit (0);
@@ -544,15 +541,15 @@ next_sequence:
                                printf("failed to load palette for file: '%s'\n", filename);
                                exit (0);
                        }
-                       sprintf(name, filename, index + count);
+                       sprintf(name, filename, index + z);
                        /* always load full unscaled image, then scale down */
-                       if (load_marked(mark_buffer + width*height*count, w, h, name) == 0) {
-                               scale_mark(mark_buffer + width*height*count, w, h, scale);
-                               ptr = darrayGetPr(cI) + width*height*3*count;
-                               for (i = 0; i < height; i++) {
-                                       for (j = 0; j < width; j++) {
+                       if (load_marked(mark_buffer + w*h*z, load_w, load_h, name) == 0) {
+                               scale_mark(mark_buffer + w*h*z, load_w, load_h, scale);
+                               ptr = darrayGetPr(cI) + w*h*3*z;
+                               for (y = 0; y < h; y++) {
+                                       for (x = 0; x < w; x++) {
                                                /* do not apply mask on index 0 */
-                                               c = mark_buffer[i*width+j + width*height*count];
+                                               c = mark_buffer[y*w+x + w*h*z];
                                                if (c == 0)
                                                        continue;
                                                /* check for any brightness/contrast change */
@@ -561,62 +558,62 @@ next_sequence:
                                                /* do not apply white pixles, this meas: keep original color */
                                                if (mark_palette[c-1].r == 255 && mark_palette[c-1].g == 255 && mark_palette[c-1].b == 255)
                                                        continue;
-                                               ptr[i*width+j] = mark_palette[c-1].r / 255.0F;
-                                               ptr[i*width+j + width*height] = mark_palette[c-1].g / 255.0F;
-                                               ptr[i*width+j + width*height*2] = mark_palette[c-1].b / 255.0F;
+                                               ptr[y*w+x] = mark_palette[c-1].r / 255.0F;
+                                               ptr[y*w+x + w*h] = mark_palette[c-1].g / 255.0F;
+                                               ptr[y*w+x + w*h*2] = mark_palette[c-1].b / 255.0F;
                                        }
                                }
                        } else
-                               memset(mark_buffer + width*height*count, 0, width*height);
+                               memset(mark_buffer + w*h*z, 0, w*h);
                } else {
                        // load marked image and convert their RGB components to double YUV array
-                       rc = load_img(-1, &img_buffer, &w, &h, argv[2], index + count);
+                       rc = load_img(-1, &img_buffer, &load_w, &load_h, argv[2], index + z);
                        if (!rc) {
-                               if (w/scale != width || h/scale != height) {
+                               if (load_w/scale != w || load_h/scale != h) {
                                        fprintf(stderr, "Error: All input images must have equal dimenstions.\n");
                                        return 0;
                                }
-                               ptr = darrayGetPr(cI) + width*height*3*count;
-                               img2array(img_buffer, width, height, ptr, width, height);
+                               ptr = darrayGetPr(cI) + w*h*3*z;
+                               img2array(img_buffer, w, h, ptr, w, h);
                                free(img_buffer);
                        } else {
                                fprintf(stderr, "Failed to load marked image, omitting...\n");
-                               memcpy(darrayGetPr(cI) + width*height*3*count, darrayGetPr(gI) + width*height*3*count, width*height*3 * sizeof(double));
+                               memcpy(darrayGetPr(cI) + w*h*3*z, darrayGetPr(gI) + w*h*3*z, w*h*3 * sizeof(double));
                        }
                }
        }
 
        /* create color mask and ntsc arrays for the colorization process */
-       dims[0] = width; dims[1] = height; dims[2] = k;
+       dims[0] = w; dims[1] = h; dims[2] = k;
        markIm = darrayCreate(3, dims);
        if (!markIm) {
                printf("failed to create mark array\n");
                exit (0);
        }
        features = (change_bc) ? 4 : 2;
-       dims[0] = width; dims[1] = height; dims[2] = features+1; dims[3] = k;
+       dims[0] = w; dims[1] = h; dims[2] = features+1; dims[3] = k;
        ntscIm = darrayCreate(4, dims);
        if (!ntscIm) {
                printf("failed to create ntsc array\n");
                exit (0);
        }
 
-       for (count = 0; count < k; count++) {
+       for (z = 0; z < k; z++) {
                if (sequence || !strcmp(argv[2], "marked")) {
                        unsigned char c;
-                       ptr = darrayGetPr(markIm) + width*height*count;
+                       ptr = darrayGetPr(markIm) + w*h*z;
                        // use marked mask to fill markIm
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
-                                       if (j < width && i < height) {
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
+                                       if (x < w && y < h) {
                                                /* do not apply mask on index 0 */
-                                               c = mark_buffer[i*width+j + width*height*count];
+                                               c = mark_buffer[y*w+x + w*h*z];
                                        } else
                                                c = 0;
                                        if (c)
-                                               ptr[i*width+j] = 1.0F;
+                                               ptr[y*w+x] = 1.0F;
                                        else
-                                               ptr[i*width+j] = 0.0F;
+                                               ptr[y*w+x] = 0.0F;
                                }
                        }
                } else {
@@ -626,10 +623,10 @@ next_sequence:
                        // - sum all components to get grey image
                        // - apply threshold (pixle is 1F, if the absolute difference is > 0.01F)
                        // original code: markIm=(sum(abs(gI-cI),3)>0.01);
-                       ptr = darrayGetPr(gI) + width*height*3*count;
-                       ptr2 = darrayGetPr(cI) + width*height*3*count;
-                       ptr3 = darrayGetPr(markIm) + width*height*count;
-                       for (i = 0, ii = width * height; i < ii; i++) {
+                       ptr = darrayGetPr(gI) + w*h*3*z;
+                       ptr2 = darrayGetPr(cI) + w*h*3*z;
+                       ptr3 = darrayGetPr(markIm) + w*h*z;
+                       for (i = 0, ii = w * h; i < ii; i++) {
                                diff = 0;
                                sum = ptr[i] - ptr2[i];
                                if (sum < 0)
@@ -654,142 +651,94 @@ next_sequence:
                }
 
                // convert original image into YUV
-               ptr = darrayGetPr(gI) + width*height*3*count;
-               rgb2yuv(ptr, ptr, width, height);
+               ptr = darrayGetPr(gI) + w*h*3*z;
+               rgb2yuv(ptr, ptr, w, h);
 
                // convert maked image into YUV
-               ptr = darrayGetPr(cI) + width*height*3*count;
-               rgb2yuv(ptr, ptr, width, height);
+               ptr = darrayGetPr(cI) + w*h*3*z;
+               rgb2yuv(ptr, ptr, w, h);
 
                if (test != BC_ONLY) {
                        if (sequence || !strcmp(argv[2], "marked")) {
                                unsigned char c;
                                // generate NTSC image: use luminance from original image and chrominance from original or marked image
-                               ptr = darrayGetPr(gI) + width*height*3*count;
-                               ptr2 = darrayGetPr(cI) + width*height*3*count;
-                               ptr3 = darrayGetPr(ntscIm) + width*height*(features+1)*count;
+                               ptr = darrayGetPr(gI) + w*h*3*z;
+                               ptr2 = darrayGetPr(cI) + w*h*3*z;
+                               ptr3 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
                                /* use original y component */
-                               memcpy(ptr3, ptr, width * height * sizeof(double));
+                               memcpy(ptr3, ptr, w * h * sizeof(double));
                                /* apply new uv components */
-                               for (i = 0; i < height; i++) {
-                                       for (j = 0; j < width; j++) {
-                                               c = mark_buffer[i*width+j + width*height*count];
+                               for (y = 0; y < h; y++) {
+                                       for (x = 0; x < w; x++) {
+                                               c = mark_buffer[y*w+x + w*h*z];
                                                if (c == 0) {
-                                                       ptr3[width * height + width * i + j] = ptr[width * height + width * i + j];
-                                                       ptr3[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i + j];
+                                                       ptr3[w * h + w * y + x] = ptr[w * h + w * y + x];
+                                                       ptr3[w * h * 2 + w * y + x] = ptr[w * h * 2 + w * y + x];
                                                } else {
-                                                       ptr3[width * height + width * i + j] = ptr2[width * height + width * i + j];
-                                                       ptr3[width * height * 2 + width * i + j] = ptr2[width * height * 2 + width * i + j];
+                                                       ptr3[w * h + w * y + x] = ptr2[w * h + w * y + x];
+                                                       ptr3[w * h * 2 + w * y + x] = ptr2[w * h * 2 + w * y + x];
                                                }
                                        }
                                }
                        } else {
                                // generate NTSC image: use luminance from original image and chrominance from maked image
-                               ptr = darrayGetPr(gI) + width*height*3*count;
-                               ptr2 = darrayGetPr(cI) + width*height*3*count;
-                               ptr3 = darrayGetPr(ntscIm) + width*height*(features+1)*count;
-                               memcpy(ptr3, ptr, width * height * sizeof(double));
-                               memcpy(ptr3 + width * height, ptr2 + width * height, width * height * sizeof(double));
-                               memcpy(ptr3 + width * height * 2, ptr2 + width * height * 2, width * height * sizeof(double));
+                               ptr = darrayGetPr(gI) + w*h*3*z;
+                               ptr2 = darrayGetPr(cI) + w*h*3*z;
+                               ptr3 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                               memcpy(ptr3, ptr, w * h * sizeof(double));
+                               memcpy(ptr3 + w * h, ptr2 + w * h, w * h * sizeof(double));
+                               memcpy(ptr3 + w * h * 2, ptr2 + w * h * 2, w * h * sizeof(double));
                        }
                } else {
                        /* use grey image as result if BC_ONLY test is selected */
-                       ptr = darrayGetPr(gI) + width*height*3*count;
-                       ptr2 = darrayGetPr(ntscIm) + width*height*(features+1)*count;
-                       memcpy(ptr2, ptr, width * height * 3 * sizeof(double));
+                       ptr = darrayGetPr(gI) + w*h*3*z;
+                       ptr2 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       memcpy(ptr2, ptr, w * h * 3 * sizeof(double));
                }
 
                if (change_bc) {
                        unsigned char c;
-                       ptr2 = darrayGetPr(ntscIm) + width*height*(features+1)*count;
+                       ptr2 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
                        /* use original y component */
-                       memcpy(ptr2, ptr, width * height * sizeof(double));
+                       memcpy(ptr2, ptr, w * h * sizeof(double));
                        /* apply brightness and contrast from makred pixles to uv components of grey image */
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
                                        /* use unchanged brightness and contrast on index 0 */
-                                       c = mark_buffer[i*width+j + width*height*count];
+                                       c = mark_buffer[y*w+x + w*h*z];
                                        if (c == 0) {
-                                               ptr2[i*width+j + width*height*3] = 0;
-                                               ptr2[i*width+j + width*height*4] = 0.1;
+                                               ptr2[y*w+x + w*h*3] = 0;
+                                               ptr2[y*w+x + w*h*4] = 0.1;
                                        } else {
-                                               ptr2[i*width+j + width*height*3] = mark_palette[c-1].bright / 10.0;
-                                               ptr2[i*width+j + width*height*4] = mark_palette[c-1].contrast / 10.0;
+                                               ptr2[y*w+x + w*h*3] = mark_palette[c-1].bright / 10.0;
+                                               ptr2[y*w+x + w*h*4] = mark_palette[c-1].contrast / 10.0;
                                        }
                                }
                        }
                }
        }
 
-       /* create flow vectors */
-       dims[0] = width; dims[1] = height; dims[2] = k - 1;
-       dx = darrayCreate(3, dims);
-       if (!dx) {
-               printf("failed to create array\n");
-               exit (0);
-       }
-       dy = darrayCreate(3, dims);
-       if (!dy) {
-               printf("failed to create array\n");
-               exit (0);
-       }
-       idx = darrayCreate(3, dims);
-       if (!idx) {
-               printf("failed to create array\n");
-               exit (0);
-       }
-       idy = darrayCreate(3, dims);
-       if (!idy) {
-               printf("failed to create array\n");
-               exit (0);
-       }
 #ifdef WITH_OPENCV
        if (k > 1 && flow_enable) {
-               printf("Calculating optical flow for %d frames: steps=%d win=%d max=%d\n", k, flow_steps/scale, flow_win/scale, flow_max/scale);
-               flow_map_x = (int *)malloc(sizeof(int)*width*height);
-               if (!flow_map_x) {
-                       printf("failed to alloc array\n");
+               /* create flow vectors */
+               dims[0] = w; dims[1] = h; dims[2] = k - 1; dims[3] = 2;
+               flow = darrayCreate(4, dims);
+               if (!flow) {
+                       printf("failed to create array\n");
                        exit (0);
                }
-               flow_map_y = (int *)malloc(sizeof(int)*width*height);
-               if (!flow_map_y) {
-                       printf("failed to alloc array\n");
+               flow_i = darrayCreate(4, dims);
+               if (!flow_i) {
+                       printf("failed to create array\n");
                        exit (0);
                }
+               printf("Calculating optical flow for %d frames: window=%d\n", k, flow_window/scale);
        }
-       for (count = 0; count < k-1; count++) {
-               if (flow_map_x && flow_map_y)
-                       create_flow_maps(NULL, img_buffer_all + width*height*3*(count+1), img_buffer_all + width*height*3*count, width, height, flow_steps/scale, flow_win/scale, flow_max/scale, NULL, NULL, flow_map_x, flow_map_y, NULL);
-               ptr = darrayGetPr(dx) + width*height*count;
-               if (flow_map_x) {
-                       for (j = 0; j < height; j++) {
-                               for (i = 0; i < width; i++)
-                                       ptr[i + width * j] = flow_map_x[width*j + i];
-                       }
-               }
-               ptr = darrayGetPr(dy) + width*height*count;
-               if (flow_map_y) {
-                       for (j = 0; j < height; j++) {
-                               for (i = 0; i < width; i++)
-                                       ptr[i + width * j] = flow_map_y[width*j + i];
-                       }
-               }
-               if (flow_map_x && flow_map_y)
-                       create_flow_maps(img_buffer_all + width*height*3*count, NULL, img_buffer_all + width*height*3*(count+1), width, height, flow_steps/scale, flow_win/scale, flow_max/scale, flow_map_x, flow_map_y, NULL, NULL, NULL);
-               ptr = darrayGetPr(idx) + width*height*count;
-               if (flow_map_x) {
-                       for (j = 0; j < height; j++) {
-                               for (i = 0; i < width; i++)
-                                       ptr[i + width * j] = flow_map_x[width*j + i];
-                       }
-               }
-               ptr = darrayGetPr(idy) + width*height*count;
-               if (flow_map_y) {
-                       for (j = 0; j < height; j++) {
-                               for (i = 0; i < width; i++)
-                                       ptr[i + width * j] = flow_map_y[width*j + i];
-                       }
-               }
+       for (z = 0; z < k-1; z++) {
+               if (flow)
+                       create_flow_maps(NULL, img_buffer_all + w*h*3*(z+1), img_buffer_all + w*h*3*z, w, h, flow_window/scale, 0, NULL, NULL, darrayGetPr(flow) + w*h*z, darrayGetPr(flow) + w*h*z + w*h*(k-1), NULL);
+               if (flow_i)
+                       create_flow_maps(img_buffer_all + w*h*3*z, NULL, img_buffer_all + w*h*3*(z+1), w, h, flow_window/scale, 0, darrayGetPr(flow_i) + w*h*z, darrayGetPr(flow_i) + w*h*z + w*h*(k-1), NULL, NULL, NULL);
        }
 #endif
        free(img_buffer_all);
@@ -802,67 +751,63 @@ next_sequence:
 
        if (test != FLOW_NEXT && test != FLOW_PREV && test != MARKED && test != MASK && test != MASK_COLOR && test != BC_ONLY && test != BC_IMAGE) {
                printf("Colorizing %d frames, please wait...\n", k);
-               resultIm = colorize(ntscIm, markIm, dx, dy, idx, idy, in_itr_num, out_itr_num, scalexyz);
-               if (!resultIm) {
+               rc = colorize(ntscIm, markIm, flow, flow_i, in_itr_num, out_itr_num, scalexyz);
+               if (rc < 0) {
                        if (k > 1)
                                printf("No memory! Use smaller frames or less frames between key frames or add more memory.");
                        else
                                printf("No memory! Use smaller image or add more memory.");
                        exit(-1);
                }
-       } else
-               resultIm = darrayClone(ntscIm);
+       }
 
-       /* if we have a change, we apply brightness+contrast from bcIm */
+       /* if we have a change, we apply brightness+contrast from ntscIm */
        if (change_bc) {
-               for (count = 0; count < k; count++) {
-                       ptr = darrayGetPr(ntscIm) + width*height*(features+1)*count;
-                       ptr2 = darrayGetPr(resultIm) + width*height*(features+1)*count;
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
+               for (z = 0; z < k; z++) {
+                       ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       ptr2 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
                                        /* apply contrast */
-                                       ptr2[width * i + j] = (ptr2[width * i + j] - 0.5) * ptr[width * height * 4 + width * i + j] * 10.0 + 0.5;
+                                       ptr2[w * y + x] = (ptr2[w * y + x] - 0.5) * ptr[w * h * 4 + w * y + x] * 10.0 + 0.5;
                                        /* apply brightness */
-                                       ptr2[width * i + j] += ptr[width * height * 3 + width * i + j] * 10.0;
-                                       if (ptr2[width * i + j] < 0)
-                                               ptr2[width * i + j] = 0;
-                                       if (ptr2[width * i + j] > 1)
-                                               ptr2[width * i + j] = 1;
+                                       ptr2[w * y + x] += ptr[w * h * 3 + w * y + x] * 10.0;
+                                       if (ptr2[w * y + x] < 0)
+                                               ptr2[w * y + x] = 0;
+                                       if (ptr2[w * y + x] > 1)
+                                               ptr2[w * y + x] = 1;
 #if 0
 #warning TEST: show brightness and contrast change as uv vectors on a grey array */
-ptr2[width * i + j] = 0.5;
-ptr2[width * height + width * i + j] = ptr[width * height + width * i + j] * 10;
-ptr2[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i + j] * 10 - 1;
+ptr2[w * y + x] = 0.5;
+ptr2[w * h + w * y + x] = ptr[w * h + w * y + x] * 10;
+ptr2[w * h * 2 + w * y + x] = ptr[w * h * 2 + w * y + x] * 10 - 1;
 #endif
                                }
                        }
                }
        }
 
-       darrayDestroy(ntscIm);
-       ntscIm = NULL;
-
 #ifdef WITH_OPENCV
        if (test == FLOW_NEXT || test == FLOW_PREV) {
                /* apply flow planes to result image as u and y vector */
-               for (count = 0; count < k; count++) {
-                       ptr = darrayGetPr(resultIm) + width*height*(features+1)*count;
+               for (z = 0; z < k; z++) {
+                       ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
                        if (test == FLOW_NEXT) {
-                               ptr2 = darrayGetPr(dx) + width*height*count;
-                               ptr3 = darrayGetPr(dy) + width*height*count;
+                               ptr2 = darrayGetPr(flow) + w*h*z;
+                               ptr3 = darrayGetPr(flow) + w*h*z*(k-1);
                        } else {
-                               ptr2 = darrayGetPr(idx) + width*height*count;
-                               ptr3 = darrayGetPr(idy) + width*height*count;
+                               ptr2 = darrayGetPr(flow_i) + w*h*z;
+                               ptr3 = darrayGetPr(flow_i) + w*h*z*(k-1);
                        }
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
-                                       ptr[width * i + j] = 0.5;
-                                       if (count < k-1) {
-                                               ptr[width * i + j + width*height] = ptr2[width * i + j] / flow_max;
-                                               ptr[width * i + j + width*height*2] = ptr3[width * i + j] / flow_max;
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
+                                       ptr[w * y + x] = 0.5;
+                                       if (z < k-1) {
+                                               ptr[w * y + x + w*h] = ptr2[w * y + x] / 50;
+                                               ptr[w * y + x + w*h*2] = ptr3[w * y + x] / 50;
                                        } else {
-                                               ptr[width * i + j + width*height] = 0;
-                                               ptr[width * i + j + width*height*2] = 0;
+                                               ptr[w * y + x + w*h] = 0;
+                                               ptr[w * y + x + w*h*2] = 0;
                                        }
                                }
                        }
@@ -872,14 +817,14 @@ ptr2[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i +
 
        if (test == MASK) {
                /* apply maked mask as image */
-               for (count = 0; count < k; count++) {
-                       ptr = darrayGetPr(resultIm) + width*height*(features+1)*count;
-                       ptr2 = darrayGetPr(markIm) + width*height*count;
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
-                                       ptr[width * i + j] = ptr2[width * i + j];
-                                       ptr[width * i + j + width*height] = ptr2[width * i + j];
-                                       ptr[width * i + j + width*height*2] = ptr2[width * i + j];
+               for (z = 0; z < k; z++) {
+                       ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       ptr2 = darrayGetPr(markIm) + w*h*z;
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
+                                       ptr[w * y + x] = ptr2[w * y + x];
+                                       ptr[w * y + x + w*h] = ptr2[w * y + x];
+                                       ptr[w * y + x + w*h*2] = ptr2[w * y + x];
                                }
                        }
                }
@@ -887,16 +832,16 @@ ptr2[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i +
 
        if (test == MASK_COLOR) {
                /* apply maked mask on grey image */
-               for (count = 0; count < k; count++) {
-                       ptr = darrayGetPr(resultIm) + width*height*(features+1)*count;
-                       ptr2 = darrayGetPr(markIm) + width*height*count;
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
+               for (z = 0; z < k; z++) {
+                       ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       ptr2 = darrayGetPr(markIm) + w*h*z;
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
                                        /* darken unmarked areas, make maked areas with uniformed brightness */
-                                       if (ptr2[width * i + j] < 0.5F)
-                                               ptr[width * i + j] = ptr[width * i + j] / 4;
+                                       if (ptr2[w * y + x] < 0.5F)
+                                               ptr[w * y + x] = ptr[w * y + x] / 4;
                                        else
-                                               ptr[width * i + j] = 0.5F;
+                                               ptr[w * y + x] = 0.5F;
                                }
                        }
                }
@@ -904,49 +849,49 @@ ptr2[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i +
 
        if (test == BC_IMAGE) {
                /* apply bc image as result image */
-               for (count = 0; count < k; count++) {
-                       ptr = darrayGetPr(resultIm) + width*height*(features+1)*count;
-                       ptr2 = darrayGetPr(resultIm) + width*height*(features+1)*count;
+               for (z = 0; z < k; z++) {
+                       ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                       ptr2 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
                        /* use uniformed brightness as y component */
-                       for (i = 0; i < height; i++) {
-                               for (j = 0; j < width; j++) {
-                                       ptr[width * i + j] = 0.5F;
+                       for (y = 0; y < h; y++) {
+                               for (x = 0; x < w; x++) {
+                                       ptr[w * y + x] = 0.5F;
                                }
                        }
                        if (change_bc) {
-                               ptr2 = darrayGetPr(resultIm) + width*height*(features+1)*count;
-                               memcpy(ptr+width*height, ptr2+width*height*3, width * height * sizeof(double));
-                               memcpy(ptr+width*height*2, ptr2+width*height*4, width * height * sizeof(double));
+                               ptr2 = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+                               memcpy(ptr+w*h, ptr2+w*h*3, w * h * sizeof(double));
+                               memcpy(ptr+w*h*2, ptr2+w*h*4, w * h * sizeof(double));
                        } else {
-                               memset(ptr+width*height, 0, width * height * sizeof(double));
-                               memset(ptr+width*height*2, 0, width * height * sizeof(double));
+                               memset(ptr+w*h, 0, w * h * sizeof(double));
+                               memset(ptr+w*h*2, 0, w * h * sizeof(double));
                        }
                }
        }
 
        // save result YUV array to image with RGB components
-       img_buffer = (unsigned char *)malloc(width * height * 3);
+       img_buffer = (unsigned char *)malloc(w * h * 3);
        if (!img_buffer) {
                fprintf(stderr, "Failed to allocate image buffer\n");
                return 0;
        }
-       for (count = 0; count < k; count++) {
-               ptr = darrayGetPr(resultIm) + width*height*(features+1)*count;
-               yuv2rgb(ptr, ptr, width, height);
-               array2img(ptr, width, height, img_buffer, width, height);
+       for (z = 0; z < k; z++) {
+               ptr = darrayGetPr(ntscIm) + w*h*(features+1)*z;
+               yuv2rgb(ptr, ptr, w, h);
+               array2img(ptr, w, h, img_buffer, w, h);
                if (sequence) {
                        static char name[256], *p, *q;
-                       p = sequence[count].filename;
+                       p = sequence[z].filename;
                        while((q = strchr(p, DIR_SEPERATOR)))
                                p = q + 1;
-                       strcpy(name, sequence[count].filename);
-                       name[p - sequence[count].filename] = '\0';
+                       strcpy(name, sequence[z].filename);
+                       name[p - sequence[z].filename] = '\0';
                        strcat(name, "colorized_");
                        strcat(name, p);
                        filename = name;
                } else
                        filename = argv[3];
-               save_img(img_buffer, width, height, filename, index + count);
+               save_img(img_buffer, w, h, filename, index + z);
        }
        free(img_buffer);
 
@@ -954,22 +899,12 @@ ptr2[width * height * 2 + width * i + j] = ptr[width * height * 2 + width * i +
        printf("Elapsed time: %d minutes, %d seconds\n", (int)(end-start)/60, (int)(end-start)%60);
 
        // destroy
-       darrayDestroy(resultIm);
-       resultIm = NULL;
-       darrayDestroy(dx);
-       dx = NULL;
-       darrayDestroy(dy);
-       dy = NULL;
-       darrayDestroy(idx);
-       idx = NULL;
-       darrayDestroy(idy);
-       idy = NULL;
-#ifdef WITH_OPENCV
-       free(flow_map_x);
-       flow_map_x = NULL;
-       free(flow_map_y);
-       flow_map_y = NULL;
-#endif
+       darrayDestroy(ntscIm);
+       ntscIm = NULL;
+       darrayDestroy(flow);
+       flow = NULL;
+       darrayDestroy(flow_i);
+       flow_i = NULL;
        darrayDestroy(markIm);
        markIm = NULL;
 
index dcd0f07..2de9f56 100644 (file)
@@ -1,18 +1,19 @@
 #include <stdio.h>
 #include "ppm.h"
 #include "opticalflow.h"
+#include "yuv.h"
 #include <opencv2/core/core_c.h>
 #include <opencv2/imgproc/imgproc_c.h>
 #include <opencv2/video/tracking.hpp>
 #include "dir_seperator.h"
 
-int flow_enable = 0, flow_steps, flow_win, flow_max;
+int flow_enable = 0, flow_window, flow_view_vector, flow_view_uv;
 
 void flow_default(void)
 {
-       flow_steps = 8;
-       flow_win = 30;
-       flow_max = 16;
+       flow_window = 15;
+       flow_view_vector = 16;
+       flow_view_uv = 0; /* disabled */
 }
 
 /* save flow setting */
@@ -28,9 +29,9 @@ void save_flow(const char *filename)
                return;
        }
        fprintf(fp, "flow_enable %d\n", flow_enable);
-       fprintf(fp, "flow_steps %d\n", flow_steps);
-       fprintf(fp, "flow_win %d\n", flow_win);
-       fprintf(fp, "flow_max %d\n", flow_max);
+       fprintf(fp, "flow_window %d\n", flow_window);
+       fprintf(fp, "flow_view_vector %d\n", flow_view_vector);
+       fprintf(fp, "flow_view_uv %d\n", flow_view_uv);
        fclose(fp);
 }
 
@@ -60,19 +61,19 @@ int load_flow(const char *filename)
                        sscanf(p, "%d", &flow_enable);
                        continue;
                }
-               if (!strncmp(line, "flow_steps ", 11)) {
-                       p = line + 11;
-                       sscanf(p, "%d", &flow_steps);
+               if (!strncmp(line, "flow_window ", 12)) {
+                       p = line + 12;
+                       sscanf(p, "%d", &flow_window);
                        continue;
                }
-               if (!strncmp(line, "flow_win ", 9)) {
-                       p = line + 9;
-                       sscanf(p, "%d", &flow_win);
+               if (!strncmp(line, "flow_view_vector ", 17)) {
+                       p = line + 17;
+                       sscanf(p, "%d", &flow_view_vector);
                        continue;
                }
-               if (!strncmp(line, "flow_max ", 9)) {
-                       p = line + 9;
-                       sscanf(p, "%d", &flow_max);
+               if (!strncmp(line, "flow_view_uv ", 13)) {
+                       p = line + 13;
+                       sscanf(p, "%d", &flow_view_uv);
                        continue;
                }
        }
@@ -80,156 +81,37 @@ int load_flow(const char *filename)
        return 0;
 }
 
-/* smooth image, it is safe if image pointer is NULL pointer */
-static IplImage *smooth_image(IplImage *image, int width, int height, int dim)
-{
-       IplImage *temp;
-       CvSize size;
-
-       if (!image)
-               return NULL;
-
-       size.width = width;
-       size.height = height;
-
-       temp = cvCreateImage(size, IPL_DEPTH_8U, 1);
-       if (!temp)
-               return image;
-
-       /* use odd dimenstions of convolution matrix */
-       cvSmooth(image, temp, CV_GAUSSIAN, dim | 1, dim | 1, 0, 0);
-       cvReleaseImage(&image);
-       return temp;
-}
-
-/* remove bad features and replace them by average flow */
-static int filter_max_flow(int num, char *features_found, CvPoint2D32f *corners_from, CvPoint2D32f *corners_to, int max_flow)
+/* draw flow matrix as lines */
+static void draw_lines_from_matrix(CvMat *flow, int width, int height, int steps, int color, IplImage *image)
 {
-       int i;
-       int n = 0;
-       int average_x = 0;
-       int average_y = 0;
-
-       for (i = 0; i < num; i++) {
-               if (features_found[i]) {
-                       if (max_flow
-                        && sqrt(pow(corners_from[i].x-corners_to[i].x, 2) +
-                                pow(corners_from[i].y-corners_to[i].y, 2)) > max_flow) {
-                               features_found[i] = 0;
-                               continue;
-                       }
-                       average_x += corners_to[i].x - corners_from[i].x;
-                       average_y += corners_to[i].y - corners_from[i].y;
-                       n++;
-               }
-       }
-
-       if (n > 1) {
-               average_x /= n;
-               average_y /= n;
-       }
-       /* replace bad vectors by average */
-       for (i = 0; i < num; i++) {
-               if (!features_found[i]) {
-                       corners_to[i].x = corners_from[i].x + average_x;
-                       corners_to[i].y = corners_from[i].y + average_y;
-               }
-       }
-       return n;
-}
-
-/* draw motion vectors as lines */
-static void draw_lines(int num, char *features_found, CvPoint2D32f *corners_from, CvPoint2D32f *corners_to, int color, IplImage *image)
-{
-       int i;
-
-       for (i = 0; i < num; i++) {
-               if (features_found[i]) {
-                       cvLine(image, cvPoint(cvRound(corners_from[i].x), cvRound(corners_from[i].y)), cvPoint(cvRound(corners_to[i].x), cvRound(corners_to[i].y)), CV_RGB(color,color,color), 1, 8, 0);
+       int flow_width = width / steps;
+       int flow_height = height / steps;
+       int i, j, x, y;
+       const float *f;
+
+       for (i = 0, y = steps/2; i < flow_height && y < height; i++, y += steps) {
+               f = (const float *)(flow->data.ptr + flow->step * y);
+               for (j = 0, x = steps/2; j < flow_width && x < width; j++, x += steps) {
+                       cvLine(image,
+                               cvPoint(cvRound(x), cvRound(y)),
+                               cvPoint(cvRound(f[2*x] + x), cvRound(f[2*x+1] + y)),
+                               CV_RGB(color,color,color),
+                               1, 8, 0);
                }
        }
 }
 
-/* take vectors to generate a map
- *
- * the vectors are linear interpolated to fill all pixles */
-static void scale_flow_map(int *step2_x, int *step2_y, CvPoint2D32f *corners_from, CvPoint2D32f *corners_to, int width, int height, int flow_width, int flow_height, int steps)
+/* take flow matrix to generate a map */
+static void copy_matrix_to_map(CvMat *flow, int width, int height, double *flow_map_x, double *flow_map_y)
 {
-       /* step 1: interpolate horizontally */
-       float step1_x[width*flow_height];
-       float step1_y[width*flow_height];
-       int i, j, k, pos;
-       float from_x, from_y, to_x, to_y;
-
-       /* interpolate horizontally */
-       for (i = 0; i < flow_height; i++) { // count vertical vectors
-               /* get first vector position */
-               pos = (width - (steps * (flow_width - 1))) / 2;
-               from_x = corners_to[i * flow_width + 0].x - corners_from[i * flow_width + 0].x;
-               from_y = corners_to[i * flow_width + 0].y - corners_from[i * flow_width + 0].y;
-               /* use left most vector to fill left of first vector */
-               for (k = 0; k < pos; k++) { // count pixles before first vector
-                       step1_x[i * width + k] = from_x;
-                       step1_y[i * width + k] = from_y;
-               }
-               /* interpolate between vectors */
-               for (j = 0; j < flow_width - 1; j++) { // count horizontal vectors
-                       to_x = corners_to[i * flow_width + j+1].x - corners_from[i * flow_width + j+1].x;
-                       to_y = corners_to[i * flow_width + j+1].y - corners_from[i * flow_width + j+1].y;
-                       for (k = 0; k < steps; k++) { // count from vector to next vector
-                               step1_x[i * width + pos] = from_x * (steps-k) / steps  +  to_x * k / steps;
-                               step1_y[i * width + pos] = from_y * (steps-k) / steps  +  to_y * k / steps;
-                               pos++;
-                       }
-                       /* shift vector for next loop */
-                       from_x = to_x;
-                       from_y = to_y;
-               }
-               /* use right most vector to fill right of last vector */
-               for (k = pos; k < width; k++) { // count pixles after last vector
-                       step1_x[i * width + k] = from_x;
-                       step1_y[i * width + k] = from_y;
-                       pos++;
-               }
-               if (pos != width) {
-                       printf("pos != with\n");
-                       abort();
-               }
-       }
-
-       /* interpolate vertically */
-       for (i = 0; i < width; i++) { // count horizontal vectors
-               /* get first vector position */
-               pos = (height - (steps * (flow_height - 1))) / 2;
-               from_x = step1_x[0 * width + i];
-               from_y = step1_y[0 * width + i];
-               /* use top most vector to fill above first vector */
-               for (k = 0; k < pos; k++) { // count pixles above first vector
-                       step2_x[k * width + i] = from_x;
-                       step2_y[k * width + i] = from_y;
-               }
-               /* interpolate between vectors */
-               for (j = 0; j < flow_height - 1; j++) { // count vertical vectors
-                       to_x = step1_x[(j+1) * width + i];
-                       to_y = step1_y[(j+1) * width + i];
-                       for (k = 0; k < steps; k++) { // count from vector to next vector
-                               step2_x[pos * width + i] = from_x * (steps-k) / steps  +  to_x * k / steps;
-                               step2_y[pos * width + i] = from_y * (steps-k) / steps  +  to_y * k / steps;
-                               pos++;
-                       }
-                       /* shift vector for next loop */
-                       from_x = to_x;
-                       from_y = to_y;
-               }
-               /* use bottom vector to fill below of last vector */
-               for (k = pos; k < height; k++) { // count pixles after last vector
-                       step2_x[pos * width + i] = from_x;
-                       step2_y[pos * width + i] = from_y;
-                       pos++;
-               }
-               if (pos != height) {
-                       printf("pos != height\n");
-                       abort();
+       int x, y;
+       const float *f;
+
+       for (y = 0; y < height; y++) {
+               f = (const float *)(flow->data.ptr + flow->step * y);
+               for (x = 0; x < width; x++) {
+                       *flow_map_x++ = *f++;
+                       *flow_map_y++ = *f++;
                }
        }
 }
@@ -250,21 +132,27 @@ static void scale_flow_map(int *step2_x, int *step2_y, CvPoint2D32f *corners_fro
  * flow_map_y_next = flow map for delta y of next image (result is stored here)
  *
  */
-void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char *img_next_buffer, unsigned char *img_buffer, int width, int height, int steps, int win_size, int max_flow, int *flow_map_x_prev, int *flow_map_y_prev, int *flow_map_x_next, int *flow_map_y_next, void *_image_preview)
+void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char *img_next_buffer, unsigned char *img_buffer, int width, int height, int win_size, int steps, double *flow_map_x_prev, double *flow_map_y_prev, double *flow_map_x_next, double *flow_map_y_next, void *_image_preview)
 {
        IplImage *image_prev = NULL, *image = NULL, *image_next = NULL, *image_preview = _image_preview;
-       IplImage *pyr_a = NULL, *pyr_b = NULL;
-       CvSize size, pyr_sz;
-       int flow_width = width / steps;
-       int flow_height = height / steps;
-       CvPoint2D32f corners_input[flow_width * flow_height];
-       CvPoint2D32f corners_prev[flow_width * flow_height];
-       CvPoint2D32f corners_next[flow_width * flow_height];
-       char features_found_prev[flow_width * flow_height];
-       char features_found_next[flow_width * flow_height];
-       int input_n;
+       CvMat *flow;
+       CvSize size;
        int i, j;
 
+       const float pyrScale = 0.5;
+       const float levels = 3;
+       const float winsize = win_size;
+       const float iterations = 8;
+       const float polyN = 5;
+       const float polySigma = 1.2;
+       const int flags = 0;
+
+       flow = cvCreateMat(height, width, CV_32FC2);
+       if (!flow) {
+               fprintf(stderr, "failed to allocate opencv matrix\n");
+               return NULL;
+       }
+               
        /* apply images to opencv image */
        size.width = width;
        size.height = height;
@@ -305,59 +193,25 @@ void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char
                }
        }
 
-#if 1
-       /* smooth image, it is safe if image pointer is NULL pointer */
-       image_prev = smooth_image(image_prev, width, height, steps);
-       image = smooth_image(image, width, height, steps);
-       image_next = smooth_image(image_next, width, height, steps);
-       if (image_preview)
-               image_preview = smooth_image(image_preview, width, height, steps);
-#endif
-
        /* optical flow and compose previous and next image with current image */
-       pyr_sz.width = image->width;// + 8;
-       pyr_sz.height = image->height;// / 3;
-
-       pyr_a = cvCreateImage(pyr_sz, IPL_DEPTH_8U, 1 );
-       pyr_b = cvCreateImage(pyr_sz, IPL_DEPTH_8U, 1 );
-
        if (image) {
-               /* define corners */
-               input_n = 0;
-               if ((image_prev || image_next)) {
-                       for (i = 0; i < flow_height; i++) {
-                               for (j = 0; j < flow_width; j++) {
-                                       corners_input[input_n].x = (width - (steps * (flow_width - 1))) / 2 + (j * steps);
-                                       corners_input[input_n].y = (height - (steps * (flow_height - 1))) / 2 + (i * steps);
-                                       input_n++;
-                               }
-                       }
-               }
-
                /* calc flow */
                if (image_prev) {
-                       cvCalcOpticalFlowPyrLK(image, image_prev, pyr_a, pyr_b, corners_input, corners_prev, input_n, cvSize(win_size | 1, win_size | 1), 5, features_found_prev, NULL, cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, 0.3), 0);
-                       filter_max_flow(input_n, features_found_prev, corners_input, corners_prev, max_flow);
-                       if (image_preview)
-                               draw_lines(input_n, features_found_prev, corners_input, corners_prev, 255, image_preview);
+                       cvCalcOpticalFlowFarneback(image, image_prev, flow, pyrScale, levels, winsize, iterations, polyN, polySigma, flags);
+                       if (image_preview && steps > 1)
+                               draw_lines_from_matrix(flow, width, height, steps, 255, image_preview);
                        if (flow_map_x_prev && flow_map_y_prev)
-                               scale_flow_map(flow_map_x_prev, flow_map_y_prev, corners_input, corners_prev, width, height, flow_width, flow_height, steps);
+                               copy_matrix_to_map(flow, width, height, flow_map_x_prev, flow_map_y_prev);
                }
                if (image_next) {
-                       cvCalcOpticalFlowPyrLK(image, image_next, pyr_a, pyr_b, corners_input, corners_next, input_n, cvSize(win_size | 1, win_size | 1), 5, features_found_next, NULL, cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, 0.3), 0);
-                       filter_max_flow(input_n, features_found_next, corners_input, corners_next, max_flow);
-                       if (image_preview)
-                               draw_lines(input_n, features_found_next, corners_input, corners_next, 192, image_preview);
+                       cvCalcOpticalFlowFarneback(image, image_next, flow, pyrScale, levels, winsize, iterations, polyN, polySigma, flags);
+                       if (image_preview && steps > 1)
+                               draw_lines_from_matrix(flow, width, height, steps, 192, image_preview);
                        if (flow_map_x_next && flow_map_y_next)
-                               scale_flow_map(flow_map_x_next, flow_map_y_next, corners_input, corners_next, width, height, flow_width, flow_height, steps);
+                               copy_matrix_to_map(flow, width, height, flow_map_x_next, flow_map_y_next);
                }
        }
 
-       if (pyr_a)
-               cvReleaseImage(&pyr_a);
-       if (pyr_b)
-               cvReleaseImage(&pyr_b);
-
        if (image)
                cvReleaseImage(&image);
        if (image_prev)
@@ -365,6 +219,8 @@ void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char
        if (image_next)
                cvReleaseImage(&image_next);
 
+       cvReleaseMat(&flow);
+
        return image_preview;
 }
 
@@ -375,18 +231,19 @@ void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char
  * filename_next = name of next image
  * img_buffer = current image buffer
  * width, height = dimensions of image buffer
- * steps = pixles between each feature = size of grid to track
  * win_size = flow algorithm window size
- * max_flow = maximum length of allowed flow vector (to avoid errors)
+ * steps = pixles between each feature to visualize vectors (1 to disable)
+ * uv_scale = apply vector (pointing to next plane) by applying uv-color (0 to disable)
  *
  */
-void create_flow_view(const char *filename_prev, const char *filename_next, unsigned char *img_buffer, int width, int height, int steps, int win_size, int max_flow)
+void create_flow_view(const char *filename_prev, const char *filename_next, unsigned char *img_buffer, int width, int height, int win_size, int steps, int uv_scale)
 {
        unsigned char *img_prev_buffer = NULL, *img_next_buffer = NULL;
-       int *flow_map_x_prev = NULL, *flow_map_y_prev = NULL, *flow_map_x_next = NULL, *flow_map_y_next = NULL;
+       double *flow_map_x_prev = NULL, *flow_map_y_prev = NULL, *flow_map_x_next = NULL, *flow_map_y_next = NULL;
        IplImage *image_preview = NULL;
        CvSize size;
-       unsigned char c, r, g, b;
+       unsigned char c;
+       double r, g, b;
        int w, h, rc;
        int i, j;
 
@@ -426,13 +283,13 @@ void create_flow_view(const char *filename_prev, const char *filename_next, unsi
                image_preview = cvCreateImage(size, IPL_DEPTH_8U, 1);
        }
 
-       flow_map_x_prev = (int *)malloc(sizeof(int)*width*height);
-       flow_map_y_prev = (int *)malloc(sizeof(int)*width*height);
-       flow_map_x_next = (int *)malloc(sizeof(int)*width*height);
-       flow_map_y_next = (int *)malloc(sizeof(int)*width*height);
+       flow_map_x_prev = (double *)malloc(sizeof(double)*width*height);
+       flow_map_y_prev = (double *)malloc(sizeof(double)*width*height);
+       flow_map_x_next = (double *)malloc(sizeof(double)*width*height);
+       flow_map_y_next = (double *)malloc(sizeof(double)*width*height);
 
        if (image_preview)
-               image_preview = create_flow_maps(img_prev_buffer, img_next_buffer, img_buffer, width, height, steps, win_size, max_flow, flow_map_x_prev, flow_map_y_prev, flow_map_x_next, flow_map_y_next, image_preview);
+               image_preview = create_flow_maps(img_prev_buffer, img_next_buffer, img_buffer, width, height, win_size, steps, flow_map_x_prev, flow_map_y_prev, flow_map_x_next, flow_map_y_next, image_preview);
 
        if (img_buffer && image_preview) {
                /* paint vectors from image_preview to img_buffer */
@@ -440,30 +297,31 @@ void create_flow_view(const char *filename_prev, const char *filename_next, unsi
                        for (j = 0; j < width; j++) {
                                c = CV_IMAGE_ELEM(image_preview, uchar, i, j);
                                if (c > 224) {
-                                       r = 255;
+                                       r = 1;
                                        g = 0;
                                        b = 0;
                                } else if (c > 160) {
                                        r = 0;
-                                       g = 255;
+                                       g = 1;
                                        b = 0;
                                } else {
-#if 0
-                                       if (img_next_buffer) {
-                                               r = g = b = c / 2 + 64;
-                                               r += flow_map_x_next[width * i + j] * 64 / max_flow;
-                                               b += flow_map_y_next[width * i + j] * 64 / max_flow;
-                                       } else
-#endif
-                                       {
-                                               r = c;
-                                               g = c;
-                                               b = c;
+                                       if (!steps && uv_scale)
+                                               c += 64;
+                                       if (img_next_buffer && uv_scale) {
+                                               yuv2rgb_pixle(
+                                                       (double)c / 255.0,
+                                                       flow_map_x_next[width * i + j] / uv_scale,
+                                                       flow_map_y_next[width * i + j] / uv_scale,
+                                                       &r, &g, &b);
+                                       } else {
+                                               r = (double)c / 255.0;
+                                               g = (double)c / 255.0;
+                                               b = (double)c / 255.0;
                                        }
                                }
-                               img_buffer[(i*width+j)*3 + 0] = r;
-                               img_buffer[(i*width+j)*3 + 1] = g;
-                               img_buffer[(i*width+j)*3 + 2] = b;
+                               img_buffer[(i*width+j)*3 + 0] = r * 255.0;
+                               img_buffer[(i*width+j)*3 + 1] = g * 255.0;
+                               img_buffer[(i*width+j)*3 + 2] = b * 255.0;
                        }
                }
        }
index f1dae7e..8d7cf91 100644 (file)
@@ -1,9 +1,9 @@
-extern int flow_enable, flow_steps, flow_win, flow_max;
+extern int flow_enable, flow_window, flow_levels, flow_iter, flow_view_vector, flow_view_uv;
 
 void flow_default(void);
 void save_flow(const char *filename);
 int load_flow(const char *filename);
 
-void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char *img_next_buffer, unsigned char *img_buffer, int width, int height, int steps, int win_size, int max_flow, int *flow_map_x_prev, int *flow_map_y_prev, int *flow_map_x_next, int *flow_map_y_next, void *_image_preview);
+void *create_flow_maps(const unsigned char *img_prev_buffer, const unsigned char *img_next_buffer, unsigned char *img_buffer, int width, int height, int win_size, int steps, double *flow_map_x_prev, double *flow_map_y_prev, double *flow_map_x_next, double *flow_map_y_next, void *_image_preview);
 
-void create_flow_view(const char *filename_prev, const char *filename_next, unsigned char *img_buffer, int width, int height, int steps, int win_size, int max_flow);
+void create_flow_view(const char *filename_prev, const char *filename_next, unsigned char *img_buffer, int width, int height, int win_size, int steps, int uv_scale);
index 75dad67..e3c63c6 100644 (file)
--- a/version.h
+++ b/version.h
@@ -1 +1 @@
-"0.1"
+"0.2"