96bd385260e569589b961cece0e8c784097deda4
[colorize.git] / lib / colorize.c
1 #include <stdio.h>
2 #include <string.h>
3 #include <stdlib.h>
4 #include <math.h>
5 #include "darray.h"
6 #include "multigrid.h"
7 #include "colorize.h"
8
9 /* colorize image array using image_mark
10  *
11  * image:
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
16  *
17  * image_mark:
18  *      first dimansion = width
19  *      second dimension = height
20  *      fourth dimension = time
21  */
22 int colorize_prepare(struct colorize *col, int alloc)
23 {
24         int dimensions[4];
25         int ww, hh, kk, last_kk, d, temp, z;
26         darray_t *luminance = NULL;
27         int rc = -1;
28
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)");
35
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");
40
41         /* Calculate maximum depth
42          * Scale down dimensions until at least 2x2 remains
43          */
44         temp = min(col->w, col->h);
45         for (col->max_depth = 0; temp > 2; col->max_depth++)
46                 temp=(temp+1)>>1;
47         ASSERT(col->max_depth >= 2, "image too small");
48
49         if (alloc)
50                 col->nb_list = calloc(col->max_depth, sizeof(struct neighbors *));
51         if (!col->nb_list) {
52                 rc = -1;
53                 goto error;
54         }
55         if (alloc)
56                 col->marks = calloc(col->max_depth, sizeof(darray_t *));
57         if (!col->marks) {
58                 rc = -1;
59                 goto error;
60         }
61         if (alloc)
62                 col->values = calloc(col->max_depth, sizeof(darray_t *));
63         if (!col->values) {
64                 rc = -1;
65                 goto error;
66         }
67         if (alloc)
68                 col->flows = calloc(col->max_depth, sizeof(darray_t *));
69         if (!col->flows) {
70                 rc = -1;
71                 goto error;
72         }
73         if (alloc)
74                 col->flows_i = calloc(col->max_depth, sizeof(darray_t *));
75         if (!col->flows_i) {
76                 rc = -1;
77                 goto error;
78         }
79         if (col->scalexyz) {
80                 /* the flows_xy are scaled in xy in the current level */
81                 if (alloc)
82                         col->flows_xy = calloc(col->max_depth, sizeof(darray_t *));
83                 if (!col->flows_xy) {
84                         rc = -1;
85                         goto error;
86                 }
87         }
88         dimensions[0] = col->w;
89         dimensions[1] = col->h;
90         dimensions[2] = col->k;
91         luminance = darrayCreate(3, dimensions);
92         if (!luminance) {
93                 rc = -1;
94                 goto error;
95         }
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
102                  */
103                 printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
104                 dimensions[0] = ww;
105                 dimensions[1] = hh;
106                 dimensions[2] = last_kk;
107                 if (alloc)
108                         col->marks[d] = darrayCreate(3, dimensions);
109                 if (!col->marks[d]) {
110                         rc = -1;
111                         goto error;
112                 }
113                 dimensions[0] = ww;
114                 dimensions[1] = hh;
115                 dimensions[2] = last_kk-1;
116                 dimensions[3] = 2;
117                 if (alloc)
118                         col->flows[d] = darrayCreate(4, dimensions);
119                 if (!col->flows[d]) {
120                         rc = -1;
121                         goto error;
122                 }
123                 if (alloc)
124                         col->flows_i[d] = darrayCreate(4, dimensions);
125                 if (!col->flows_i[d]) {
126                         rc = -1;
127                         goto error;
128                 }
129                 if (d < col->scalexyz) {
130                         dimensions[0] = (ww+1)>>1;
131                         dimensions[1] = (hh+1)>>1;
132                         dimensions[2] = kk-1;
133                         dimensions[3] = 2;
134                         if (alloc)
135                                 col->flows_xy[d] = darrayCreate(4, dimensions);
136                         if (!col->flows_xy[d]) {
137                                 rc = -1;
138                                 goto error;
139                         }
140                 }
141                 if (d == 0) {
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));
146                         }
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));
154                         }
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);
160                 } else {
161                         rc = shrink_xy(luminance, NULL);
162                         if (rc < 0)
163                                 goto error;
164                         rc = shrink_xy(col->marks[d-1], col->marks[d]);
165                         if (rc < 0)
166                                 goto error;
167                         rc = shrink_xy(col->flows[d-1], col->flows[d]);
168                         if (rc < 0)
169                                 goto error;
170                         rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
171                         if (rc < 0)
172                                 goto error;
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)
176                          */
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 */
180                                 /* shrink time */
181                                 rc = shrink_z(luminance, NULL, col->flows[d]);
182                                 if (rc < 0)
183                                         goto error;
184                                 rc = shrink_z(col->marks[d], NULL, col->flows[d]);
185                                 if (rc < 0)
186                                         goto error;
187                                 rc = shrink_flow_z(col->flows[d], NULL);
188                                 if (rc < 0)
189                                         goto error;
190                                 rc = shrink_flow_z(col->flows_i[d], NULL);
191                                 if (rc < 0)
192                                         goto error;
193                         }
194                         /* If any of the restricted pixles is 0 (marked),
195                          * the result must be 0 (marked) too.
196                          */
197                         quantize(col->marks[d], 0.9999, 1, 0);
198                 }
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]);
202                         if (rc < 0)
203                                 goto error;
204                         multiply(col->flows_xy[d], 0.5);
205                 }
206
207                 /* generate values array */
208                 dimensions[0] = ww;
209                 dimensions[1] = hh;
210                 dimensions[2] = kk;
211                 if (alloc)
212                         col->values[d] = darrayCreate(3, dimensions);
213                 if (!col->values[d]) {
214                         rc = -1;
215                         goto error;
216                 }
217
218                 /* generate array of neighbors */
219 #ifndef TEST_NO_RENDER
220                 if (alloc)
221                         col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
222                 if (!col->nb_list[d]) {
223                         rc = -1;
224                         goto error;
225                 }
226                 /* weighten */
227                 if (alloc)
228                         weighten(luminance, col->nb_list[d]);
229 #endif
230
231                 last_kk = kk;
232                 if (d < col->scalexyz)
233                         kk = (kk+1)>>1;
234         }
235         darrayDestroy(luminance);
236         luminance = NULL;
237         printf("\n");
238
239         dimensions[0] = col->w;
240         dimensions[1] = col->h;
241         dimensions[2] = col->k;
242         if (alloc)
243                 col->init = darrayCreate(3, dimensions);
244         if (!col->init) {
245                 rc = -1;
246                 goto error;
247         }
248
249 #ifdef TEST_PYRAMID
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));
261 #ifdef TEST_FLOW
262                 draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
263 #endif
264 #ifdef TEST_MARKED
265                 paste(col->marks[0], col->values[0], 0, 0, 0);
266 #endif
267
268         }
269         for (d = 1; d < col->max_depth; d++) {
270                 rc = shrink_xy(luminance, NULL);
271                 if (rc < 0)
272                         goto error;
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);
276 #ifdef TEST_FLOW
277                 draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
278 #endif
279 #ifdef TEST_MARKED
280                 paste(col->marks[d], col->values[0], 0, 0, 0);
281 #endif
282         }
283         darrayDestroy(luminance);
284         luminance = NULL;
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));
290         }
291 #endif
292
293         rc = 0;
294
295 error:
296         darrayDestroy(luminance);
297
298         return rc;
299 }
300
301 int colorize_solve(struct colorize *col, int quick, char **feat_names)
302 {
303         int z, d, c, ww, hh, kk;
304         int rc = -1;
305
306         if (!quick)
307                 printf("-> apply component to grid and solve");
308         for (c = 1; c < col->comp; c++) {
309                 int pixles, j;
310                 double *value_ptr, *init_ptr, *mark_ptr;
311
312                 printf(" #%d(%s)", c, feat_names[c-1]); fflush(stdout);
313                 /* apply component
314                  * - copy component into value array (level 0), otherwise 0
315                  * - use maked colors for init array, otherwise use 0
316                  */
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));
321                 }
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 */
324                                 kk=(kk+1)>>1;   
325                         value_ptr = darrayGetPr(col->values[d]);
326                         pixles = ww*hh*kk;
327                         for (j = 0; j < pixles; j++) {
328                                 darrayCheckPr(col->values[d], value_ptr, 1);
329                                 *value_ptr++ = 0.0;
330                         }
331                 }
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;
342                         } else {
343                                 darrayCheckPr(col->init, init_ptr, 1);
344                                 *init_ptr++ = 0.0;
345                         }
346                         value_ptr++;
347                 }
348
349 #ifndef TEST_NO_RENDER
350                 /* colorize component */
351                 if (quick) {
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);
353                         if (rc < 0)
354                                 goto error;
355                 } else {
356                         int i;
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);
359                                 if (rc < 0)
360                                         goto error;
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);
362                                 if (rc < 0)
363                                         goto error;
364                                 rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
365                                 if (rc < 0)
366                                         goto error;
367                                 rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
368                                 if (rc < 0)
369                                         goto error;
370                         }
371                 }
372 #endif
373
374 #ifdef TEST_PYRAMID
375                 for (d = 1; d < col->max_depth; d++)
376                         paste(col->values[d], col->values[0], 0, 0, 0);
377 #endif
378
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));
384                 }
385         }
386         if (!quick)
387                 printf("\n");
388
389 #ifdef TEST_SHRINK_XY
390         for (c = 0; c < col->comp; c++) {
391                 darray_t *temp;
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);
402                 darrayDestroy(temp);
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));
405         }
406 #endif
407
408 #ifdef TEST_SHRINK_Z
409         for (c = 0; c < col->comp; c++) {
410                 darray_t *temp;
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);
420                 darrayDestroy(temp);
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));
423         }
424 #endif
425
426 #ifdef TEST_EXPAND_XY
427         for (c = 0; c < col->comp; c++) {
428                 darray_t *temp;
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]);
439                 darrayDestroy(temp);
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));
442         }
443 #endif
444
445 #ifdef TEST_EXPAND_Z
446         for (c = 0; c < col->comp; c++) {
447                 darray_t *temp;
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");
458                 darrayDestroy(temp);
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));
461         }
462 #endif
463
464         rc = 0;
465
466 error:
467         return rc;
468 }
469
470 void colorize_free(struct colorize *col)
471 {
472         int d;
473
474         /* free memory */
475         for (d = 0; d < col->max_depth; d++) {
476                 if (col->nb_list)
477                         free(col->nb_list[d]);
478                 if (col->marks)
479                         darrayDestroy(col->marks[d]);
480                 if (col->values)
481                         darrayDestroy(col->values[d]);
482                 if (col->flows_xy)
483                         darrayDestroy(col->flows_xy[d]);
484                 if (col->flows)
485                         darrayDestroy(col->flows[d]);
486                 if (col->flows_i)
487                         darrayDestroy(col->flows_i[d]);
488         }
489         free(col->nb_list);
490         free(col->marks);
491         free(col->values);
492         free(col->flows_xy);
493         free(col->flows);
494         free(col->flows_i);
495         darrayDestroy(col->init);
496 }
497
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)
499 {
500         struct colorize *col;
501         int rc = -1;
502         int alloc = 0;
503
504         if (!_col || !(*_col)) {
505                 col = calloc(sizeof(*col), 1);
506                 memset(col, 0, sizeof(*col));
507                 col->image = image;
508                 col->image_mark = image_mark;
509                 col->flow = flow;
510                 col->flow_i = flow_i;
511                 col->inner_iter = inner_iter;
512                 col->outer_iter = outer_iter;
513                 col->scalexyz = scalexyz;
514                 alloc = 1;
515         }
516         /* use or assign *_col */
517         if (_col) {
518                 if (!(*_col))
519                         *_col = col;
520                 else
521                         col = *_col;
522         }
523
524         rc = colorize_prepare(col, alloc);
525         if (rc < 0)
526                 goto error;
527
528         rc = colorize_solve(col, quick, feat_names);
529         if (rc < 0)
530                 goto error;
531
532         rc = 0;
533
534 error:
535         if (!_col)
536                 colorize_free(col);
537
538         return rc;
539 }