scripts/startup/nodeitems_builtins.py

@@ -335,6 +335,7 @@ compositor_node_categories = [
         NodeItem("CompositorNodeSunBeams"),
         NodeItem("CompositorNodeDenoise"),
         NodeItem("CompositorNodeAntiAliasing"),
+        NodeItem("CompositorNodeKuwahara"),
     ]),
     CompositorNodeCategory("CMP_OP_VECTOR", "Vector", items=[
         NodeItem("CompositorNodeNormal"),

source/blender/blenkernel/BKE_node.h

@@ -1061,6 +1061,7 @@ void BKE_nodetree_remove_layer_n(struct bNodeTree *ntree, struct Scene *scene, i
 #define CMP_NODE_INPAINT 272
 #define CMP_NODE_DESPECKLE 273
 #define CMP_NODE_ANTIALIASING 274
+#define CMP_NODE_KUWAHARA 275
 
 #define CMP_NODE_GLARE 301
 #define CMP_NODE_TONEMAP 302

source/blender/compositor/CMakeLists.txt

@@ -322,6 +322,8 @@ if(WITH_COMPOSITOR_CPU)
     nodes/COM_FilterNode.h
     nodes/COM_InpaintNode.cc
     nodes/COM_InpaintNode.h
+    nodes/COM_KuwaharaNode.h
+    nodes/COM_KuwaharaNode.cc
     nodes/COM_PosterizeNode.cc
     nodes/COM_PosterizeNode.h
 
@@ -349,6 +351,10 @@ if(WITH_COMPOSITOR_CPU)
     operations/COM_GaussianXBlurOperation.h
     operations/COM_GaussianYBlurOperation.cc
     operations/COM_GaussianYBlurOperation.h
+    operations/COM_KuwaharaClassicOperation.h
+    operations/COM_KuwaharaClassicOperation.cc
+    operations/COM_KuwaharaAnisotropicOperation.h
+    operations/COM_KuwaharaAnisotropicOperation.cc
     operations/COM_MovieClipAttributeOperation.cc
     operations/COM_MovieClipAttributeOperation.h
     operations/COM_MovieDistortionOperation.cc

source/blender/compositor/intern/COM_Converter.cc

@@ -62,6 +62,7 @@
 #include "COM_InvertNode.h"
 #include "COM_KeyingNode.h"
 #include "COM_KeyingScreenNode.h"
+#include "COM_KuwaharaNode.h"
 #include "COM_LensDistortionNode.h"
 #include "COM_LuminanceMatteNode.h"
 #include "COM_MapRangeNode.h"
@@ -436,6 +437,9 @@ Node *COM_convert_bnode(bNode *b_node)
     case CMP_NODE_COMBINE_XYZ:
       node = new CombineXYZNode(b_node);
       break;
+    case CMP_NODE_KUWAHARA:
+      node = new KuwaharaNode(b_node);
+      break;
   }
   return node;
 }

source/blender/compositor/nodes/COM_KuwaharaNode.cc

@@ -0,0 +1,103 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "COM_KuwaharaNode.h"
+
+#include "COM_FastGaussianBlurOperation.h"
+#include "COM_KuwaharaAnisotropicOperation.h"
+#include "COM_KuwaharaClassicOperation.h"
+#include "COM_MathBaseOperation.h"
+#include "COM_SetValueOperation.h"
+
+namespace blender::compositor {
+
+void KuwaharaNode::convert_to_operations(NodeConverter &converter,
+                                         const CompositorContext & /*context*/) const
+{
+  const bNode *node = this->get_bnode();
+  const NodeKuwaharaData *data = (const NodeKuwaharaData *)node->storage;
+
+  switch (data->variation) {
+    case CMP_NODE_KUWAHARA_CLASSIC: {
+      KuwaharaClassicOperation *operation = new KuwaharaClassicOperation();
+      operation->set_kernel_size(data->size);
+
+      converter.add_operation(operation);
+      converter.map_input_socket(get_input_socket(0), operation->get_input_socket(0));
+      converter.map_output_socket(get_output_socket(0), operation->get_output_socket());
+      break;
+    }
+
+    case CMP_NODE_KUWAHARA_ANISOTROPIC: {
+      /* Edge detection */
+      auto const_fact = new SetValueOperation();
+      const_fact->set_value(1.0f);
+      converter.add_operation(const_fact);
+
+      auto sobel_x = new ConvolutionFilterOperation();
+      sobel_x->set3x3Filter(1, 0, -1, 2, 0, -2, 1, 0, -1);
+      converter.add_operation(sobel_x);
+      converter.map_input_socket(get_input_socket(0), sobel_x->get_input_socket(0));
+      converter.add_link(const_fact->get_output_socket(0), sobel_x->get_input_socket(1));
+
+      auto sobel_y = new ConvolutionFilterOperation();
+      sobel_y->set3x3Filter(1, 2, 1, 0, 0, 0, -1, -2, -1);
+      converter.add_operation(sobel_y);
+      converter.map_input_socket(get_input_socket(0), sobel_y->get_input_socket(0));
+      converter.add_link(const_fact->get_output_socket(0), sobel_y->get_input_socket(1));
+
+      /* Compute intensity of edges */
+      auto sobel_xx = new MathMultiplyOperation();
+      auto sobel_yy = new MathMultiplyOperation();
+      auto sobel_xy = new MathMultiplyOperation();
+      converter.add_operation(sobel_xx);
+      converter.add_operation(sobel_yy);
+      converter.add_operation(sobel_xy);
+
+      converter.add_link(sobel_x->get_output_socket(0), sobel_xx->get_input_socket(0));
+      converter.add_link(sobel_x->get_output_socket(0), sobel_xx->get_input_socket(1));
+
+      converter.add_link(sobel_y->get_output_socket(0), sobel_yy->get_input_socket(0));
+      converter.add_link(sobel_y->get_output_socket(0), sobel_yy->get_input_socket(1));
+
+      converter.add_link(sobel_x->get_output_socket(0), sobel_xy->get_input_socket(0));
+      converter.add_link(sobel_y->get_output_socket(0), sobel_xy->get_input_socket(1));
+
+      /* Blurring for more robustness. */
+      const int sigma = data->smoothing;
+
+      auto blur_sobel_xx = new FastGaussianBlurOperation();
+      auto blur_sobel_yy = new FastGaussianBlurOperation();
+      auto blur_sobel_xy = new FastGaussianBlurOperation();
+
+      blur_sobel_yy->set_size(sigma, sigma);
+      blur_sobel_xx->set_size(sigma, sigma);
+      blur_sobel_xy->set_size(sigma, sigma);
+
+      converter.add_operation(blur_sobel_xx);
+      converter.add_operation(blur_sobel_yy);
+      converter.add_operation(blur_sobel_xy);
+
+      converter.add_link(sobel_xx->get_output_socket(0), blur_sobel_xx->get_input_socket(0));
+      converter.add_link(sobel_yy->get_output_socket(0), blur_sobel_yy->get_input_socket(0));
+      converter.add_link(sobel_xy->get_output_socket(0), blur_sobel_xy->get_input_socket(0));
+
+      /* Apply anisotropic Kuwahara filter. */
+      KuwaharaAnisotropicOperation *aniso = new KuwaharaAnisotropicOperation();
+      aniso->set_kernel_size(data->size + 4);
+      converter.map_input_socket(get_input_socket(0), aniso->get_input_socket(0));
+      converter.add_operation(aniso);
+
+      converter.add_link(blur_sobel_xx->get_output_socket(0), aniso->get_input_socket(1));
+      converter.add_link(blur_sobel_yy->get_output_socket(0), aniso->get_input_socket(2));
+      converter.add_link(blur_sobel_xy->get_output_socket(0), aniso->get_input_socket(3));
+
+      converter.map_output_socket(get_output_socket(0), aniso->get_output_socket(0));
+
+      break;
+    }
+  }
+}
+
+}  // namespace blender::compositor

source/blender/compositor/nodes/COM_KuwaharaNode.h

@@ -0,0 +1,23 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "COM_Node.h"
+
+namespace blender::compositor {
+
+/**
+ * \brief KuwaharaNode
+ * \ingroup Node
+ */
+
+class KuwaharaNode : public Node {
+ public:
+  KuwaharaNode(bNode *editor_node) : Node(editor_node) {}
+  void convert_to_operations(NodeConverter &converter,
+                             const CompositorContext &context) const override;
+};
+
+}  // namespace blender::compositor

source/blender/compositor/operations/COM_FastGaussianBlurOperation.cc

@@ -11,6 +11,7 @@ namespace blender::compositor {
 FastGaussianBlurOperation::FastGaussianBlurOperation() : BlurBaseOperation(DataType::Color)
 {
   iirgaus_ = nullptr;
+  data_.filtertype = R_FILTER_FAST_GAUSS;
 }
 
 void FastGaussianBlurOperation::execute_pixel(float output[4], int x, int y, void *data)
@@ -68,6 +69,15 @@ void FastGaussianBlurOperation::deinit_execution()
   BlurBaseOperation::deinit_mutex();
 }
 
+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 (i.e. data_). We could use factory pattern to solve this problem. */
+  data_.sizex = size_x;
+  data_.sizey = size_y;
+  sizeavailable_ = true;
+}
+
 void *FastGaussianBlurOperation::initialize_tile_data(rcti *rect)
 {
   lock_mutex();

source/blender/compositor/operations/COM_FastGaussianBlurOperation.h

@@ -28,6 +28,8 @@ class FastGaussianBlurOperation : public BlurBaseOperation {
   void deinit_execution() override;
   void init_execution() override;
 
+  void set_size(int size_x, int size_y);
+
   void get_area_of_interest(int input_idx, const rcti &output_area, rcti &r_input_area) override;
   void update_memory_buffer_started(MemoryBuffer *output,
                                     const rcti &area,

source/blender/compositor/operations/COM_KuwaharaAnisotropicOperation.cc

@@ -0,0 +1,295 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "COM_KuwaharaAnisotropicOperation.h"
+
+#include "BLI_math_base.hh"
+#include "BLI_vector.hh"
+#include "IMB_colormanagement.h"
+
+namespace blender::compositor {
+
+KuwaharaAnisotropicOperation::KuwaharaAnisotropicOperation()
+{
+  this->add_input_socket(DataType::Color);
+  this->add_input_socket(DataType::Color);
+  this->add_input_socket(DataType::Color);
+  this->add_input_socket(DataType::Color);
+
+  this->add_output_socket(DataType::Color);
+
+  this->n_div_ = 8;
+  this->set_kernel_size(5);
+
+  this->flags_.is_fullframe_operation = true;
+}
+
+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()
+{
+  image_reader_ = nullptr;
+}
+
+void KuwaharaAnisotropicOperation::execute_pixel_sampled(float output[4],
+                                                         float x,
+                                                         float y,
+                                                         PixelSampler sampler)
+{
+  const int width = this->get_width();
+  const int height = this->get_height();
+
+  BLI_assert(width == s_xx_reader_->get_width());
+  BLI_assert(height == s_xx_reader_->get_height());
+  BLI_assert(width == s_yy_reader_->get_width());
+  BLI_assert(height == s_yy_reader_->get_height());
+  BLI_assert(width == s_xy_reader_->get_width());
+  BLI_assert(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);
+
+  Vector<double> mean(n_div_);
+  Vector<double> sum(n_div_);
+  Vector<double> var(n_div_);
+  Vector<double> weight(n_div_);
+
+  for (int ch = 0; ch < 3; ch++) {
+    mean.fill(0.0);
+    sum.fill(0.0);
+    var.fill(0.0);
+    weight.fill(0.0);
+
+    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));
+
+        /* Clamp image to avoid artefacts at borders. */
+        const int xx = math::clamp(int(x) + dx2, 0, width - 1);
+        const int yy = math::clamp(int(y) + dy2, 0, height - 1);
+
+        const double ddx2 = double(dx2);
+        const double ddy2 = double(dy2);
+        const double theta = atan2(ddy2, ddx2) + M_PI;
+        const int t = 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++) {
+      double weight_inv = 1.0 / weight[i];
+      mean[i] = weight[i] != 0 ? mean[i] * weight_inv : 0.0;
+      sum[i] = weight[i] != 0 ? sum[i] * weight_inv : 0.0;
+      var[i] = weight[i] != 0 ? var[i] * weight_inv : 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;
+    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)
+{
+  /* 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()
+{
+  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)
+{
+  /* Implementation based on Kyprianidis, Jan & Kang, Henry & Döllner, Jürgen. (2009).
+   * "Image and Video Abstraction by Anisotropic Kuwahara Filtering".
+   * Comput. Graph. Forum. 28. 1955-1963. 10.1111/j.1467-8659.2009.01574.x.
+   * Used reference implementation from lime image processing library (MIT license). */
+
+  MemoryBuffer *image = inputs[0];
+  MemoryBuffer *s_xx = inputs[1];
+  MemoryBuffer *s_yy = inputs[2];
+  MemoryBuffer *s_xy = inputs[3];
+
+  const int width = image->get_width();
+  const int height = image->get_height();
+
+  BLI_assert(width == s_xx->get_width());
+  BLI_assert(height == s_xx->get_height());
+  BLI_assert(width == s_yy->get_width());
+  BLI_assert(height == s_yy->get_height());
+  BLI_assert(width == s_xy->get_width());
+  BLI_assert(height == s_xy->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 (BuffersIterator<float> it = output->iterate_with(inputs, area); !it.is_end(); ++it) {
+    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. */
+    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);
+
+    /* 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);
+
+    Vector<double> mean(n_div_);
+    Vector<double> sum(n_div_);
+    Vector<double> var(n_div_);
+    Vector<double> weight(n_div_);
+
+    for (int ch = 0; ch < 3; ch++) {
+      mean.fill(0.0);
+      sum.fill(0.0);
+      var.fill(0.0);
+      weight.fill(0.0);
+
+      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));
+
+          /* Clamp image to avoid artefacts at borders. */
+          const int xx = math::clamp(x + dx2, 0, width - 1);
+          const int yy = math::clamp(y + dy2, 0, height - 1);
+
+          const double ddx2 = double(dx2);
+          const double ddy2 = double(dy2);
+          const double theta = atan2(ddy2, ddx2) + M_PI;
+          const int t = int(floor(theta / angle)) % n_div_;
+          double d2 = dx2 * dx2 + dy2 * dy2;
+          double g = exp(-d2 / (2.0 * kernel_size_));
+          const 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. */
+      double de = 0.0;
+      double nu = 0.0;
+      for (int i = 0; i < n_div_; i++) {
+        double weight_inv = 1.0 / weight[i];
+        mean[i] = weight[i] != 0 ? mean[i] * weight_inv : 0.0;
+        sum[i] = weight[i] != 0 ? sum[i] * weight_inv : 0.0;
+        var[i] = weight[i] != 0 ? var[i] * weight_inv : 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;
+      it.out[ch] = val;
+    }
+
+    /* No changes for alpha channel. */
+    it.out[3] = image->get_value(x, y, 3);
+  }
+}
+
+}  // namespace blender::compositor

source/blender/compositor/operations/COM_KuwaharaAnisotropicOperation.h

@@ -0,0 +1,36 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "COM_MultiThreadedOperation.h"
+
+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();
+
+  void init_execution() override;
+  void deinit_execution() override;
+  void execute_pixel_sampled(float output[4], float x, float y, PixelSampler sampler) override;
+
+  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,
+                                    Span<MemoryBuffer *> inputs) override;
+};
+
+}  // namespace blender::compositor

