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)
  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         col->nb_list = calloc(col->max_depth, sizeof(struct neighbors *));
  50         if (!col->nb_list)
  51                 goto error;
  52         col->marks = calloc(col->max_depth, sizeof(darray_t *));
  53         if (!col->marks)
  54                 goto error;
  55         col->values = calloc(col->max_depth, sizeof(darray_t *));
  56         if (!col->values)
  57                 goto error;
  58         col->flows = calloc(col->max_depth, sizeof(darray_t *));
  59         if (!col->flows)
  60                 goto error;
  61         col->flows_i = calloc(col->max_depth, sizeof(darray_t *));
  62         if (!col->flows_i)
  63                 goto error;
  64         if (col->scalexyz) {
  65                 /* the flows_xy are scaled in xy in the current level */
  66                 col->flows_xy = calloc(col->max_depth, sizeof(darray_t *));
  67                 if (!col->flows_xy)
  68                         goto error;
  69         }
  70         dimensions[0] = col->w;
  71         dimensions[1] = col->h;
  72         dimensions[2] = col->k;
  73         luminance = darrayCreate(3, dimensions);
  74         if (!luminance)
  75                 goto error;
  76         /* generate sparse multigrid */
  77         last_kk = kk = col->k;
  78         printf("-> prepare sparse grid");
  79         for (d = 0, ww=col->w, hh=col->h; d < col->max_depth; d++, ww=(ww+1)>>1, hh=(hh+1)>>1) {
  80                 /* create temporary arrays
  81                  * copy fist array to level 0 or shrink previous level
  82                  */
  83                 printf(" #%d(w=%d h=%d k=%d)", d, ww, hh, kk); fflush(stdout);
  84                 dimensions[0] = ww;
  85                 dimensions[1] = hh;
  86                 dimensions[2] = last_kk;
  87                 col->marks[d] = darrayCreate(3, dimensions);
  88                 if (!col->marks[d])
  89                         goto error;
  90                 dimensions[0] = ww;
  91                 dimensions[1] = hh;
  92                 dimensions[2] = last_kk-1;
  93                 dimensions[3] = 2;
  94                 col->flows[d] = darrayCreate(4, dimensions);
  95                 if (!col->flows[d])
  96                         goto error;
  97                 col->flows_i[d] = darrayCreate(4, dimensions);
  98                 if (!col->flows_i[d])
  99                         goto error;
 100                 if (d < col->scalexyz) {
 101                         dimensions[0] = (ww+1)>>1;
 102                         dimensions[1] = (hh+1)>>1;
 103                         dimensions[2] = kk-1;
 104                         dimensions[3] = 2;
 105                         col->flows_xy[d] = darrayCreate(4, dimensions);
 106                         if (!col->flows_xy[d])
 107                                 goto error;
 108                 }
 109                 if (d == 0) {
 110                         for (z = 0; z < kk; z++) {
 111                                 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
 112                                 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
 113                                 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h*sizeof(double));
 114                         }
 115                         if (col->flow && col->flow_i) {
 116                                 darrayCheckPr(col->flows[d], darrayGetPr(col->flows[d]), col->w*col->h*(kk-1)*2);
 117                                 darrayCheckPr(col->flow, darrayGetPr(col->flow), col->w*col->h*(kk-1)*2);
 118                                 memcpy(darrayGetPr(col->flows[d]), darrayGetPr(col->flow), col->w*col->h*(kk-1)*2*sizeof(double));
 119                                 darrayCheckPr(col->flows_i[d], darrayGetPr(col->flows_i[d]), col->w*col->h*(kk-1)*2);
 120                                 darrayCheckPr(col->flow_i, darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2);
 121                                 memcpy(darrayGetPr(col->flows_i[d]), darrayGetPr(col->flow_i), col->w*col->h*(kk-1)*2*sizeof(double));
 122                         }
 123                         darrayCheckPr(col->marks[d], darrayGetPr(col->marks[d]), col->w*col->h*kk);
 124                         darrayCheckPr(col->image_mark, darrayGetPr(col->image_mark), col->w*col->h*kk);
 125                         memcpy(darrayGetPr(col->marks[d]), darrayGetPr(col->image_mark), col->w*col->h*kk*sizeof(double));
 126                         /* invert, so 0 = marked */
 127                         quantize(col->marks[d], 0.5, 0, 1);
 128                 } else {
 129                         rc = shrink_xy(luminance, NULL);
 130                         if (rc < 0)
 131                                 goto error;
 132                         rc = shrink_xy(col->marks[d-1], col->marks[d]);
 133                         if (rc < 0)
 134                                 goto error;
 135                         rc = shrink_xy(col->flows[d-1], col->flows[d]);
 136                         if (rc < 0)
 137                                 goto error;
 138                         rc = shrink_xy(col->flows_i[d-1], col->flows_i[d]);
 139                         if (rc < 0)
 140                                 goto error;
 141                         /* shrink length of flow vector on smaller grid.
 142                          * This is not performed when comparing with
 143                          * the original algorithm. (it is a bug there)
 144                          */
 145                         multiply(col->flows[d], 0.5);
 146                         multiply(col->flows_i[d], 0.5);
 147                         if (d-1 < col->scalexyz) { /* process if previous d shall be shrinked in z direction */
 148                                 /* shrink time */
 149                                 rc = shrink_z(luminance, NULL, col->flows[d]);
 150                                 if (rc < 0)
 151                                         goto error;
 152                                 rc = shrink_z(col->marks[d], NULL, col->flows[d]);
 153                                 if (rc < 0)
 154                                         goto error;
 155                                 rc = shrink_flow_z(col->flows[d], NULL);
 156                                 if (rc < 0)
 157                                         goto error;
 158                                 rc = shrink_flow_z(col->flows_i[d], NULL);
 159                                 if (rc < 0)
 160                                         goto error;
 161                         }
 162                         /* If any of the restricted pixles is 0 (marked),
 163                          * the result must be 0 (marked) too.
 164                          */
 165                         quantize(col->marks[d], 0.9999, 1, 0);
 166                 }
 167                 if (d < col->scalexyz) {
 168                         /* create x and y shrinked flow vectors for restriction and prolonging process */
 169                         rc = shrink_xy(col->flows[d], col->flows_xy[d]);
 170                         if (rc < 0)
 171                                 goto error;
 172                         multiply(col->flows_xy[d], 0.5);
 173                 }
 174
 175                 /* generate values array */
 176                 dimensions[0] = ww;
 177                 dimensions[1] = hh;
 178                 dimensions[2] = kk;
 179                 col->values[d] = darrayCreate(3, dimensions);
 180                 if (!col->values[d])
 181                         goto error;
 182
 183                 /* generate array of neighbors */
 184 #ifndef TEST_NO_RENDER
 185                 col->nb_list[d] = gen_neighbor(ww, hh, kk, col->flows[d], col->flows_i[d]);
 186                 if (!col->nb_list[d])
 187                         goto error;
 188                 /* weighten */
 189                 weighten(luminance, col->nb_list[d]);
 190 #endif
 191
 192                 last_kk = kk;
 193                 if (d < col->scalexyz)
 194                         kk = (kk+1)>>1;
 195         }
 196         darrayDestroy(luminance);
 197         luminance = NULL;
 198         printf("\n");
 199
 200         dimensions[0] = col->w;
 201         dimensions[1] = col->h;
 202         dimensions[2] = col->k;
 203         col->init = darrayCreate(3, dimensions);
 204         if (!col->init)
 205                 goto error;
 206
 207 #ifdef TEST_PYRAMID
 208         dimensions[0] = col->w;
 209         dimensions[1] = col->h;
 210         dimensions[2] = col->k;
 211         luminance = darrayCreate(3, dimensions);
 212         for (z = 0; z < col->k; z++) {
 213                 darrayCheckPr(luminance, darrayGetPr(luminance) + z*col->w*col->h, col->w*col->h);
 214                 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
 215                 memcpy(darrayGetPr(luminance) + z*col->w*col->h, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h*sizeof(double));
 216                 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
 217                 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*comp, col->w*col->h);
 218                 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));
 219 #ifdef TEST_FLOW
 220                 draw_flow(col->values[0], col->flows[0], col->flows_i[0]);
 221 #endif
 222 #ifdef TEST_MARKED
 223                 paste(col->marks[0], col->values[0], 0, 0, 0);
 224 #endif
 225
 226         }
 227         for (d = 1; d < col->max_depth; d++) {
 228                 rc = shrink_xy(luminance, NULL);
 229                 if (rc < 0)
 230                         goto error;
 231                 if (d-1 < col->scalexyz)
 232                         shrink_z(luminance, NULL, col->flows_xy[d-1]);
 233                 paste(luminance, col->values[0], 0, 0, 0);
 234 #ifdef TEST_FLOW
 235                 draw_flow(col->values[0], col->flows_x[d], col->flows_ix[d]);
 236 #endif
 237 #ifdef TEST_MARKED
 238                 paste(col->marks[d], col->values[0], 0, 0, 0);
 239 #endif
 240         }
 241         darrayDestroy(luminance);
 242         luminance = NULL;
 243         /* copy test pyramid image as Y-component */
 244         for (z = 0; z < col->k; z++) {
 245                 darrayCheckPr(col->image, darrayGetPr(col->image) + z*col->w*col->h*col->comp, col->w*col->h);
 246                 darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
 247                 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));
 248         }
 249 #endif
 250
 251         rc = 0;
 252
 253 error:
 254         darrayDestroy(luminance);
 255
 256         return rc;
 257 }
 258
 259 int colorize_solve(struct colorize *col, int quick)
 260 {
 261         int z, d, c, ww, hh, kk;
 262         int rc = -1;
 263
 264         if (!quick)
 265                 printf("-> apply component to grid and solve");
 266         for (c = 1; c < col->comp; c++) {
 267                 int pixles, j;
 268                 double *value_ptr, *init_ptr, *mark_ptr;
 269
 270                 printf(" #%d", c); fflush(stdout);
 271                 /* apply component
 272                  * - copy component into value array (level 0), otherwise 0
 273                  * - use maked colors for init array, otherwise use 0
 274                  */
 275                 for (z = 0; z < col->k; z++) {
 276                         darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
 277                         darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
 278                         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));
 279                 }
 280                 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) {
 281                         if (d-1 < col->scalexyz) /* scale kk, if last array was scaled */
 282                                 kk=(kk+1)>>1;
 283                         value_ptr = darrayGetPr(col->values[d]);
 284                         pixles = ww*hh*kk;
 285                         for (j = 0; j < pixles; j++) {
 286                                 darrayCheckPr(col->values[d], value_ptr, 1);
 287                                 *value_ptr++ = 0.0;
 288                         }
 289                 }
 290                 value_ptr = darrayGetPr(col->values[0]);
 291                 init_ptr = darrayGetPr(col->init);
 292                 mark_ptr = darrayGetPr(col->image_mark);
 293                 pixles = col->w*col->h*col->k;
 294                 for (j = 0; j < pixles; j++) {
 295                         darrayCheckPr(col->image_mark, mark_ptr, 1);
 296                         if (*mark_ptr++ > 0.5) {
 297                                 darrayCheckPr(col->init, init_ptr, 1);
 298                                 darrayCheckPr(col->values[0], value_ptr, 1);
 299                                 *init_ptr++ = *value_ptr;
 300                         } else {
 301                                 darrayCheckPr(col->init, init_ptr, 1);
 302                                 *init_ptr++ = 0.0;
 303                         }
 304                         value_ptr++;
 305                 }
 306
 307 #ifndef TEST_NO_RENDER
 308                 /* colorize component */
 309                 if (quick) {
 310                         rc = solve_mg(col->max_depth, col->inner_iter, 1, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
 311                         if (rc < 0)
 312                                 goto error;
 313                 } else {
 314                         int i;
 315                         for (i = 0; i < col->outer_iter; i++) {
 316                                 rc = solve_mg(col->max_depth, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
 317                                 if (rc < 0)
 318                                         goto error;
 319                                 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);
 320                                 if (rc < 0)
 321                                         goto error;
 322                                 rc = solve_mg(2, col->inner_iter, 2, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
 323                                 if (rc < 0)
 324                                         goto error;
 325                                 rc = solve_mg(1, col->inner_iter, 4, col->values, col->marks, col->init, col->nb_list, col->flows_xy, col->scalexyz);
 326                                 if (rc < 0)
 327                                         goto error;
 328                         }
 329                 }
 330 #endif
 331
 332 #ifdef TEST_PYRAMID
 333                 for (d = 1; d < col->max_depth; d++)
 334                         paste(col->values[d], col->values[0], 0, 0, 0);
 335 #endif
 336
 337                 /* apply colorized component */
 338                 for (z = 0; z < col->k; z++) {
 339                         darrayCheckPr(col->image, darrayGetPr(col->image) + c*col->w*col->h + z*col->w*col->h*col->comp, col->w*col->h);
 340                         darrayCheckPr(col->values[0], darrayGetPr(col->values[0]) + z*col->w*col->h, col->w*col->h);
 341                         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));
 342                 }
 343         }
 344         if (!quick)
 345                 printf("\n");
 346
 347 #ifdef TEST_SHRINK_XY
 348         for (c = 0; c < col->comp; c++) {
 349                 darray_t *temp;
 350                 int dimensions[3] = {col->w, col->h, col->k};
 351                 for (z = 0; z < col->k; z++)
 352                         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));
 353                 dimensions[0] = (col->w+1)>>1;
 354                 dimensions[1] = (col->h+1)>>1;
 355                 temp = darrayCreate(3, dimensions);
 356                 ASSERT(temp, "no memory");
 357                 rc = shrink_xy(col->values[0], temp);
 358                 ASSERT(rc >= 0, "no memory");
 359                 paste(temp, col->values[0], 0, 0, 0);
 360                 darrayDestroy(temp);
 361                 for (z = 0; z < col->k; z++)
 362                         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));
 363         }
 364 #endif
 365
 366 #ifdef TEST_SHRINK_Z
 367         for (c = 0; c < col->comp; c++) {
 368                 darray_t *temp;
 369                 int dimensions[3] = {col->w, col->h, col->k};
 370                 for (z = 0; z < col->k; z++)
 371                         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));
 372                 dimensions[2] = (col->k+1)>>1;
 373                 temp = darrayCreate(3, dimensions);
 374                 ASSERT(temp, "no memory");
 375                 rc = shrink_z(col->values[0], temp, col->flows[0]);
 376                 ASSERT(rc >= 0, "no memory");
 377                 paste(temp, col->values[0], 0, 0, 0);
 378                 darrayDestroy(temp);
 379                 for (z = 0; z < col->k; z++)
 380                         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));
 381         }
 382 #endif
 383
 384 #ifdef TEST_EXPAND_XY
 385         for (c = 0; c < col->comp; c++) {
 386                 darray_t *temp;
 387                 int dimensions[3] = {col->w, col->h, col->k};
 388                 for (z = 0; z < col->k; z++)
 389                         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));
 390                 dimensions[0] = (col->w+1)>>1;
 391                 dimensions[1] = (col->h+1)>>1;
 392                 temp = darrayCreate(3, dimensions);
 393                 ASSERT(temp, "no memory");
 394                 rc = shrink_xy(col->values[0], temp);
 395                 ASSERT(rc >= 0, "no memory");
 396                 expand_xy(temp, col->values[0]);
 397                 darrayDestroy(temp);
 398                 for (z = 0; z < col->k; z++)
 399                         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         }
 401 #endif
 402
 403 #ifdef TEST_EXPAND_Z
 404         for (c = 0; c < col->comp; c++) {
 405                 darray_t *temp;
 406                 int dimensions[3] = {col->w, col->h, col->k};
 407                 for (z = 0; z < col->k; z++)
 408                         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));
 409                 dimensions[2] = (col->k+1)>>1;
 410                 temp = darrayCreate(3, dimensions);
 411                 ASSERT(temp, "no memory");
 412                 rc = shrink_z(col->values[0], temp, col->flows[0]);
 413                 ASSERT(rc >= 0, "no memory");
 414                 rc = expand_z(temp, col->values[0], col->flows[0]);
 415                 ASSERT(rc >= 0, "no memory");
 416                 darrayDestroy(temp);
 417                 for (z = 0; z < col->k; z++)
 418                         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));
 419         }
 420 #endif
 421
 422         rc = 0;
 423
 424 error:
 425         return rc;
 426 }
 427
 428 void colorize_free(struct colorize *col)
 429 {
 430         int d;
 431
 432         /* free memory */
 433         for (d = 0; d < col->max_depth; d++) {
 434                 if (col->nb_list)
 435                         free(col->nb_list[d]);
 436                 if (col->marks)
 437                         darrayDestroy(col->marks[d]);
 438                 if (col->values)
 439                         darrayDestroy(col->values[d]);
 440                 if (col->flows_xy)
 441                         darrayDestroy(col->flows_xy[d]);
 442                 if (col->flows)
 443                         darrayDestroy(col->flows[d]);
 444                 if (col->flows_i)
 445                         darrayDestroy(col->flows_i[d]);
 446         }
 447         free(col->nb_list);
 448         free(col->marks);
 449         free(col->values);
 450         free(col->flows_xy);
 451         free(col->flows);
 452         free(col->flows_i);
 453         darrayDestroy(col->init);
 454 }
 455
 456 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)
 457 {
 458         struct colorize col;
 459         int rc = -1;
 460
 461         memset(&col, 0, sizeof(col));
 462         col.image = image;
 463         col.image_mark = image_mark;
 464         col.flow = flow;
 465         col.flow_i = flow_i;
 466         col.inner_iter = inner_iter;
 467         col.outer_iter = outer_iter;
 468         col.scalexyz = scalexyz;
 469
 470         rc = colorize_prepare(&col);
 471         if (rc < 0)
 472                 goto error;
 473
 474         rc = colorize_solve(&col, 0);
 475         if (rc < 0)
 476                 goto error;
 477
 478         rc = 0;
 479
 480 error:
 481         colorize_free(&col);
 482
 483         return rc;
 484 }