source/blender/compositor/nodes/COM_KuwaharaNode.cc

@@ -82,7 +82,7 @@ void KuwaharaNode::convert_to_operations(NodeConverter &converter,
 
       /* Apply anisotropic Kuwahara filter */
       KuwaharaAnisotropicOperation *aniso = new KuwaharaAnisotropicOperation();
-      aniso->set_kernel_size(data->size);
+      aniso->set_kernel_size(data->size + 4);
       converter.map_input_socket(get_input_socket(0), aniso->get_input_socket(0));
       converter.add_operation(aniso);
 

source/blender/compositor/nodes/COM_KuwaharaNode.h

@@ -1,5 +1,5 @@
 /* SPDX-License-Identifier: GPL-2.0-or-later
- * Copyright 2011 Blender Foundation. */
+ * Copyright 2023 Blender Foundation. */
 
 #pragma once
 

source/blender/compositor/operations/COM_FastGaussianBlurOperation.cc

@@ -70,8 +70,8 @@ void FastGaussianBlurOperation::deinit_execution()
 
 void FastGaussianBlurOperation::set_size(int size_x, int size_y)
 {
-  // todo: there should be a better way to use the operation without knowing specifics of the blur
-  // node (data_) Could use factory pattern to solve this problem.
+  /* TODO: there should be a better way to use the operation without knowing specifics of the blur
+   * node (i.e. data_). We could use factory pattern to solve this problem. */
   data_.sizex = size_x;
   data_.sizey = size_y;
   sizeavailable_ = true;

source/blender/compositor/operations/COM_KuwaharaAnisotropicOperation.cc

@@ -17,14 +17,10 @@ KuwaharaAnisotropicOperation::KuwaharaAnisotropicOperation()
   this->add_input_socket(DataType::Color);
   this->add_input_socket(DataType::Color);
 
-  this->add_output_socket(DataType::Color);
-  // output for debug
-  // todo: remove
-  this->add_output_socket(DataType::Color);
-  this->add_output_socket(DataType::Color);
   this->add_output_socket(DataType::Color);
 
-  this->set_kernel_size(8);
+  this->n_div_ = 8;
+  this->set_kernel_size(5);
 
   this->flags_.is_fullframe_operation = true;
 }
@@ -32,6 +28,9 @@ KuwaharaAnisotropicOperation::KuwaharaAnisotropicOperation()
 void KuwaharaAnisotropicOperation::init_execution()
 {
   image_reader_ = this->get_input_socket_reader(0);
+  s_xx_reader_ = this->get_input_socket_reader(1);
+  s_yy_reader_ = this->get_input_socket_reader(2);
+  s_xy_reader_ = this->get_input_socket_reader(3);
 }
 
 void KuwaharaAnisotropicOperation::deinit_execution()
@@ -44,12 +43,124 @@ void KuwaharaAnisotropicOperation::execute_pixel_sampled(float output[4],
                                                          float y,
                                                          PixelSampler sampler)
 {
-  /* Not implemented */
+  BLI_assert(this->get_width() == s_xx_reader_->get_width());
+  BLI_assert(this->get_height() == s_xx_reader_->get_height());
+  BLI_assert(this->get_width() == s_yy_reader_->get_width());
+  BLI_assert(this->get_height() == s_yy_reader_->get_height());
+  BLI_assert(this->get_width() == s_xy_reader_->get_width());
+  BLI_assert(this->get_height() == s_xy_reader_->get_height());
+
+  /* Values recommended by authors in original paper. */
+  const float angle = 2.0 * M_PI / n_div_;
+  const float q = 3.0;
+  const float EPS = 1.0e-10;
+
+  /* For now use green channel to compute orientation */
+  /* todo: convert to HSV and compute orientation and strength on luminance channel */
+  float tmp[4];
+  s_xx_reader_->read(tmp, x, y, nullptr);
+  const float a = tmp[1];
+  s_xy_reader_->read(tmp, x, y, nullptr);
+  const float b = tmp[1];
+  s_yy_reader_->read(tmp, x, y, nullptr);
+  const float c = tmp[1];
+
+  /* Compute egenvalues of structure tensor */
+  const double tr = a + c;
+  const double discr = sqrt((a - b) * (a - b) + 4 * b * c);
+  const double lambda1 = (tr + discr) / 2;
+  const double lambda2 = (tr - discr) / 2;
+
+  /* Compute orientation and its strength based on structure tensor */
+  const double orientation = 0.5 * atan2(2 * b, a - c);
+  const double strength = (lambda1 == 0 && lambda2 == 0) ?
+                              0 :
+                              (lambda1 - lambda2) / (lambda1 + lambda2);
+
+  for (int ch = 0; ch < 3; ch++) {
+    Vector<double> mean(n_div_, 0.0f);
+    Vector<double> sum(n_div_, 0.0f);
+    Vector<double> var(n_div_, 0.0f);
+    Vector<double> weight(n_div_, 0.0f);
+
+    double sx = 1.0f / (strength + 1.0f);
+    double sy = (1.0f + strength) / 1.0f;
+    double theta = -orientation;
+
+    for (int dy = -kernel_size_; dy <= kernel_size_; dy++) {
+      for (int dx = -kernel_size_; dx <= kernel_size_; dx++) {
+        if (dx == 0 && dy == 0)
+          continue;
+
+        /* Rotate and scale the kernel. This is the "anisotropic" part. */
+        int dx2 = (int)(sx * (cos(theta) * dx - sin(theta) * dy));
+        int dy2 = (int)(sy * (sin(theta) * dx + cos(theta) * dy));
+        int xx = x + dx2;
+        int yy = y + dy2;
+
+        /* Clamp image to avoid artefacts at borders. */
+        if (xx < 0) {
+          xx = 0;
+        }
+        if (xx >= this->get_width()) {
+          xx = this->get_width() - 1;
+        }
+        if (yy < 0) {
+          yy = 0;
+        }
+        if (yy >= this->get_height()) {
+          yy = this->get_height() - 1;
+        }
+
+        double ddx2 = (double)dx2;
+        double ddy2 = (double)dy2;
+        double theta = atan2(ddy2, ddx2) + M_PI;
+        int t = static_cast<int>(floor(theta / angle)) % n_div_;
+        double d2 = dx2 * dx2 + dy2 * dy2;
+        double g = exp(-d2 / (2.0 * kernel_size_));
+        float color[4];
+        image_reader_->read(color, xx, yy, nullptr);
+        const double v = color[ch];
+        /* todo(zazizizou): only compute lum once per region */
+        const float lum = IMB_colormanagement_get_luminance(color);
+        /* todo(zazizizou): only compute mean for the selected region */
+        mean[t] += g * v;
+        sum[t] += g * lum;
+        var[t] += g * lum * lum;
+        weight[t] += g;
+      }
+    }
+
+    /* Calculate weighted average */
+    double de = 0.0;
+    double nu = 0.0;
+    for (int i = 0; i < n_div_; i++) {
+      mean[i] = weight[i] != 0 ? mean[i] / weight[i] : 0.0;
+      sum[i] = weight[i] != 0 ? sum[i] / weight[i] : 0.0;
+      var[i] = weight[i] != 0 ? var[i] / weight[i] : 0.0;
+      var[i] = var[i] - sum[i] * sum[i];
+      var[i] = var[i] > FLT_EPSILON ? sqrt(var[i]) : FLT_EPSILON;
+      double w = powf(var[i], -q);
+
+      de += mean[i] * w;
+      nu += w;
+    }
+
+    double val = nu > EPS ? de / nu : 0.0;
+    CLAMP_MAX(val, 1.0f);
+    output[ch] = val;
+  }
+
+  /* No changes for alpha channel */
+  image_reader_->read_sampled(tmp, x, y, sampler);
+  output[3] = tmp[3];
 }
 
 void KuwaharaAnisotropicOperation::set_kernel_size(int kernel_size)
 {
-  kernel_size_ = kernel_size;
+  /* Filter will be split into n_div.
+   * Add n_div / 2 to avoid artefacts such as random black pixels in image. */
+  kernel_size_ = kernel_size + n_div_ / 2;
 }
 
 int KuwaharaAnisotropicOperation::get_kernel_size()
@@ -57,6 +168,11 @@ int KuwaharaAnisotropicOperation::get_kernel_size()
   return kernel_size_;
 }
 