source/blender/compositor/operations/COM_KuwaharaClassicOperation.cc

@@ -0,0 +1,199 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "COM_KuwaharaClassicOperation.h"
+
+#include "IMB_colormanagement.h"
+
+namespace blender::compositor {
+
+KuwaharaClassicOperation::KuwaharaClassicOperation()
+{
+  this->add_input_socket(DataType::Color);
+  this->add_output_socket(DataType::Color);
+  this->set_kernel_size(4);
+
+  this->flags_.is_fullframe_operation = true;
+}
+
+void KuwaharaClassicOperation::init_execution()
+{
+  image_reader_ = this->get_input_socket_reader(0);
+}
+
+void KuwaharaClassicOperation::deinit_execution()
+{
+  image_reader_ = nullptr;
+}
+
+void KuwaharaClassicOperation::execute_pixel_sampled(float output[4],
+                                                     float x,
+                                                     float y,
+                                                     PixelSampler sampler)
+{
+  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)
+{
+  kernel_size_ = kernel_size;
+}
+
+int KuwaharaClassicOperation::get_kernel_size()
+{
+  return kernel_size_;
+}
+
+void KuwaharaClassicOperation::update_memory_buffer_partial(MemoryBuffer *output,
+                                                            const rcti &area,
+                                                            Span<MemoryBuffer *> inputs)
+{
+  MemoryBuffer *image = inputs[0];
+
+  for (BuffersIterator<float> it = output->iterate_with(inputs, area); !it.is_end(); ++it) {
+    const int x = it.x;
+    const int y = it.y;
+    it.out[3] = image->get_value(x, y, 3);
+
+    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 < image->get_width() && yy < image->get_height()) {
+            const float v = image->get_value(xx, yy, ch);
+            float color[4];
+            image->read_elem(xx, yy, color);
+            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->get_value(x, y, ch) = mean[min_index];
+    }
+  }
+}
+
+}  // namespace blender::compositor

