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 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;
98 printf("-> prepare sparse grid");
99 for (d = 0, ww=col->w, hh=col->h; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
100 /* create temporary arrays
101 * copy fist array to level 0 or shrink previous level
103 printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
106 dimensions[2] = last_kk;
108 col->marks[d] = darrayCreate(3, dimensions);
109 if (!col->marks[d]) {
115 dimensions[2] = last_kk-1;
118 col->flows[d] = darrayCreate(4, dimensions);
119 if (!col->flows[d]) {
124 col->flows_i[d] = darrayCreate(4, dimensions);
125 if (!col->flows_i[d]) {
129 if (d < col->scalexyz) {
130 dimensions[0] = (ww+1)>>1;
131 dimensions[1] = (hh+1)>>1;
132 dimensions[2] = kk-1;
135 col->flows_xy[d] = darrayCreate(4, dimensions);
136 if (!col->flows_xy[d]) {
142 for (z = 0; z < kk; z++) {
143 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
144 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
145 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
147 if (col->flow && col->flow_i) {
148 darrayCheckPr(col->flows[d], darrayGetPr(col->flows[d]), col->w*col->h*(kk-1)*2);
149 darrayCheckPr(col->flow, darrayGetPr(col->flow), col->w*col->h*(kk-1)*2);
150 memcpy(darrayGetPr(col->flows[d]), darrayGetPr(col->flow), col->w*col->h*(kk-1)*2*sizeof(double));
151 darrayCheckPr(col->flows_i[d], darrayGetPr(col->flows_i[d]), col->w*col->h*(kk-1)*2);
152 darrayCheckPr(col->flow_i, darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2);
153 memcpy(darrayGetPr(col->flows_i[d]), darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2*sizeof(double));
155 darrayCheckPr(col->marks[d], darrayGetPr(col->marks[d]), col->w*col->h*kk);
156 darrayCheckPr(col->image_mark, darrayGetPr(col->image_mark), col->w*col->h*kk);
157 memcpy(darrayGetPr(col->marks[d]), darrayGetPr(col->image_mark), col->w*col->h*kk*sizeof(double));
158 /* invert, so 0 = marked */
159 quantize(col->marks[d], 0.5, 0, 1);
161 rc = shrink_xy(luminance, NULL);
164 rc = shrink_xy(col->marks[d-1], col->marks[d]);
167 rc = shrink_xy(col->flows[d-1], col->flows[d]);
170 rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
173 /* shrink length of flow vector on smaller grid.
174 * This is not performed when comparing with
175 * the original algorithm. (it is a bug there)
177 multiply(col->flows[d], 0.5);
178 multiply(col->flows_i[d], 0.5);
179 if (d-1 < col->scalexyz) { /* process if previous d shall be shrinked in z direction */
181 rc = shrink_z(luminance, NULL, col->flows[d]);
184 rc = shrink_z(col->marks[d], NULL, col->flows[d]);
187 rc = shrink_flow_z(col->flows[d], NULL);
190 rc = shrink_flow_z(col->flows_i[d], NULL);
194 /* If any of the restricted pixles is 0 (marked),
195 * the result must be 0 (marked) too.
197 quantize(col->marks[d], 0.9999, 1, 0);
199 if (d < col->scalexyz) {
200 /* create x and y shrinked flow vectors for restriction and prolonging process */
201 rc = shrink_xy(col->flows[d], col->flows_xy[d]);
204 multiply(col->flows_xy[d], 0.5);
207 /* generate values array */
212 col->values[d] = darrayCreate(3, dimensions);
213 if (!col->values[d]) {
218 /* generate array of neighbors */
219 #ifndef TEST_NO_RENDER
221 col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
222 if (!col->nb_list[d]) {
228 weighten(luminance, col->nb_list[d]);
232 if (d < col->scalexyz)
235 darrayDestroy(luminance);
239 dimensions[0] = col->w;
240 dimensions[1] = col->h;
241 dimensions[2] = col->k;
243 col->init = darrayCreate(3, dimensions);
250 dimensions[0] = col->w;
251 dimensions[1] = col->h;
252 dimensions[2] = col->k;
253 luminance = darrayCreate(3, dimensions);
254 for (z = 0; z < col->k; z++) {
255 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
256 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
257 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
258 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + 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(col->values[0]) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
262 draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
265 paste(col->marks[0], col->values[0], 0, 0, 0);
269 for (d = 1; d < col->max_depth; d++) {
270 rc = shrink_xy(luminance, NULL);
273 if (d-1 < col->scalexyz)
274 shrink_z(luminance, NULL, col->flows_xy[d-1]);
275 paste(luminance, col->values[0], 0, 0, 0);
277 draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
280 paste(col->marks[d], col->values[0], 0, 0, 0);
283 darrayDestroy(luminance);
285 /* copy test pyramid image as Y-component */
286 for (z = 0; z < col->k; z++) {
287 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
288 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
289 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));
296 darrayDestroy(luminance);
301 int colorize_solve(struct colorize *col, int quick, char **feat_names)
303 int z, d, c, ww, hh, kk;
307 printf("-> apply component to grid and solve");
308 for (c = 1; c < col->comp; c++) {
310 double *value_ptr, *init_ptr, *mark_ptr;
312 printf(" #%d(%s)", c, feat_names[c-1]); fflush(stdout);
314 * - copy component into value array (level 0), otherwise 0
315 * - use maked colors for init array, otherwise use 0
317 for (z = 0; z < col->k; z++) {
318 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
319 darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
320 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));
322 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) {
323 if (d-1 < col->scalexyz) /* scale kk, if last array was scaled */
325 value_ptr = darrayGetPr(col->values[d]);
327 for (j = 0; j < pixles; j++) {
328 darrayCheckPr(col->values[d], value_ptr, 1);
332 value_ptr = darrayGetPr(col->values[0]);
333 init_ptr = darrayGetPr(col->init);
334 mark_ptr = darrayGetPr(col->image_mark);
335 pixles = col->w*col->h*col->k;
336 for (j = 0; j < pixles; j++) {
337 darrayCheckPr(col->image_mark, mark_ptr, 1);
338 if (*mark_ptr++ > 0.5) {
339 darrayCheckPr(col->init, init_ptr, 1);
340 darrayCheckPr(col->values[0], value_ptr, 1);
341 *init_ptr++ = *value_ptr;
343 darrayCheckPr(col->init, init_ptr, 1);
349 #ifndef TEST_NO_RENDER
350 /* colorize component */
352 rc = solve_mg(col->max_depth, col->inner_iter, 1, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
357 for (i = 0; i < col->outer_iter; i++) {
358 rc = solve_mg(col->max_depth, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
361 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);
364 rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
367 rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
375 for (d = 1; d < col->max_depth; d++)
376 paste(col->values[d], col->values[0], 0, 0, 0);
379 /* apply colorized component */
380 for (z = 0; z < col->k; z++) {
381 darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
382 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
383 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));
389 #ifdef TEST_SHRINK_XY
390 for (c = 0; c < col->comp; c++) {
392 int dimensions[3] = {col->w, col->h, col->k};
393 for (z = 0; z < col->k; z++)
394 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));
395 dimensions[0] = (col->w+1)>>1;
396 dimensions[1] = (col->h+1)>>1;
397 temp = darrayCreate(3, dimensions);
398 ASSERT(temp, "no memory");
399 rc = shrink_xy(col->values[0], temp);
400 ASSERT(rc >= 0, "no memory");
401 paste(temp, col->values[0], 0, 0, 0);
403 for (z = 0; z < col->k; z++)
404 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));
409 for (c = 0; c < col->comp; c++) {
411 int dimensions[3] = {col->w, col->h, col->k};
412 for (z = 0; z < col->k; z++)
413 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));
414 dimensions[2] = (col->k+1)>>1;
415 temp = darrayCreate(3, dimensions);
416 ASSERT(temp, "no memory");
417 rc = shrink_z(col->values[0], temp, col->flows[0]);
418 ASSERT(rc >= 0, "no memory");
419 paste(temp, col->values[0], 0, 0, 0);
421 for (z = 0; z < col->k; z++)
422 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));
426 #ifdef TEST_EXPAND_XY
427 for (c = 0; c < col->comp; c++) {
429 int dimensions[3] = {col->w, col->h, col->k};
430 for (z = 0; z < col->k; z++)
431 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));
432 dimensions[0] = (col->w+1)>>1;
433 dimensions[1] = (col->h+1)>>1;
434 temp = darrayCreate(3, dimensions);
435 ASSERT(temp, "no memory");
436 rc = shrink_xy(col->values[0], temp);
437 ASSERT(rc >= 0, "no memory");
438 expand_xy(temp, col->values[0]);
440 for (z = 0; z < col->k; z++)
441 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));
446 for (c = 0; c < col->comp; c++) {
448 int dimensions[3] = {col->w, col->h, col->k};
449 for (z = 0; z < col->k; z++)
450 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));
451 dimensions[2] = (col->k+1)>>1;
452 temp = darrayCreate(3, dimensions);
453 ASSERT(temp, "no memory");
454 rc = shrink_z(col->values[0], temp, col->flows[0]);
455 ASSERT(rc >= 0, "no memory");
456 rc = expand_z(temp, col->values[0], col->flows[0]);
457 ASSERT(rc >= 0, "no memory");
459 for (z = 0; z < col->k; z++)
460 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));
470 void colorize_free(struct colorize *col)
475 for (d = 0; d < col->max_depth; d++) {
477 free(col->nb_list[d]);
479 darrayDestroy(col->marks[d]);
481 darrayDestroy(col->values[d]);
483 darrayDestroy(col->flows_xy[d]);
485 darrayDestroy(col->flows[d]);
487 darrayDestroy(col->flows_i[d]);
495 darrayDestroy(col->init);
498 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 quick, int scalexyz, char **feat_names, struct colorize **_col)
500 struct colorize *col;
504 if (!_col || !(*_col)) {
505 col = calloc(sizeof(*col), 1);
506 memset(col, 0, sizeof(*col));
508 col->image_mark = image_mark;
510 col->flow_i = flow_i;
511 col->inner_iter = inner_iter;
512 col->outer_iter = outer_iter;
513 col->scalexyz = scalexyz;
516 /* use or assign *_col */
524 rc = colorize_prepare(col, alloc);
528 rc = colorize_solve(col, quick, feat_names);