+int KuwaharaAnisotropicOperation::get_n_div()
+{
+  return n_div_;
+}
+
 void KuwaharaAnisotropicOperation::update_memory_buffer_partial(MemoryBuffer *output,
                                                                 const rcti &area,
                                                                 Span<MemoryBuffer *> inputs)
@@ -80,8 +196,8 @@ void KuwaharaAnisotropicOperation::update_memory_buffer_partial(MemoryBuffer *ou
   BLI_assert(image->get_width() == s_xy->get_width());
   BLI_assert(image->get_height() == s_xy->get_height());
 
-  const int n_div = 8;  // recommended by authors in original paper
-  const float angle = 2.0 * M_PI / n_div;
+  /* Values recommended by authors in original paper. */
+  const float angle = 2.0 * M_PI / n_div_;
   const float q = 3.0;
   const float EPS = 1.0e-10;
 
@@ -89,8 +205,8 @@ void KuwaharaAnisotropicOperation::update_memory_buffer_partial(MemoryBuffer *ou
     const int x = it.x;
     const int y = it.y;
 
-    // for now use green channel to compute orientation
-    // todo: convert to HSV and compute orientation and strength on luminance channel
+    /* For now use green channel to compute orientation */
+    /* todo: convert to HSV and compute orientation and strength on luminance channel */
     const float a = s_xx->get_value(x, y, 1);
     const float b = s_xy->get_value(x, y, 1);
     const float c = s_yy->get_value(x, y, 1);
@@ -109,63 +225,75 @@ void KuwaharaAnisotropicOperation::update_memory_buffer_partial(MemoryBuffer *ou
 
     for (int ch = 0; ch < 3; ch++) {
 
-      Vector<float> mean(n_div, 0.0f);
-      Vector<float> sum(n_div, 0.0f);
-      Vector<float> var(n_div, 0.0f);
-      Vector<float> weight(n_div, 0.0f);
+      Vector<double> mean(n_div_, 0.0f);
+      Vector<double> sum(n_div_, 0.0f);
+      Vector<double> var(n_div_, 0.0f);
+      Vector<double> weight(n_div_, 0.0f);
 
-      float sx = 1.0f / (strength + 1.0f);
-      float sy = (1.0f + strength) / 1.0f;
-      float theta = -orientation;
+      double sx = 1.0f / (strength + 1.0f);
+      double sy = (1.0f + strength) / 1.0f;
+      double theta = -orientation;
 
       for (int dy = -kernel_size_; dy <= kernel_size_; dy++) {
         for (int dx = -kernel_size_; dx <= kernel_size_; dx++) {
           if (dx == 0 && dy == 0)
             continue;
 
-          // rotate and scale the kernel. This is the "anisotropic" part.
-          int dx2 = static_cast<int>(sx * (cos(theta) * dx - sin(theta) * dy));
-          int dy2 = static_cast<int>(sy * (sin(theta) * dx + cos(theta) * dy));
+          /* Rotate and scale the kernel. This is the "anisotropic" part. */
+          int dx2 = (int)(sx * (cos(theta) * dx - sin(theta) * dy));
+          int dy2 = (int)(sy * (sin(theta) * dx + cos(theta) * dy));
           int xx = x + dx2;
           int yy = y + dy2;
 
-          if (xx >= 0 && yy >= 0 && xx < image->get_width() && yy < image->get_height()) {
-            float ddx2 = (float)dx2;
-            float ddy2 = (float)dy2;
-            float theta = atan2(ddy2, ddx2) + M_PI;
-            int t = static_cast<int>(floor(theta / angle)) % n_div;
-            float d2 = dx2 * dx2 + dy2 * dy2;
-            float g = exp(-d2 / (2.0 * kernel_size_));
-            float v = image->get_value(xx, yy, ch);
-            float color[4];
-            image->read_elem(xx, yy, color);
-            // todo(zazizizou): only compute lum once per region
-            const float lum = IMB_colormanagement_get_luminance(color);
-            // todo(zazizizou): only compute mean for the selected region
-            mean[t] += g * v;
-            sum[t] += g * lum;
-            var[t] += g * lum * lum;
-            weight[t] += g;
+          /* Clamp image to avoid artefacts at borders. */
+          if (xx < 0) {
+            xx = 0;
           }
+          if (xx >= image->get_width()) {
+            xx = image->get_width() - 1;
+          }
+          if (yy < 0) {
+            yy = 0;
+          }
+          if (yy >= image->get_height()) {
+            yy = image->get_height() - 1;
+          }
+
+          double ddx2 = (double)dx2;
+          double ddy2 = (double)dy2;
+          double theta = atan2(ddy2, ddx2) + M_PI;
+          int t = static_cast<int>(floor(theta / angle)) % n_div_;
+          double d2 = dx2 * dx2 + dy2 * dy2;
+          double g = exp(-d2 / (2.0 * kernel_size_));
+          double v = image->get_value(xx, yy, ch);
+          float color[4];
+          image->read_elem(xx, yy, color);
+          /* todo(zazizizou): only compute lum once per region */
+          const float lum = IMB_colormanagement_get_luminance(color);
+          /* todo(zazizizou): only compute mean for the selected region */
+          mean[t] += g * v;
+          sum[t] += g * lum;
+          var[t] += g * lum * lum;
+          weight[t] += g;
         }
       }
 
-      // Calculate weighted average
-      float de = 0.0;
-      float nu = 0.0;
-      for (int i = 0; i < n_div; i++) {
+      /* Calculate weighted average */
+      double de = 0.0;
+      double nu = 0.0;
+      for (int i = 0; i < n_div_; i++) {
         mean[i] = weight[i] != 0 ? mean[i] / weight[i] : 0.0;
         sum[i] = weight[i] != 0 ? sum[i] / weight[i] : 0.0;
         var[i] = weight[i] != 0 ? var[i] / weight[i] : 0.0;
         var[i] = var[i] - sum[i] * sum[i];
         var[i] = var[i] > FLT_EPSILON ? sqrt(var[i]) : FLT_EPSILON;
-        float w = powf(var[i], -q);
+        double w = powf(var[i], -q);
 
         de += mean[i] * w;
         nu += w;
       }
 
-      float val = nu > EPS ? de / nu : 0.0;
+      double val = nu > EPS ? de / nu : 0.0;
       CLAMP_MAX(val, 1.0f);
       it.out[ch] = val;
     }

source/blender/compositor/operations/COM_KuwaharaAnisotropicOperation.h

@@ -9,8 +9,12 @@ namespace blender::compositor {
 
 class KuwaharaAnisotropicOperation : public MultiThreadedOperation {
   SocketReader *image_reader_;
+  SocketReader *s_xx_reader_;
+  SocketReader *s_yy_reader_;
+  SocketReader *s_xy_reader_;
 
   int kernel_size_;
+  int n_div_;
 
  public:
   KuwaharaAnisotropicOperation();
@@ -21,6 +25,7 @@ class KuwaharaAnisotropicOperation : public MultiThreadedOperation {
 
   void set_kernel_size(int kernel_size);
   int get_kernel_size();
+  int get_n_div();
 
   void update_memory_buffer_partial(MemoryBuffer *output,
                                     const rcti &area,

source/blender/compositor/operations/COM_KuwaharaClassicOperation.cc

@@ -33,7 +33,75 @@ void KuwaharaClassicOperation::execute_pixel_sampled(float output[4],
                                                      float y,
                                                      PixelSampler sampler)
 {
-  /* Not implemented */
+  for (int ch = 0; ch < 3; ch++) {
+    float sum[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+    float mean[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+    float var[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+    int cnt[4] = {0, 0, 0, 0};
+
+    /* Split surroundings of pixel into 4 overlapping regions */
+    for (int dy = -kernel_size_; dy <= kernel_size_; dy++) {
+      for (int dx = -kernel_size_; dx <= kernel_size_; dx++) {
+
+        int xx = x + dx;
+        int yy = y + dy;
+        if (xx >= 0 && yy >= 0 && xx < this->get_width() && yy < this->get_height()) {
+          float color[4];
+          image_reader_->read_sampled(color, xx, yy, sampler);
+          const float v = color[ch];
+          const float lum = IMB_colormanagement_get_luminance(color);
+
+          if (dx <= 0 && dy <= 0) {
+            mean[0] += v;
+            sum[0] += lum;
+            var[0] += lum * lum;
+            cnt[0]++;
+          }
+
+          if (dx >= 0 && dy <= 0) {
+            mean[1] += v;
+            sum[1] += lum;
+            var[1] += lum * lum;
+            cnt[1]++;
+          }
+
+          if (dx <= 0 && dy >= 0) {
+            mean[2] += v;
+            sum[2] += lum;
+            var[2] += lum * lum;
+            cnt[2]++;
+          }
+
+          if (dx >= 0 && dy >= 0) {
+            mean[3] += v;
+            sum[3] += lum;
+            var[3] += lum * lum;
+            cnt[3]++;
+          }
+        }
+      }
+    }
+
+    /* Compute region variances */
+    for (int i = 0; i < 4; i++) {
+      mean[i] = cnt[i] != 0 ? mean[i] / cnt[i] : 0.0f;
+      sum[i] = cnt[i] != 0 ? sum[i] / cnt[i] : 0.0f;
+      var[i] = cnt[i] != 0 ? var[i] / cnt[i] : 0.0f;
+      const float temp = sum[i] * sum[i];
+      var[i] = var[i] > temp ? sqrt(var[i] - temp) : 0.0f;
+    }
+
+    /* Choose the region with lowest variance */
+    float min_var = FLT_MAX;
+    int min_index = 0;
+    for (int i = 0; i < 4; i++) {
+      if (var[i] < min_var) {
+        min_var = var[i];
+        min_index = i;
+      }
+    }
+    output[ch] = mean[min_index];
+  }
 }
 
 void KuwaharaClassicOperation::set_kernel_size(int kernel_size)

source/blender/makesrna/intern/rna_nodetree.c

@@ -9256,7 +9256,7 @@ static void def_cmp_kuwahara(StructRNA *srna)
 
   prop = RNA_def_property(srna, "size", PROP_INT, PROP_NONE);
   RNA_def_property_int_sdna(prop, NULL, "size");
-  RNA_def_property_ui_range(prop, 4, 50, 1, -1);
+  RNA_def_property_ui_range(prop, 1, 100, 1, -1);
   RNA_def_property_ui_text(
       prop, "Size", "Size of filter. Larger values give stronger stylized effect");
   RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");

source/blender/nodes/composite/nodes/node_composite_kuwahara.cc

@@ -6,7 +6,6 @@
  */
 
 #include "COM_node_operation.hh"
-//#include "node_composite_util.hh"
 
 /* **************** Kuwahara ******************** */