source/blender/compositor/operations/COM_KuwaharaClassicOperation.h

@@ -0,0 +1,31 @@
+/* SPDX-FileCopyrightText: 2023 Blender Foundation
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "COM_MultiThreadedOperation.h"
+
+namespace blender::compositor {
+
+class KuwaharaClassicOperation : public MultiThreadedOperation {
+  SocketReader *image_reader_;
+
+  int kernel_size_;
+
+ public:
+  KuwaharaClassicOperation();
+
+  void init_execution() override;
+  void deinit_execution() override;
+  void execute_pixel_sampled(float output[4], float x, float y, PixelSampler sampler) override;
+
+  void set_kernel_size(int kernel_size);
+  int get_kernel_size();
+
+  void update_memory_buffer_partial(MemoryBuffer *output,
+                                    const rcti &area,
+                                    Span<MemoryBuffer *> inputs) override;
+};
+
+}  // namespace blender::compositor

source/blender/makesdna/DNA_node_types.h

@@ -881,6 +881,12 @@ typedef struct NodeBilateralBlurData {
   char _pad[2];
 } NodeBilateralBlurData;
 
+typedef struct NodeKuwaharaData {
+  short size;
+  short variation;
+  int smoothing;
+} NodeKuwaharaData;
+
 typedef struct NodeAntiAliasingData {
   float threshold;
   float contrast_limit;
@@ -2136,6 +2142,12 @@ typedef enum CMPNodeGlareType {
   CMP_NODE_GLARE_GHOST = 3,
 } CMPNodeGlareType;
 
+/* Kuwahara Node. Stored in variation */
+typedef enum CMPNodeKuwahara {
+  CMP_NODE_KUWAHARA_CLASSIC = 0,
+  CMP_NODE_KUWAHARA_ANISOTROPIC = 1,
+} CMPNodeKuwahara;
+
 /* Stabilize 2D node. Stored in custom1. */
 typedef enum CMPNodeStabilizeInterpolation {
   CMP_NODE_STABILIZE_INTERPOLATION_NEAREST = 0,

source/blender/makesrna/intern/rna_nodetree.c

@@ -9374,6 +9374,43 @@ static void def_cmp_denoise(StructRNA *srna)
   RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");
 }
 
