From 9ad18bbdf7863fdb20cbb6c0523aab2dc436d6cc Mon Sep 17 00:00:00 2001 From: Andreas Eversberg Date: Wed, 30 Apr 2014 11:03:03 +0200 Subject: [PATCH] Version 0.2 - various fixes, cleanups and improvements --- gui/colorize.c | 38 ++- gui/flow.c | 104 +++---- gui/image.c | 55 +++- gui/image.h | 2 +- gui/main.c | 53 ++-- gui/main.h | 1 + gui/timeline.c | 6 +- lib/Makefile.am | 3 +- lib/colorize.c | 484 ++++++++++++++++++++++++++++++ lib/colorize.h | 28 +- lib/darray.c | 14 + lib/darray.h | 12 + lib/multigrid.c | 874 ++++++++++++++++-------------------------------------- lib/multigrid.h | 19 ++ src/colorize.c | 407 +++++++++++-------------- src/opticalflow.c | 336 ++++++--------------- src/opticalflow.h | 6 +- version.h | 2 +- 18 files changed, 1259 insertions(+), 1185 deletions(-) create mode 100644 lib/colorize.c create mode 100644 lib/multigrid.h diff --git a/gui/colorize.c b/gui/colorize.c index 2e1edf2..61d396d 100644 --- a/gui/colorize.c +++ b/gui/colorize.c @@ -10,9 +10,11 @@ #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); diff --git a/gui/flow.c b/gui/flow.c index e1f6740..2d70331 100644 --- a/gui/flow.c +++ b/gui/flow.c @@ -10,16 +10,16 @@ #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); } diff --git a/gui/image.c b/gui/image.c index 70f5cb7..1590279 100644 --- a/gui/image.c +++ b/gui/image.c @@ -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); diff --git a/gui/image.h b/gui/image.h index 2dca9b7..59b4b4a 100644 --- a/gui/image.h +++ b/gui/image.h @@ -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); diff --git a/gui/main.c b/gui/main.c index 00b76b4..71cff5f 100644 --- a/gui/main.c +++ b/gui/main.c @@ -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 */ diff --git a/gui/main.h b/gui/main.h index 8aaa1b3..58ff8c0 100644 --- a/gui/main.h +++ b/gui/main.h @@ -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; diff --git a/gui/timeline.c b/gui/timeline.c index 081113e..5a25b86 100644 --- a/gui/timeline.c +++ b/gui/timeline.c @@ -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) { diff --git a/lib/Makefile.am b/lib/Makefile.am index 187c279..f45e7b0 100644 --- a/lib/Makefile.am +++ b/lib/Makefile.am @@ -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 index 0000000..6c8caba --- /dev/null +++ b/lib/colorize.c @@ -0,0 +1,484 @@ +#include +#include +#include +#include +#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; +} diff --git a/lib/colorize.h b/lib/colorize.h index 228e3d6..295bf36 100644 --- a/lib/colorize.h +++ b/lib/colorize.h @@ -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 diff --git a/lib/darray.c b/lib/darray.c index adf205c..aed5911 100644 --- a/lib/darray.c +++ b/lib/darray.c @@ -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 + diff --git a/lib/darray.h b/lib/darray.h index f328824..c6d4789 100644 --- a/lib/darray.h +++ b/lib/darray.h @@ -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 */ diff --git a/lib/multigrid.c b/lib/multigrid.c index 20f3312..1687d07 100644 --- a/lib/multigrid.c +++ b/lib/multigrid.c @@ -3,7 +3,11 @@ #include #include #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 @@ -33,41 +37,24 @@ /* 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>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 index 0000000..1253e64 --- /dev/null +++ b/lib/multigrid.h @@ -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 + diff --git a/src/colorize.c b/src/colorize.c index a96e607..d7f14c7 100644 --- a/src/colorize.c +++ b/src/colorize.c @@ -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; diff --git a/src/opticalflow.c b/src/opticalflow.c index dcd0f07..2de9f56 100644 --- a/src/opticalflow.c +++ b/src/opticalflow.c @@ -1,18 +1,19 @@ #include #include "ppm.h" #include "opticalflow.h" +#include "yuv.h" #include #include #include #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; } } } diff --git a/src/opticalflow.h b/src/opticalflow.h index f1dae7e..8d7cf91 100644 --- a/src/opticalflow.h +++ b/src/opticalflow.h @@ -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); diff --git a/version.h b/version.h index 75dad67..e3c63c6 100644 --- a/version.h +++ b/version.h @@ -1 +1 @@ -"0.1" +"0.2" -- 2.13.6