+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include "darray.h"
+#include "multigrid.h"
+#include "colorize.h"
+
+/* colorize image array using image_mark
+ *
+ * image:
+ * first dimansion = width
+ * second dimension = height
+ * thrid dimension = components (first component is luminance, other components are features, like u, v)
+ * fourth dimension = time
+ *
+ * image_mark:
+ * first dimansion = width
+ * second dimension = height
+ * fourth dimension = time
+ */
+int colorize_prepare(struct colorize *col)
+{
+ int dimensions[4];
+ int ww, hh, kk, last_kk, d, temp, z;
+ darray_t *luminance = NULL;
+ int rc = -1;
+
+ ASSERT(darrayGetNumberOfDimensions(col->image) == 4, "image array must have 4 dimensions");
+ col->w = darrayGetDimensions(col->image)[0];
+ col->h = darrayGetDimensions(col->image)[1];
+ col->comp = darrayGetDimensions(col->image)[2];
+ col->k = darrayGetDimensions(col->image)[3];
+ ASSERT(col->comp > 0, "image array requires at least one component (Y-component)");
+
+ ASSERT(darrayGetNumberOfDimensions(col->image_mark) == 3, "image_mark array must have 3 dimensions");
+ ASSERT(col->w == darrayGetDimensions(col->image_mark)[0], "image_mark and image must both have same x-size");
+ ASSERT(col->h == darrayGetDimensions(col->image_mark)[1], "image_mark and image must both have same y-size");
+ ASSERT(col->k == darrayGetDimensions(col->image_mark)[2], "image_mark and image must both have same z-size");
+
+ /* Calculate maximum depth
+ * Scale down dimensions until at least 2x2 remains
+ */
+ temp = min(col->w, col->h);
+ for (col->max_depth = 0; temp > 2; col->max_depth++)
+ temp=(temp+1)>>1;
+ ASSERT(col->max_depth >= 2, "image too small");
+
+ col->nb_list = calloc(col->max_depth, sizeof(struct neighbors *));
+ if (!col->nb_list)
+ goto error;
+ col->marks = calloc(col->max_depth, sizeof(darray_t *));
+ if (!col->marks)
+ goto error;
+ col->values = calloc(col->max_depth, sizeof(darray_t *));
+ if (!col->values)
+ goto error;
+ col->flows = calloc(col->max_depth, sizeof(darray_t *));
+ if (!col->flows)
+ goto error;
+ col->flows_i = calloc(col->max_depth, sizeof(darray_t *));
+ if (!col->flows_i)
+ goto error;
+ if (col->scalexyz) {
+ /* the flows_xy are scaled in xy in the current level */
+ col->flows_xy = calloc(col->max_depth, sizeof(darray_t *));
+ if (!col->flows_xy)
+ goto error;
+ }
+ dimensions[0] = col->w;
+ dimensions[1] = col->h;
+ dimensions[2] = col->k;
+ luminance = darrayCreate(3, dimensions);
+ if (!luminance)
+ goto error;
+ /* generate sparse multigrid */
+ last_kk = kk = col->k;
+ printf("-> prepare sparse grid");
+ for (d = 0, ww=col->w, hh=col->h; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
+ /* create temporary arrays
+ * copy fist array to level 0 or shrink previous level
+ */
+ printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
+ dimensions[0] = ww;
+ dimensions[1] = hh;
+ dimensions[2] = last_kk;
+ col->marks[d] = darrayCreate(3, dimensions);
+ if (!col->marks[d])
+ goto error;
+ dimensions[0] = ww;
+ dimensions[1] = hh;
+ dimensions[2] = last_kk-1;
+ dimensions[3] = 2;
+ col->flows[d] = darrayCreate(4, dimensions);
+ if (!col->flows[d])
+ goto error;
+ col->flows_i[d] = darrayCreate(4, dimensions);
+ if (!col->flows_i[d])
+ goto error;
+ if (d < col->scalexyz) {
+ dimensions[0] = (ww+1)>>1;
+ dimensions[1] = (hh+1)>>1;
+ dimensions[2] = kk-1;
+ dimensions[3] = 2;
+ col->flows_xy[d] = darrayCreate(4, dimensions);
+ if (!col->flows_xy[d])
+ goto error;
+ }
+ if (d == 0) {
+ for (z = 0; z < kk; z++) {
+ darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
+ darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
+ memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
+ }
+ if (col->flow && col->flow_i) {
+ darrayCheckPr(col->flows[d], darrayGetPr(col->flows[d]), col->w*col->h*(kk-1)*2);
+ darrayCheckPr(col->flow, darrayGetPr(col->flow), col->w*col->h*(kk-1)*2);
+ memcpy(darrayGetPr(col->flows[d]), darrayGetPr(col->flow), col->w*col->h*(kk-1)*2*sizeof(double));
+ darrayCheckPr(col->flows_i[d], darrayGetPr(col->flows_i[d]), col->w*col->h*(kk-1)*2);
+ darrayCheckPr(col->flow_i, darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2);
+ memcpy(darrayGetPr(col->flows_i[d]), darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2*sizeof(double));
+ }
+ darrayCheckPr(col->marks[d], darrayGetPr(col->marks[d]), col->w*col->h*kk);
+ darrayCheckPr(col->image_mark, darrayGetPr(col->image_mark), col->w*col->h*kk);
+ memcpy(darrayGetPr(col->marks[d]), darrayGetPr(col->image_mark), col->w*col->h*kk*sizeof(double));
+ /* invert, so 0 = marked */
+ quantize(col->marks[d], 0.5, 0, 1);
+ } else {
+ rc = shrink_xy(luminance, NULL);
+ if (rc < 0)
+ goto error;
+ rc = shrink_xy(col->marks[d-1], col->marks[d]);
+ if (rc < 0)
+ goto error;
+ rc = shrink_xy(col->flows[d-1], col->flows[d]);
+ if (rc < 0)
+ goto error;
+ rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
+ if (rc < 0)
+ goto error;
+ /* shrink length of flow vector on smaller grid.
+ * This is not performed when comparing with
+ * the original algorithm. (it is a bug there)
+ */
+ multiply(col->flows[d], 0.5);
+ multiply(col->flows_i[d], 0.5);
+ if (d-1 < col->scalexyz) { /* process if previous d shall be shrinked in z direction */
+ /* shrink time */
+ rc = shrink_z(luminance, NULL, col->flows[d]);
+ if (rc < 0)
+ goto error;
+ rc = shrink_z(col->marks[d], NULL, col->flows[d]);
+ if (rc < 0)
+ goto error;
+ rc = shrink_flow_z(col->flows[d], NULL);
+ if (rc < 0)
+ goto error;
+ rc = shrink_flow_z(col->flows_i[d], NULL);
+ if (rc < 0)
+ goto error;
+ }
+ /* If any of the restricted pixles is 0 (marked),
+ * the result must be 0 (marked) too.
+ */
+ quantize(col->marks[d], 0.9999, 1, 0);
+ }
+ if (d < col->scalexyz) {
+ /* create x and y shrinked flow vectors for restriction and prolonging process */
+ rc = shrink_xy(col->flows[d], col->flows_xy[d]);
+ if (rc < 0)
+ goto error;
+ multiply(col->flows_xy[d], 0.5);
+ }
+
+ /* generate values array */
+ dimensions[0] = ww;
+ dimensions[1] = hh;
+ dimensions[2] = kk;
+ col->values[d] = darrayCreate(3, dimensions);
+ if (!col->values[d])
+ goto error;
+
+ /* generate array of neighbors */
+#ifndef TEST_NO_RENDER
+ col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
+ if (!col->nb_list[d])
+ goto error;
+ /* weighten */
+ weighten(luminance, col->nb_list[d]);
+#endif
+
+ last_kk = kk;
+ if (d < col->scalexyz)
+ kk = (kk+1)>>1;
+ }
+ darrayDestroy(luminance);
+ luminance = NULL;
+ printf("\n");
+
+ dimensions[0] = col->w;
+ dimensions[1] = col->h;
+ dimensions[2] = col->k;
+ col->init = darrayCreate(3, dimensions);
+ if (!col->init)
+ goto error;
+
+#ifdef TEST_PYRAMID
+ dimensions[0] = col->w;
+ dimensions[1] = col->h;
+ dimensions[2] = col->k;
+ luminance = darrayCreate(3, dimensions);
+ for (z = 0; z < col->k; z++) {
+ darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
+ darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
+ memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
+ darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+ darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
+#ifdef TEST_FLOW
+ draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
+#endif
+#ifdef TEST_MARKED
+ paste(col->marks[0], col->values[0], 0, 0, 0);
+#endif
+
+ }
+ for (d = 1; d < col->max_depth; d++) {
+ rc = shrink_xy(luminance, NULL);
+ if (rc < 0)
+ goto error;
+ if (d-1 < col->scalexyz)
+ shrink_z(luminance, NULL, col->flows_xy[d-1]);
+ paste(luminance, col->values[0], 0, 0, 0);
+#ifdef TEST_FLOW
+ draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
+#endif
+#ifdef TEST_MARKED
+ paste(col->marks[d], col->values[0], 0, 0, 0);
+#endif
+ }
+ darrayDestroy(luminance);
+ luminance = NULL;
+ /* copy test pyramid image as Y-component */
+ for (z = 0; z < col->k; z++) {
+ darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
+ darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+#endif
+
+ rc = 0;
+
+error:
+ darrayDestroy(luminance);
+
+ return rc;
+}
+
+int colorize_solve(struct colorize *col, int quick)
+{
+ int z, d, c, ww, hh, kk;
+ int rc = -1;
+
+ if (!quick)
+ printf("-> apply component to grid and solve");
+ for (c = 1; c < col->comp; c++) {
+ int pixles, j;
+ double *value_ptr, *init_ptr, *mark_ptr;
+
+ printf(" #%d", c); fflush(stdout);
+ /* apply component
+ * - copy component into value array (level 0), otherwise 0
+ * - use maked colors for init array, otherwise use 0
+ */
+ for (z = 0; z < col->k; z++) {
+ darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+ darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
+ }
+ for (d = 1, ww=(col->w+1)>>1, hh=(col->h+1)>>1, kk=col->k; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
+ if (d-1 < col->scalexyz) /* scale kk, if last array was scaled */
+ kk=(kk+1)>>1;
+ value_ptr = darrayGetPr(col->values[d]);
+ pixles = ww*hh*kk;
+ for (j = 0; j < pixles; j++) {
+ darrayCheckPr(col->values[d], value_ptr, 1);
+ *value_ptr++ = 0.0;
+ }
+ }
+ value_ptr = darrayGetPr(col->values[0]);
+ init_ptr = darrayGetPr(col->init);
+ mark_ptr = darrayGetPr(col->image_mark);
+ pixles = col->w*col->h*col->k;
+ for (j = 0; j < pixles; j++) {
+ darrayCheckPr(col->image_mark, mark_ptr, 1);
+ if (*mark_ptr++ > 0.5) {
+ darrayCheckPr(col->init, init_ptr, 1);
+ darrayCheckPr(col->values[0], value_ptr, 1);
+ *init_ptr++ = *value_ptr;
+ } else {
+ darrayCheckPr(col->init, init_ptr, 1);
+ *init_ptr++ = 0.0;
+ }
+ value_ptr++;
+ }
+
+#ifndef TEST_NO_RENDER
+ /* colorize component */
+ if (quick) {
+ rc = solve_mg(col->max_depth, col->inner_iter, 1, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+ if (rc < 0)
+ goto error;
+ } else {
+ int i;
+ for (i = 0; i < col->outer_iter; i++) {
+ rc = solve_mg(col->max_depth, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+ if (rc < 0)
+ goto error;
+ rc = solve_mg((col->max_depth+1) / 2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+ if (rc < 0)
+ goto error;
+ rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+ if (rc < 0)
+ goto error;
+ rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
+ if (rc < 0)
+ goto error;
+ }
+ }
+#endif
+
+#ifdef TEST_PYRAMID
+ for (d = 1; d < col->max_depth; d++)
+ paste(col->values[d], col->values[0], 0, 0, 0);
+#endif
+
+ /* apply colorized component */
+ for (z = 0; z < col->k; z++) {
+ darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
+ darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+ }
+ if (!quick)
+ printf("\n");
+
+#ifdef TEST_SHRINK_XY
+ for (c = 0; c < col->comp; c++) {
+ darray_t *temp;
+ int dimensions[3] = {col->w, col->h, col->k};
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+ dimensions[0] = (col->w+1)>>1;
+ dimensions[1] = (col->h+1)>>1;
+ temp = darrayCreate(3, dimensions);
+ ASSERT(temp, "no memory");
+ rc = shrink_xy(col->values[0], temp);
+ ASSERT(rc >= 0, "no memory");
+ paste(temp, col->values[0], 0, 0, 0);
+ darrayDestroy(temp);
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+#endif
+
+#ifdef TEST_SHRINK_Z
+ for (c = 0; c < col->comp; c++) {
+ darray_t *temp;
+ int dimensions[3] = {col->w, col->h, col->k};
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+ dimensions[2] = (col->k+1)>>1;
+ temp = darrayCreate(3, dimensions);
+ ASSERT(temp, "no memory");
+ rc = shrink_z(col->values[0], temp, col->flows[0]);
+ ASSERT(rc >= 0, "no memory");
+ paste(temp, col->values[0], 0, 0, 0);
+ darrayDestroy(temp);
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+#endif
+
+#ifdef TEST_EXPAND_XY
+ for (c = 0; c < col->comp; c++) {
+ darray_t *temp;
+ int dimensions[3] = {col->w, col->h, col->k};
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+ dimensions[0] = (col->w+1)>>1;
+ dimensions[1] = (col->h+1)>>1;
+ temp = darrayCreate(3, dimensions);
+ ASSERT(temp, "no memory");
+ rc = shrink_xy(col->values[0], temp);
+ ASSERT(rc >= 0, "no memory");
+ expand_xy(temp, col->values[0]);
+ darrayDestroy(temp);
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+#endif
+
+#ifdef TEST_EXPAND_Z
+ for (c = 0; c < col->comp; c++) {
+ darray_t *temp;
+ int dimensions[3] = {col->w, col->h, col->k};
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, col->w*col->h*sizeof(double));
+ dimensions[2] = (col->k+1)>>1;
+ temp = darrayCreate(3, dimensions);
+ ASSERT(temp, "no memory");
+ rc = shrink_z(col->values[0], temp, col->flows[0]);
+ ASSERT(rc >= 0, "no memory");
+ rc = expand_z(temp, col->values[0], col->flows[0]);
+ ASSERT(rc >= 0, "no memory");
+ darrayDestroy(temp);
+ for (z = 0; z < col->k; z++)
+ memcpy(darrayGetPr(col->image) + z*col->w*col->h*col->comp + c*col->w*col->h, darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h*sizeof(double));
+ }
+#endif
+
+ rc = 0;
+
+error:
+ return rc;
+}
+
+void colorize_free(struct colorize *col)
+{
+ int d;
+
+ /* free memory */
+ for (d = 0; d < col->max_depth; d++) {
+ if (col->nb_list)
+ free(col->nb_list[d]);
+ if (col->marks)
+ darrayDestroy(col->marks[d]);
+ if (col->values)
+ darrayDestroy(col->values[d]);
+ if (col->flows_xy)
+ darrayDestroy(col->flows_xy[d]);
+ if (col->flows)
+ darrayDestroy(col->flows[d]);
+ if (col->flows_i)
+ darrayDestroy(col->flows_i[d]);
+ }
+ free(col->nb_list);
+ free(col->marks);
+ free(col->values);
+ free(col->flows_xy);
+ free(col->flows);
+ free(col->flows_i);
+ darrayDestroy(col->init);
+}
+
+int colorize(const darray_t *image, const darray_t *image_mark, const darray_t *flow, const darray_t *flow_i, int inner_iter, int outer_iter, int scalexyz)
+{
+ struct colorize col;
+ int rc = -1;
+
+ memset(&col, 0, sizeof(col));
+ col.image = image;
+ col.image_mark = image_mark;
+ col.flow = flow;
+ col.flow_i = flow_i;
+ col.inner_iter = inner_iter;
+ col.outer_iter = outer_iter;
+ col.scalexyz = scalexyz;
+
+ rc = colorize_prepare(&col);
+ if (rc < 0)
+ goto error;
+
+ rc = colorize_solve(&col, 0);
+ if (rc < 0)
+ goto error;
+
+ rc = 0;
+
+error:
+ colorize_free(&col);
+
+ return rc;
+}