9 /* colorize image array using image_mark
12 * first dimansion = width
13 * second dimension = height
14 * thrid dimension = components (first component is luminance, other components are features, like u, v)
15 * fourth dimension = time
18 * first dimansion = width
19 * second dimension = height
20 * fourth dimension = time
22 int colorize_prepare(struct colorize *col, int quick, int alloc)
25 int ww, hh, kk, last_kk, d, temp, z;
26 darray_t *luminance = NULL;
29 ASSERT(darrayGetNumberOfDimensions(col->image) == 4, "image array must have 4 dimensions");
30 col->w = darrayGetDimensions(col->image)[0];
31 col->h = darrayGetDimensions(col->image)[1];
32 col->comp = darrayGetDimensions(col->image)[2];
33 col->k = darrayGetDimensions(col->image)[3];
34 ASSERT(col->comp > 0, "image array requires at least one component (Y-component)");
36 ASSERT(darrayGetNumberOfDimensions(col->image_mark) == 3, "image_mark array must have 3 dimensions");
37 ASSERT(col->w == darrayGetDimensions(col->image_mark)[0], "image_mark and image must both have same x-size");
38 ASSERT(col->h == darrayGetDimensions(col->image_mark)[1], "image_mark and image must both have same y-size");
39 ASSERT(col->k == darrayGetDimensions(col->image_mark)[2], "image_mark and image must both have same z-size");
41 /* Calculate maximum depth
42 * Scale down dimensions until at least 2x2 remains
44 temp = min(col->w, col->h);
45 for (col->max_depth = 0; temp > 2; col->max_depth++)
47 ASSERT(col->max_depth >= 2, "image too small");
50 col->nb_list = calloc(col->max_depth, sizeof(struct neighbors *));
56 col->marks = calloc(col->max_depth, sizeof(darray_t *));
62 col->values = calloc(col->max_depth, sizeof(darray_t *));
68 col->flows = calloc(col->max_depth, sizeof(darray_t *));
74 col->flows_i = calloc(col->max_depth, sizeof(darray_t *));
80 /* the flows_xy are scaled in xy in the current level */
82 col->flows_xy = calloc(col->max_depth, sizeof(darray_t *));
88 dimensions[0] = col->w;
89 dimensions[1] = col->h;
90 dimensions[2] = col->k;
91 luminance = darrayCreate(3, dimensions);
96 /* generate sparse multigrid */
97 last_kk = kk = col->k;
99 printf("-> prepare sparse grid");
100 for (d = 0, ww=col->w, hh=col->h; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
101 /* create temporary arrays
102 * copy fist array to level 0 or shrink previous level
105 printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
108 dimensions[2] = last_kk;
110 col->marks[d] = darrayCreate(3, dimensions);
111 if (!col->marks[d]) {
117 dimensions[2] = last_kk-1;
120 col->flows[d] = darrayCreate(4, dimensions);
121 if (!col->flows[d]) {
126 col->flows_i[d] = darrayCreate(4, dimensions);
127 if (!col->flows_i[d]) {
131 if (d < col->scalexyz) {
132 dimensions[0] = (ww+1)>>1;
133 dimensions[1] = (hh+1)>>1;
134 dimensions[2] = kk-1;
137 col->flows_xy[d] = darrayCreate(4, dimensions);
138 if (!col->flows_xy[d]) {
144 for (z = 0; z < kk; z++) {
145 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
146 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
147 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
149 if (col->flow && col->flow_i) {
150 darrayCheckPr(col->flows[d], darrayGetPr(col->flows[d]), col->w*col->h*(kk-1)*2);
151 darrayCheckPr(col->flow, darrayGetPr(col->flow), col->w*col->h*(kk-1)*2);
152 memcpy(darrayGetPr(col->flows[d]), darrayGetPr(col->flow), col->w*col->h*(kk-1)*2*sizeof(double));
153 darrayCheckPr(col->flows_i[d], darrayGetPr(col->flows_i[d]), col->w*col->h*(kk-1)*2);
154 darrayCheckPr(col->flow_i, darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2);
155 memcpy(darrayGetPr(col->flows_i[d]), darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2*sizeof(double));
157 darrayCheckPr(col->marks[d], darrayGetPr(col->marks[d]), col->w*col->h*kk);
158 darrayCheckPr(col->image_mark, darrayGetPr(col->image_mark), col->w*col->h*kk);
159 memcpy(darrayGetPr(col->marks[d]), darrayGetPr(col->image_mark), col->w*col->h*kk*sizeof(double));
160 /* invert, so 0 = marked */
161 quantize(col->marks[d], 0.5, 0, 1);
163 rc = shrink_xy(luminance, NULL);
166 rc = shrink_xy(col->marks[d-1], col->marks[d]);
169 rc = shrink_xy(col->flows[d-1], col->flows[d]);
172 rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
175 /* shrink length of flow vector on smaller grid.
176 * This is not performed when comparing with
177 * the original algorithm. (it is a bug there)
179 multiply(col->flows[d], 0.5);
180 multiply(col->flows_i[d], 0.5);
181 if (d-1 < col->scalexyz) { /* process if previous d shall be shrinked in z direction */
183 rc = shrink_z(luminance, NULL, col->flows[d]);
186 rc = shrink_z(col->marks[d], NULL, col->flows[d]);
189 rc = shrink_flow_z(col->flows[d], NULL);
192 rc = shrink_flow_z(col->flows_i[d], NULL);
196 /* If any of the restricted pixles is 0 (marked),
197 * the result must be 0 (marked) too.
199 quantize(col->marks[d], 0.9999, 1, 0);
201 if (d < col->scalexyz) {
202 /* create x and y shrinked flow vectors for restriction and prolonging process */
203 rc = shrink_xy(col->flows[d], col->flows_xy[d]);
206 multiply(col->flows_xy[d], 0.5);
209 /* generate values array */
214 col->values[d] = darrayCreate(3, dimensions);
215 if (!col->values[d]) {
220 /* generate array of neighbors */
221 #ifndef TEST_NO_RENDER
223 col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
224 if (!col->nb_list[d]) {
230 weighten(luminance, col->nb_list[d]);
234 if (d < col->scalexyz)
237 darrayDestroy(luminance);
242 dimensions[0] = col->w;
243 dimensions[1] = col->h;
244 dimensions[2] = col->k;
246 col->init = darrayCreate(3, dimensions);
253 dimensions[0] = col->w;
254 dimensions[1] = col->h;
255 dimensions[2] = col->k;
256 luminance = darrayCreate(3, dimensions);
257 for (z = 0; z < col->k; z++) {
258 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
259 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
260 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
261 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
262 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
263 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));
265 draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
268 paste(col->marks[0], col->values[0], 0, 0, 0);
272 for (d = 1; d < col->max_depth; d++) {
273 rc = shrink_xy(luminance, NULL);
276 if (d-1 < col->scalexyz)
277 shrink_z(luminance, NULL, col->flows_xy[d-1]);
278 paste(luminance, col->values[0], 0, 0, 0);
280 draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
283 paste(col->marks[d], col->values[0], 0, 0, 0);
286 darrayDestroy(luminance);
288 /* copy test pyramid image as Y-component */
289 for (z = 0; z < col->k; z++) {
290 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
291 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
292 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));
299 darrayDestroy(luminance);
304 int colorize_solve(struct colorize *col, int quick, char **feat_names)
306 int z, d, c, ww, hh, kk;
310 printf("-> apply component to grid and solve");
311 for (c = 1; c < col->comp; c++) {
313 double *value_ptr, *init_ptr, *mark_ptr;
316 printf(" #%d(%s)", c, feat_names[c-1]); fflush(stdout);
318 * - copy component into value array (level 0), otherwise 0
319 * - use maked colors for init array, otherwise use 0
321 for (z = 0; z < col->k; z++) {
322 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
323 darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
324 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));
326 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) {
327 if (d-1 < col->scalexyz) /* scale kk, if last array was scaled */
329 value_ptr = darrayGetPr(col->values[d]);
331 for (j = 0; j < pixles; j++) {
332 darrayCheckPr(col->values[d], value_ptr, 1);
336 value_ptr = darrayGetPr(col->values[0]);
337 init_ptr = darrayGetPr(col->init);
338 mark_ptr = darrayGetPr(col->image_mark);
339 pixles = col->w*col->h*col->k;
340 for (j = 0; j < pixles; j++) {
341 darrayCheckPr(col->image_mark, mark_ptr, 1);
342 if (*mark_ptr++ > 0.5) {
343 darrayCheckPr(col->init, init_ptr, 1);
344 darrayCheckPr(col->values[0], value_ptr, 1);
345 *init_ptr++ = *value_ptr;
347 darrayCheckPr(col->init, init_ptr, 1);
353 #ifndef TEST_NO_RENDER
354 /* colorize component */
356 rc = solve_mg(col->max_depth, col->inner_iter, 1, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz, 0.0);
359 } else /*if (col->target_residual_change)*/ {
360 rc = solve_mg(col->max_depth, col->inner_iter, col->outer_iter, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz, col->target_residual_change);
365 /* original algorithm: handling of outer iteration */
367 for (i = 0; i < col->outer_iter; i++) {
368 rc = solve_mg(col->max_depth, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz, 0.0);
371 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, 0.0);
374 rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz, 0.0);
377 rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz, 0.0);
386 for (d = 1; d < col->max_depth; d++)
387 paste(col->values[d], col->values[0], 0, 0, 0);
390 /* apply colorized component */
391 for (z = 0; z < col->k; z++) {
392 darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
393 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
394 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));
400 #ifdef TEST_SHRINK_XY
401 for (c = 0; c < col->comp; c++) {
403 int dimensions[3] = {col->w, col->h, col->k};
404 for (z = 0; z < col->k; z++)
405 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));
406 dimensions[0] = (col->w+1)>>1;
407 dimensions[1] = (col->h+1)>>1;
408 temp = darrayCreate(3, dimensions);
409 ASSERT(temp, "no memory");
410 rc = shrink_xy(col->values[0], temp);
411 ASSERT(rc >= 0, "no memory");
412 paste(temp, col->values[0], 0, 0, 0);
414 for (z = 0; z < col->k; z++)
415 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));
420 for (c = 0; c < col->comp; c++) {
422 int dimensions[3] = {col->w, col->h, col->k};
423 for (z = 0; z < col->k; z++)
424 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));
425 dimensions[2] = (col->k+1)>>1;
426 temp = darrayCreate(3, dimensions);
427 ASSERT(temp, "no memory");
428 rc = shrink_z(col->values[0], temp, col->flows[0]);
429 ASSERT(rc >= 0, "no memory");
430 paste(temp, col->values[0], 0, 0, 0);
432 for (z = 0; z < col->k; z++)
433 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));
437 #ifdef TEST_EXPAND_XY
438 for (c = 0; c < col->comp; c++) {
440 int dimensions[3] = {col->w, col->h, col->k};
441 for (z = 0; z < col->k; z++)
442 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));
443 dimensions[0] = (col->w+1)>>1;
444 dimensions[1] = (col->h+1)>>1;
445 temp = darrayCreate(3, dimensions);
446 ASSERT(temp, "no memory");
447 rc = shrink_xy(col->values[0], temp);
448 ASSERT(rc >= 0, "no memory");
449 expand_xy(temp, col->values[0]);
451 for (z = 0; z < col->k; z++)
452 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));
457 for (c = 0; c < col->comp; c++) {
459 int dimensions[3] = {col->w, col->h, col->k};
460 for (z = 0; z < col->k; z++)
461 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));
462 dimensions[2] = (col->k+1)>>1;
463 temp = darrayCreate(3, dimensions);
464 ASSERT(temp, "no memory");
465 rc = shrink_z(col->values[0], temp, col->flows[0]);
466 ASSERT(rc >= 0, "no memory");
467 rc = expand_z(temp, col->values[0], col->flows[0]);
468 ASSERT(rc >= 0, "no memory");
470 for (z = 0; z < col->k; z++)
471 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));
481 void colorize_free(struct colorize *col)
486 for (d = 0; d < col->max_depth; d++) {
488 free(col->nb_list[d]);
490 darrayDestroy(col->marks[d]);
492 darrayDestroy(col->values[d]);
494 darrayDestroy(col->flows_xy[d]);
496 darrayDestroy(col->flows[d]);
498 darrayDestroy(col->flows_i[d]);
506 darrayDestroy(col->init);
509 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, double target_residual_change, int quick, int scalexyz, char **feat_names, struct colorize **_col)
511 struct colorize *col;
515 if (!_col || !(*_col)) {
516 col = calloc(sizeof(*col), 1);
517 memset(col, 0, sizeof(*col));
519 col->image_mark = image_mark;
521 col->flow_i = flow_i;
522 col->inner_iter = inner_iter;
523 col->outer_iter = outer_iter;
524 col->target_residual_change = target_residual_change;
525 col->scalexyz = scalexyz;
528 /* use or assign *_col */
536 rc = colorize_prepare(col, quick, alloc);
540 rc = colorize_solve(col, quick, feat_names);