+static void def_cmp_kuwahara(StructRNA *srna)
+{
+  PropertyRNA *prop;
+
+  RNA_def_struct_sdna_from(srna, "NodeKuwaharaData", "storage");
+
+  static const EnumPropertyItem variation_items[] = {
+      {0, "CLASSIC", 0, "Classic", "Fast but less accurate variation"},
+      {1, "ANISOTROPIC", 0, "Anisotropic", "Accurate but slower variation"},
+      {0, NULL, 0, NULL, NULL},
+  };
+
+  prop = RNA_def_property(srna, "size", PROP_INT, PROP_NONE);
+  RNA_def_property_int_sdna(prop, NULL, "size");
+  RNA_def_property_range(prop, 1.0, 100.0);
+  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");
+
+  prop = RNA_def_property(srna, "variation", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_sdna(prop, NULL, "variation");
+  RNA_def_property_enum_items(prop, variation_items);
+  RNA_def_property_ui_text(prop, "", "Variation of Kuwahara filter to use");
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");
+
+  prop = RNA_def_property(srna, "smoothing", PROP_INT, PROP_NONE);
+  RNA_def_property_int_sdna(prop, NULL, "smoothing");
+  RNA_def_property_range(prop, 0.0, 50.0);
+  RNA_def_property_ui_range(prop, 0, 50, 1, -1);
+  RNA_def_property_ui_text(prop,
+                           "Smoothing",
+                           "Smoothing degree before applying filter. Higher values remove details "
+                           "and give smoother edges");
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");
+}
+
 static void def_cmp_antialiasing(StructRNA *srna)
 {
   PropertyRNA *prop;

source/blender/nodes/NOD_static_types.h

@@ -224,6 +224,7 @@ DefNode(CompositorNode, CMP_NODE_COMBINE_XYZ,    0,                      "COMBIN
 DefNode(CompositorNode, CMP_NODE_SEPARATE_XYZ,   0,                      "SEPARATE_XYZ",   SeparateXYZ,      "Separate XYZ",      ""              )
 DefNode(CompositorNode, CMP_NODE_SEPARATE_COLOR, def_cmp_combsep_color,  "SEPARATE_COLOR", SeparateColor,    "Separate Color",    ""              )
 DefNode(CompositorNode, CMP_NODE_COMBINE_COLOR,  def_cmp_combsep_color,  "COMBINE_COLOR",  CombineColor,     "Combine Color",     ""              )
+DefNode(CompositorNode, CMP_NODE_KUWAHARA,       def_cmp_kuwahara,       "KUWAHARA",       Kuwahara,         "Kuwahara",          ""              )
 
 DefNode(TextureNode,    TEX_NODE_OUTPUT,         def_tex_output,         "OUTPUT",         Output,           "Output",            ""              )
 DefNode(TextureNode,    TEX_NODE_CHECKER,        0,                      "CHECKER",        Checker,          "Checker",           ""              )

source/blender/nodes/composite/CMakeLists.txt

@@ -77,6 +77,7 @@ set(SRC
   nodes/node_composite_invert.cc
   nodes/node_composite_keying.cc
   nodes/node_composite_keyingscreen.cc
+  nodes/node_composite_kuwahara.cc
   nodes/node_composite_lensdist.cc
   nodes/node_composite_levels.cc
   nodes/node_composite_luma_matte.cc

source/blender/nodes/composite/node_composite_register.cc

@@ -64,6 +64,7 @@ void register_composite_nodes()
   register_node_type_cmp_invert();
   register_node_type_cmp_keying();
   register_node_type_cmp_keyingscreen();
+  register_node_type_cmp_kuwahara();
   register_node_type_cmp_lensdist();
   register_node_type_cmp_luma_matte();
   register_node_type_cmp_map_range();

source/blender/nodes/composite/node_composite_register.hh

@@ -60,6 +60,7 @@ void register_node_type_cmp_inpaint();
 void register_node_type_cmp_invert();
 void register_node_type_cmp_keying();
 void register_node_type_cmp_keyingscreen();
+void register_node_type_cmp_kuwahara();
 void register_node_type_cmp_lensdist();
 void register_node_type_cmp_luma_matte();
 void register_node_type_cmp_map_range();

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

@@ -0,0 +1,87 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2023 Blender Foundation */
+
+/** \file
+ * \ingroup cmpnodes
+ */
+
+#include "COM_node_operation.hh"
+
+/* **************** Kuwahara ******************** */
+
+namespace blender::nodes::node_composite_kuwahara_cc {
+
+NODE_STORAGE_FUNCS(NodeKuwaharaData)
+
+static void cmp_node_kuwahara_declare(NodeDeclarationBuilder &b)
+{
+  b.add_input<decl::Color>(N_("Image"))
+      .default_value({1.0f, 1.0f, 1.0f, 1.0f})
+      .compositor_domain_priority(0);
+  b.add_output<decl::Color>(N_("Image"));
+}
+
+static void node_composit_init_kuwahara(bNodeTree * /*ntree*/, bNode *node)
+{
+  NodeKuwaharaData *data = MEM_cnew<NodeKuwaharaData>(__func__);
+  node->storage = data;
+
+  /* Set defaults. */
+  data->size = 4;
+  data->smoothing = 2;
+}
+
+static void node_composit_buts_kuwahara(uiLayout *layout, bContext * /*C*/, PointerRNA *ptr)
+{
+  uiLayout *col;
+
+  col = uiLayoutColumn(layout, false);
+
+  uiItemR(col, ptr, "variation", 0, nullptr, ICON_NONE);
+  uiItemR(col, ptr, "size", 0, nullptr, ICON_NONE);
+
+  const int variation = RNA_enum_get(ptr, "variation");
+
+  if (variation == CMP_NODE_KUWAHARA_ANISOTROPIC) {
+    uiItemR(col, ptr, "smoothing", 0, nullptr, ICON_NONE);
+  }
+}
+
+using namespace blender::realtime_compositor;
+
+class ConvertKuwaharaOperation : public NodeOperation {
+ public:
+  using NodeOperation::NodeOperation;
+
+  void execute() override
+  {
+    get_input("Image").pass_through(get_result("Image"));
+    context().set_info_message("Viewport compositor setup not fully supported");
+  }
+};
+
+static NodeOperation *get_compositor_operation(Context &context, DNode node)
+{
+  return new ConvertKuwaharaOperation(context, node);
+}
+
+}  // namespace blender::nodes::node_composite_kuwahara_cc
+
+void register_node_type_cmp_kuwahara()
+{
+  namespace file_ns = blender::nodes::node_composite_kuwahara_cc;
+
+  static bNodeType ntype;
+
+  cmp_node_type_base(&ntype, CMP_NODE_KUWAHARA, "Kuwahara", NODE_CLASS_OP_FILTER);
+  ntype.declare = file_ns::cmp_node_kuwahara_declare;
+  ntype.draw_buttons = file_ns::node_composit_buts_kuwahara;
+  ntype.initfunc = file_ns::node_composit_init_kuwahara;
+  node_type_storage(
+      &ntype, "NodeKuwaharaData", node_free_standard_storage, node_copy_standard_storage);
+  ntype.get_compositor_operation = file_ns::get_compositor_operation;
+  ntype.realtime_compositor_unsupported_message = N_(
+      "Node not supported in the Viewport compositor");
+
+  nodeRegisterType(&ntype);
+}