operations/common/convolution-matrix.c

   1 /* This file is an image processing operation for GEGL
   2  *
   3  * GEGL is free software; you can redistribute it and/or
   4  * modify it under the terms of the GNU Lesser General Public
   5  * License as published by the Free Software Foundation; either
   6  * version 3 of the License, or (at your option) any later version.
   7  *
   8  * GEGL is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11  * Lesser General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU Lesser General Public
  14  * License along with GEGL; if not, see <https://www.gnu.org/licenses/>.
  15  *
  16  * Copyright (C) 1997 Lauri Alanko <la@iki.fi>
  17  * Copyright 2011 Robert Sasu (sasu.robert@gmail.com)
  18  */
  19
  20 #include "config.h"
  21 #include <glib/gi18n-lib.h>
  22
  23 #ifdef GEGL_PROPERTIES
  24
  25 property_double (a1, _("(1,1)"), 0.0)
  26 property_double (a2, _("(1,2)"), 0.0)
  27 property_double (a3, _("(1,3)"), 0.0)
  28 property_double (a4, _("(1,4)"), 0.0)
  29 property_double (a5, _("(1,5)"), 0.0)
  30 property_double (b1, _("(2,1)"), 0.0)
  31 property_double (b2, _("(2,2)"), 0.0)
  32 property_double (b3, _("(2,3)"), 0.0)
  33 property_double (b4, _("(2,4)"), 0.0)
  34 property_double (b5, _("(2,5)"), 0.0)
  35 property_double (c1, _("(3,1)"), 0.0)
  36 property_double (c2, _("(3,2)"), 0.0)
  37 property_double (c3, _("(3,3)"), 1.0)
  38 property_double (c4, _("(3,4)"), 0.0)
  39 property_double (c5, _("(3,5)"), 0.0)
  40 property_double (d1, _("(4,1)"), 0.0)
  41 property_double (d2, _("(4,2)"), 0.0)
  42 property_double (d3, _("(4,3)"), 0.0)
  43 property_double (d4, _("(4,4)"), 0.0)
  44 property_double (d5, _("(4,5)"), 0.0)
  45 property_double (e1, _("(5,1)"), 0.0)
  46 property_double (e2, _("(5,2)"), 0.0)
  47 property_double (e3, _("(5,3)"), 0.0)
  48 property_double (e4, _("(5,4)"), 0.0)
  49 property_double (e5, _("(5,5)"), 0.0)
  50
  51 property_double (divisor, _("Divisor"), 1.0)
  52     ui_range    (-1000.0, 1000.0)
  53     ui_meta     ("sensitive", "! normalize")
  54
  55 property_double (offset, _("Offset"), 0.0)
  56     value_range (-1.0, 1.0)
  57     ui_meta     ("sensitive", "! normalize")
  58
  59 property_boolean (red,   _("Red channel"),   TRUE)
  60 property_boolean (green, _("Green channel"), TRUE)
  61 property_boolean (blue,  _("Blue channel"),  TRUE)
  62 property_boolean (alpha, _("Alpha channel"), TRUE)
  63
  64 property_boolean (normalize,    _("Normalize"),       TRUE)
  65 property_boolean (alpha_weight, _("Alpha-weighting"), TRUE)
  66
  67 property_enum (border, _("Border"),
  68                GeglAbyssPolicy, gegl_abyss_policy,
  69                GEGL_ABYSS_CLAMP)
  70
  71 #else
  72
  73 #define GEGL_OP_AREA_FILTER
  74 #define GEGL_OP_NAME     convolution_matrix
  75 #define GEGL_OP_C_SOURCE convolution-matrix.c
  76
  77 #include "gegl-op.h"
  78 #include <stdio.h>
  79
  80 #define MAX_MATRIX_SIZE 5
  81
  82 static gboolean
  83 enough_with_3x3 (GeglProperties *o)
  84 {
  85   if (o->a1 == 0.0 &&
  86       o->a2 == 0.0 &&
  87       o->a3 == 0.0 &&
  88       o->a4 == 0.0 &&
  89       o->a5 == 0.0 &&
  90       o->b1 == 0.0 &&
  91       o->b5 == 0.0 &&
  92       o->c1 == 0.0 &&
  93       o->c5 == 0.0 &&
  94       o->d1 == 0.0 &&
  95       o->d5 == 0.0 &&
  96       o->e1 == 0.0 &&
  97       o->e2 == 0.0 &&
  98       o->e3 == 0.0 &&
  99       o->e4 == 0.0 &&
 100       o->e5 == 0.0)
 101   {
 102     return TRUE;
 103   }
 104   return FALSE;
 105 }
 106
 107 static void
 108 prepare (GeglOperation *operation)
 109 {
 110   const Babl *space = gegl_operation_get_source_space (operation, "input");
 111
 112   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
 113
 114   GeglProperties          *o       = GEGL_PROPERTIES (operation);
 115   if (enough_with_3x3 (o))
 116     op_area->left = op_area->right = op_area->top = op_area->bottom = 1; /* 3 */
 117   else
 118     op_area->left = op_area->right = op_area->top = op_area->bottom = 2; /* 5 */
 119
 120   gegl_operation_set_format (operation, "output",
 121                              babl_format_with_space ("RGBA float", space));
 122 }
 123
 124 static void
 125 make_matrix (GeglProperties  *o,
 126              gfloat           matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 127              gint             matrix_size)
 128 {
 129   if (matrix_size == 3)
 130   {
 131     matrix[0][0] = o->b2;
 132     matrix[0][1] = o->b3;
 133     matrix[0][2] = o->b4;
 134
 135     matrix[1][0] = o->c2;
 136     matrix[1][1] = o->c3;
 137     matrix[1][2] = o->c4;
 138
 139     matrix[2][0] = o->d2;
 140     matrix[2][1] = o->d3;
 141     matrix[2][2] = o->d4;
 142   }
 143   else
 144   {
 145     matrix[0][0] = o->a1;
 146     matrix[0][1] = o->a2;
 147     matrix[0][2] = o->a3;
 148     matrix[0][3] = o->a4;
 149     matrix[0][4] = o->a5;
 150
 151     matrix[1][0] = o->b1;
 152     matrix[1][1] = o->b2;
 153     matrix[1][2] = o->b3;
 154     matrix[1][3] = o->b4;
 155     matrix[1][4] = o->b5;
 156
 157     matrix[2][0] = o->c1;
 158     matrix[2][1] = o->c2;
 159     matrix[2][2] = o->c3;
 160     matrix[2][3] = o->c4;
 161     matrix[2][4] = o->c5;
 162
 163     matrix[3][0] = o->d1;
 164     matrix[3][1] = o->d2;
 165     matrix[3][2] = o->d3;
 166     matrix[3][3] = o->d4;
 167     matrix[3][4] = o->d5;
 168
 169     matrix[4][0] = o->e1;
 170     matrix[4][1] = o->e2;
 171     matrix[4][2] = o->e3;
 172     matrix[4][3] = o->e4;
 173     matrix[4][4] = o->e5;
 174   }
 175 }
 176
 177 static gboolean
 178 normalize_div_off (gfloat  matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 179                    gint    matrix_size,
 180                    gfloat *divisor,
 181                    gfloat *offset)
 182 {
 183   gint      x, y;
 184   gboolean  valid = FALSE;
 185   gfloat    sum   = 0.0;
 186
 187   for (y = 0; y < matrix_size; y++)
 188     for (x = 0; x < matrix_size; x++)
 189       {
 190         sum += matrix[x][y];
 191         if (matrix[x][y] != 0.0)
 192           valid = TRUE;
 193       }
 194
 195   if (sum > 0)
 196     {
 197       *offset  = 0.0;
 198       *divisor = sum;
 199     }
 200   else if (sum < 0)
 201     {
 202       *offset  = 1.0;
 203       *divisor = -sum;
 204     }
 205   else
 206     {
 207       *offset  = 0.5;
 208       *divisor = 1;
 209     }
 210
 211   return valid;
 212 }
 213
 214 static gint inline
 215 matrix_center_offset (const GeglRectangle *extended,
 216                       gint                 matrix_size)
 217 {
 218   return (extended->width + 1) * (matrix_size / 2);
 219 }
 220
 221 static void inline
 222 convolve_pixel_componentwise (GeglProperties       *o,
 223                 gfloat               *src_buf,
 224                 gfloat               *dst_buf,
 225                 const GeglRectangle  *result,
 226                 const GeglRectangle  *extended,
 227                 gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 228                 gint                  matrix_size,
 229                 gint                  d_offset,
 230                 gint                  ss_offset,
 231                 gint                  xx,
 232                 gint                  yy,
 233                 gfloat                matrixsum,
 234                 gfloat                inv_divisor,
 235                 gfloat                offset)
 236 {
 237   gint   i;
 238   gint   s_stride = (extended->width - matrix_size) * 4;
 239   for (i = 0; i < 4; i++)
 240     {
 241       gfloat sum = 0.0;
 242       gint s_offset = ss_offset;
 243
 244       if ((i == 0 && o->red)   ||
 245           (i == 1 && o->green) ||
 246           (i == 2 && o->blue)  ||
 247           (i == 3 && o->alpha))
 248         {
 249           gint x, y;
 250           for (y = 0; y < matrix_size; y++)
 251           {
 252             for (x = 0; x < matrix_size; x++)
 253               {
 254                 sum += matrix[x][y] * src_buf[s_offset + i];
 255                 s_offset += 4;
 256               }
 257             s_offset += s_stride;
 258           }
 259           sum = sum * inv_divisor;
 260           sum += offset;
 261           dst_buf[d_offset + i] = sum;
 262         }
 263       else
 264         {
 265           s_offset += 4 * matrix_center_offset (extended, matrix_size);
 266           dst_buf[d_offset + i] = src_buf[s_offset + i];
 267         }
 268     }
 269 }
 270
 271 static void inline
 272 convolve_pixel (GeglProperties       *o,
 273                 gfloat               *src_buf,
 274                 gfloat               *dst_buf,
 275                 const GeglRectangle  *result,
 276                 const GeglRectangle  *extended,
 277                 gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 278                 gint                  matrix_size,
 279                 gint                  d_offset,
 280                 gint                  ss_offset,
 281                 gint                  xx,
 282                 gint                  yy,
 283                 gfloat                matrixsum,
 284                 gfloat                inv_divisor,
 285                 gfloat                offset)
 286 {
 287   gint   i;
 288   gint   s_stride = (extended->width - matrix_size) * 4;
 289   for (i = 0; i < 4; i++)
 290     {
 291       gfloat sum = 0.0;
 292       gint s_offset = ss_offset;
 293       gint x, y;
 294       for (y = 0; y < matrix_size; y++)
 295       {
 296         for (x = 0; x < matrix_size; x++)
 297           {
 298             sum += matrix[x][y] * src_buf[s_offset + i];
 299             s_offset += 4;
 300           }
 301         s_offset += s_stride;
 302       }
 303       sum = sum * inv_divisor;
 304       sum += offset;
 305       dst_buf[d_offset + i] = sum;
 306     }
 307 }
 308
 309 static void inline
 310 convolve_pixel_alpha_weight (GeglProperties       *o,
 311                              gfloat               *src_buf,
 312                              gfloat               *dst_buf,
 313                              const GeglRectangle  *result,
 314                              const GeglRectangle  *extended,
 315                              gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 316                              gint                  matrix_size,
 317                              gint                  d_offset,
 318                              gint                  ss_offset,
 319                              gint                  xx,
 320                              gint                  yy,
 321                              gfloat                matrixsum,
 322                              gfloat                inv_divisor,
 323                              gfloat                offset)
 324 {
 325   gint   i;
 326   gint   s_stride = (extended->width - matrix_size) * 4;
 327
 328   for (i = 0; i < 3; i++)
 329     {
 330       gfloat sum      = 0.0;
 331       gint   s_offset = ss_offset;
 332       gint x, y;
 333       for (y = 0; y < matrix_size; y++)
 334       {
 335         for (x = 0; x < matrix_size; x++)
 336           {
 337             sum += matrix[x][y] * src_buf[s_offset + i] * src_buf[s_offset + 3];
 338             s_offset += 4;
 339           }
 340         s_offset += s_stride;
 341       }
 342       sum = sum * inv_divisor + offset;
 343       dst_buf[d_offset + i] = sum;
 344     }
 345     {
 346       gfloat sum      = 0.0;
 347       gfloat alphasum = 0.0;
 348       gint   s_offset = ss_offset;
 349       gint x, y;
 350
 351       for (y = 0; y < matrix_size; y++)
 352       {
 353         for (x = 0; x < matrix_size; x++)
 354           {
 355             float val = matrix[x][y] * src_buf[s_offset + i];
 356             sum += val;
 357             alphasum += fabsf (val);
 358             s_offset += 4;
 359           }
 360         s_offset += s_stride;
 361       }
 362
 363       if (alphasum > 0.0)
 364       {
 365         sum = sum * inv_divisor;
 366         sum = sum * matrixsum / alphasum + offset;
 367       }
 368       else
 369         sum = offset;
 370       dst_buf[d_offset + i] = sum;
 371     }
 372 }
 373
 374 static void inline
 375 convolve_pixel_alpha_weight_componentwise (GeglProperties       *o,
 376                              gfloat               *src_buf,
 377                              gfloat               *dst_buf,
 378                              const GeglRectangle  *result,
 379                              const GeglRectangle  *extended,
 380                              gfloat                matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE],
 381                              gint                  matrix_size,
 382                              gint                  d_offset,
 383                              gint                  ss_offset,
 384                              gint                  xx,
 385                              gint                  yy,
 386                              gfloat                matrixsum,
 387                              gfloat                inv_divisor,
 388                              gfloat                offset)
 389 {
 390   gint   i;
 391   gint   s_stride = (extended->width - matrix_size) * 4;
 392
 393   for (i = 0; i < 3; i++)
 394     {
 395       gfloat sum      = 0.0;
 396       gint   s_offset = ss_offset;
 397
 398       if ((i == 0 && o->red)   ||
 399           (i == 1 && o->green) ||
 400           (i == 2 && o->blue))
 401         {
 402           gint x, y;
 403           for (y = 0; y < matrix_size; y++)
 404           {
 405             for (x = 0; x < matrix_size; x++)
 406               {
 407                 sum += matrix[x][y] * src_buf[s_offset + i] * src_buf[s_offset + 3];
 408                 s_offset += 4;
 409               }
 410             s_offset += s_stride;
 411           }
 412           sum = sum * inv_divisor + offset;
 413         }
 414       else
 415         {
 416           s_offset += 4 * matrix_center_offset (extended, matrix_size);
 417           sum = src_buf[s_offset + i];
 418         }
 419       dst_buf[d_offset + i] = sum;
 420     }
 421     {
 422       gfloat sum      = 0.0;
 423       gfloat alphasum = 0.0;
 424       gint   s_offset = ss_offset;
 425
 426       if (o->alpha)
 427         {
 428           gint x, y;
 429
 430           for (y = 0; y < matrix_size; y++)
 431           {
 432             for (x = 0; x < matrix_size; x++)
 433               {
 434                 float val = matrix[x][y] * src_buf[s_offset + i];
 435                 sum += val;
 436                 alphasum += fabsf (val);
 437                 s_offset += 4;
 438               }
 439             s_offset += s_stride;
 440           }
 441
 442
 443           if (alphasum != 0)
 444           {
 445             sum = sum * inv_divisor;
 446             sum = sum * matrixsum / alphasum + offset;
 447           }
 448           else
 449             sum = offset;
 450         }
 451       else
 452         {
 453           s_offset += 4 * matrix_center_offset (extended, matrix_size);
 454           sum = src_buf[s_offset + i];
 455         }
 456       dst_buf[d_offset + i] = sum;
 457     }
 458 }
 459
 460 static gboolean
 461 process (GeglOperation       *operation,
 462          GeglBuffer          *input,
 463          GeglBuffer          *output,
 464          const GeglRectangle *result,
 465          gint                 level)
 466 {
 467   GeglProperties          *o       = GEGL_PROPERTIES (operation);
 468   GeglOperationAreaFilter *op_area = GEGL_OPERATION_AREA_FILTER (operation);
 469   const Babl              *format  = gegl_operation_get_format (operation, "output");
 470   GeglRectangle            rect;
 471   gfloat                  *src_buf;
 472   gfloat                  *dst_buf;
 473   gfloat                   matrix[MAX_MATRIX_SIZE][MAX_MATRIX_SIZE]={{0,}};
 474   gfloat                   matrixsum = 0.0;
 475   gfloat                   divisor = o->divisor;
 476   gfloat                   offset = o->offset;
 477   gfloat                   inv_divisor;
 478   gint                     x, y;
 479   gint                     matrix_size = MAX_MATRIX_SIZE;
 480
 481   if (enough_with_3x3 (o))
 482     matrix_size = 3;
 483
 484   make_matrix (o, matrix, matrix_size);
 485
 486   if (o->normalize)
 487     normalize_div_off (matrix, matrix_size, &divisor, &offset);
 488   inv_divisor = 1.0 / divisor;
 489
 490   for (x = 0; x < matrix_size; x++)
 491     for (y = 0; y < matrix_size; y++)
 492       matrixsum += fabsf (matrix[x][y]);
 493
 494   rect.x      = result->x - op_area->left;
 495   rect.width  = result->width + op_area->left + op_area->right;
 496   rect.y      = result->y - op_area->top;
 497   rect.height = result->height + op_area->top + op_area->bottom;
 498
 499   src_buf = g_new (gfloat, rect.width * rect.height * 4);
 500   dst_buf = g_new (gfloat, result->width * result->height * 4);
 501
 502   gegl_buffer_get (input, &rect, 1.0, format, src_buf,
 503                    GEGL_AUTO_ROWSTRIDE, o->border);
 504
 505   if (o->divisor != 0)
 506     {
 507       gint ss_offset = (result->y - matrix_size/2 - rect.y) * rect.width * 4 +
 508                        (result->x - matrix_size/2 - rect.x) * 4;
 509       int d_offset = 0;
 510       if (o->alpha_weight)
 511         {
 512           if (o->red == FALSE || o->green == FALSE || o->blue == FALSE || o->alpha == FALSE)
 513           {
 514             gint x, y;
 515             for (y = result->y; y < result->height + result->y; y++)
 516             {
 517               for (x = result->x; x < result->width + result->x; x++)
 518               {
 519                 convolve_pixel_alpha_weight_componentwise (o, src_buf, dst_buf, result, &rect,
 520                                            matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor, offset);
 521                 d_offset += 4;
 522                 ss_offset += 4;
 523               }
 524               ss_offset += (rect.width - result->width) * 4;
 525             }
 526           }
 527           else
 528           {
 529             gint x, y;
 530             for (y = result->y; y < result->height + result->y; y++)
 531             {
 532               for (x = result->x; x < result->width + result->x; x++)
 533               {
 534                 convolve_pixel_alpha_weight (o, src_buf, dst_buf, result, &rect,
 535                                            matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor, offset);
 536                 d_offset += 4;
 537                 ss_offset += 4;
 538               }
 539               ss_offset += (rect.width - result->width) * 4;
 540             }
 541           }
 542         }
 543       else
 544         {
 545           if (o->red == FALSE || o->green == FALSE || o->blue == FALSE || o->alpha == FALSE)
 546           {
 547             gint x, y;
 548             for (y = result->y; y < result->height + result->y; y++)
 549             {
 550               for (x = result->x; x < result->width + result->x; x++)
 551               {
 552                 convolve_pixel_componentwise (o, src_buf, dst_buf, result, &rect,
 553                                 matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor, offset);
 554                 d_offset += 4;
 555                 ss_offset += 4;
 556               }
 557               ss_offset += (rect.width - result->width) * 4;
 558             }
 559           }
 560           else
 561           {
 562             gint x, y;
 563             for (y = result->y; y < result->height + result->y; y++)
 564             {
 565               for (x = result->x; x < result->width + result->x; x++)
 566               {
 567                 convolve_pixel (o, src_buf, dst_buf, result, &rect,
 568                                 matrix, matrix_size, d_offset, ss_offset, x, y, matrixsum, inv_divisor, offset);
 569                 d_offset += 4;
 570                 ss_offset += 4;
 571               }
 572               ss_offset += (rect.width - result->width) * 4;
 573             }
 574           }
 575         }
 576         gegl_buffer_set (output, result, 0, format,
 577                          dst_buf, GEGL_AUTO_ROWSTRIDE);
 578     }
 579   else
 580     {
 581       gegl_buffer_set (output, &rect, 0, format,
 582                        src_buf, GEGL_AUTO_ROWSTRIDE);
 583     }
 584
 585   g_free (src_buf);
 586   g_free (dst_buf);
 587
 588   return TRUE;
 589 }
 590
 591 static GeglRectangle
 592 get_bounding_box (GeglOperation *operation)
 593 {
 594   GeglRectangle  result = { 0, };
 595   GeglRectangle *in_rect;
 596
 597   in_rect = gegl_operation_source_get_bounding_box (operation, "input");
 598   if (!in_rect)
 599     return result;
 600
 601   return *in_rect;
 602 }
 603
 604 static GeglAbyssPolicy
 605 get_abyss_policy (GeglOperation *operation,
 606                   const gchar   *input_pad)
 607 {
 608   GeglProperties *o = GEGL_PROPERTIES (operation);
 609
 610   return o->border;
 611 }
 612
 613 static void
 614 gegl_op_class_init (GeglOpClass *klass)
 615 {
 616   GeglOperationClass           *operation_class;
 617   GeglOperationFilterClass     *filter_class;
 618   GeglOperationAreaFilterClass *area_class;
 619
 620   operation_class = GEGL_OPERATION_CLASS (klass);
 621   filter_class    = GEGL_OPERATION_FILTER_CLASS (klass);
 622   area_class      = GEGL_OPERATION_AREA_FILTER_CLASS (klass);
 623
 624   area_class->get_abyss_policy             = get_abyss_policy;
 625   filter_class->process                    = process;
 626   operation_class->prepare                 = prepare;
 627   operation_class->get_bounding_box        = get_bounding_box;
 628
 629   gegl_operation_class_set_keys (operation_class,
 630     "categories",  "generic",
 631     "name",        "gegl:convolution-matrix",
 632     "reference-hash", "22d2d47a2da3d3e7cd402ea9fa1a3a25",
 633     "reference-hashB", "4eddc0aaa970a59ee8a813627874cdf3",
 634     "title",       _("Convolution Matrix"),
 635     "description", _("Apply a generic 5x5 convolution matrix"),
 636     NULL);
 637 }
 638
 639 #endif