Laplace operator in grids

2025-07-29 21:11 UTC o4-mini-high Open in ChatGPT ↗

Think of your 2D Unity grid as a big chessboard of scalar values (like heights, temperatures, whatever you like). The Laplacian at any given cell is just a measure of how different that cell’s value is from the values of its four immediate neighbors (up, down, left, right). Here’s how to see it step by step:

Label the grid around your cell
```
    (x, y+1)
       ↑
(x-1, y) • (x+1, y)
       ↓
    (x, y-1)
```
And let the center cell be at coordinates (x, y) with value C, and its neighbors have values:
- L = grid[x-1][y] (left)
- R = grid[x+1][y] (right)
- U = grid[x][y+1] (up)
- D = grid[x][y-1] (down)
Compute the “neighbor sum minus center”
The simplest finite‐difference Laplacian is:
```
float laplace = L + R + U + D - 4 * C;
```
What you’re doing is comparing the center cell C against the average of its four neighbors. If laplace is:
- Positive: the neighbors on average are higher than the center (so the center sits in a “valley”).
- Negative: the neighbors on average are lower than the center (so the center is on a “peak”).
- Zero: the center is exactly at the same “level” as its neighbors (flat).
Why this works
- In the continuous world, the Laplacian ∆f is “second derivative in x plus second derivative in y.”
- On a grid, second derivatives become “difference of differences,” which reduces to “sum of neighbor values minus 4× center.”
- This tells you the local curvature of your field: where it wants to smooth out (diffuse) or sharpen.

Putting it in a Unity shader or script

float ComputeLaplacian(float[,] grid, int x, int y) {
    float C = grid[x, y];
    float L = grid[x-1, y];
    float R = grid[x+1, y];
    float U = grid[x, y+1];
    float D = grid[x, y-1];
    return (L + R + U + D) - 4f * C;
}

You’d loop over all (x,y) (skipping the border or using wrap/mirror) and store or use that laplace value to drive diffusion, curvature flow, sharpening filters, etc.

Intuition in Graphics Terms

Heat diffusion: if your grid is temperature, the Laplacian tells you where heat should flow (from hot peaks into cooler valleys).
Blur/Sharpen filters: image processing uses almost the same kernel—sum of neighbors minus center—so you can blur (diffuse) or detect edges (highlight where the Laplacian is large).
Terrain smoothing: apply the Laplacian repeatedly and you’ll iron out spikes and fill holes.

So no scary ∇·∇, just a tiny 5‑point stencil on your grid: add those four neighbors, subtract four times the center, and you’ve got your Laplace!

To simulate diffusion on your 2D grid, you’ll integrate the diffusion equation in time:

$ \frac{\partial u}{\partial t} = D,\Delta u $

where

$u(x,y,t)$ is your scalar field (e.g. temperature),
$D$ is the diffusion coefficient,
$\Delta u$ is the Laplacian you’ve computed.

On a discrete grid with time step $\Delta t$, a simple explicit Euler update is:

u_new[x,y] = u_old[x,y] + D * Δt * laplacian(x,y);

Step‑by‑step

Choose your parameters

float D   = 0.1f;   // diffusion rate: how “fast” it spreads
float dt  = 0.01f;  // time step: must be small enough for stability
int   W   = width;  // grid width
int   H   = height; // grid height

Compute Laplacians
First pass: for every interior cell (x,y), compute

float C = u[x, y];
float L = u[x-1, y];
float R = u[x+1, y];
float U = u[x, y+1];
float D = u[x, y-1];
lap[x,y] = (L + R + U + D) - 4f * C;

Store these into a temporary 2D array lap[,].

Update your field
Second pass:

for (int x = 1; x < W-1; x++) {
    for (int y = 1; y < H-1; y++) {
        u_new[x,y] = u[x,y] + diffusionCoefficient * dt * lap[x,y];
    }
}

Swap u and u_new (or copy back).

Enforce boundaries
- Dirichlet: clamp edges to some fixed value.
- Neumann (zero‑flux): mirror the neighbour, e.g. u[0,y] = u[1,y].
- Periodic: wrap around, u[0,y] = u[W-2,y], etc.

Complete Example

public class Diffusion2D {
    int    W, H;
    float[,] u, uNew, lap;
    public float diffusionRate = 0.1f;
    public float dt            = 0.01f;

    public Diffusion2D(int width, int height) {
        W = width; H = height;
        u    = new float[W, H];
        uNew = new float[W, H];
        lap  = new float[W, H];
    }

    // Call once per frame (or fixed update)
    public void Step() {
        // 1) compute laplacian stencil
        for (int x = 1; x < W-1; x++) {
            for (int y = 1; y < H-1; y++) {
                float C = u[x, y];
                lap[x,y] = (u[x-1,y] + u[x+1,y] + u[x,y-1] + u[x,y+1])
                            - 4f * C;
            }
        }

        // 2) update u using Euler forward
        for (int x = 1; x < W-1; x++) {
            for (int y = 1; y < H-1; y++) {
                uNew[x,y] = u[x,y] + diffusionRate * dt * lap[x,y];
            }
        }

        // 3) handle boundaries (zero‑flux example)
        for (int x = 0; x < W; x++) {
            uNew[x,0]    = uNew[x,1];
            uNew[x,H-1]  = uNew[x,H-2];
        }
        for (int y = 0; y < H; y++) {
            uNew[0,y]    = uNew[1,y];
            uNew[W-1,y]  = uNew[W-2,y];
        }

        // 4) swap arrays
        var tmp = u;
        u    = uNew;
        uNew = tmp;
    }
}

Stability Tip

For the explicit scheme, you need $ D,\frac{\Delta t}{(\Delta x)^2} \le \tfrac12 $ (where $\Delta x$ is your grid spacing—often 1 unit). If you exceed this, you’ll see “blow‑up” (unstable oscillations). Either shrink dt, lower D, or switch to an implicit solver (which is unconditionally stable but more involved).

That’s it! Each step you’re using your Laplacian stencil to push your cell values toward the average of their neighbors, which manifests as diffusion over time.

At its heart, diffusion is just “stuff flowing downhill” from high‑concentration regions into low‑concentration ones until everything even‑s out. Adding the Laplacian to your field implements exactly that “flow downhill” rule on your grid.

1. The Laplacian as local curvature

On your grid, at each cell $(x,y)$:

If $\Delta u(x,y)$ is positive, that cell is sitting in a local valley (its neighbours sum to more than $4u(x,y)$).
If $\Delta u(x,y)$ is negative, it’s on a local peak (its neighbours sum to less).
If $\Delta u(x,y) = 0$, it’s perfectly “flat” with respect to its neighbors.

By adding a small multiple of $\Delta u$ to the cell’s value,

u_new = u_old + D·dt·Δu

you’re doing:

Valleys ($Δu>0$) get nudged up toward their neighbors.
Peaks ($Δu<0$) get pushed down toward their neighbors.
Flat spots ($Δu=0$) stay where they are.

That’s exactly how you spread (diffuse) material: peaks shrink and valleys fill in.

2. Conservation of total quantity

Notice that for each interior cell the 5‑point stencil $ L + R + U + D - 4C $ sums to zero if you sum it over the whole grid (with periodic or zero‑flux boundaries). In other words, you’re only redistributing values between neighbours, never creating or destroying total mass.

3. Continuous intuition

The continuous heat (diffusion) equation, $ \frac{\partial u}{\partial t} = D,\nabla^2 u, $ says “the rate of change at a point is proportional to its local curvature.”

Positive curvature (concave up) → value rises.
Negative curvature (concave down) → value falls.

Discretizing $\nabla^2$ by your 5‑point stencil is exactly how you translate that “curvature‐driven flow” to a grid.

4. Energy‐minimization view

Diffusion also minimizes the field’s roughness. You can think of the sum of squared differences between neighbours as an “energy”: $ E = \sum_{\langle i,j\rangle}(u_i - u_j)^2. $ Gradient‐descent on $E$ leads you to update each $u_i$ in the direction of minus the discrete Laplacian. Thus adding $+\Delta u$ drives the system toward the smoothest configuration.

Bottom line

Laplacian measures how a cell differs from its neighbours.
Adding it pushes peaks down, valleys up.
Over many small timesteps, you get the familiar “smoothing out” of diffusion.

No mysterious magic—just local averaging enforced by the Laplace stencil at each step!

Below is a side‑by‑side comparison of two simple Unity shaders—one driving diffusion by repeatedly applying the discrete Laplacian, the other performing a single‑pass Gaussian blur. You’ll see they both are just weighted sums of neighboring pixels, but with different kernels.

Diffusion Shader (Compute)	Gaussian Blur Shader (Fragment)
```hlsl	```hlsl
// Compute shader: one Jacobi iteration of diffusion	// Simple 5×5 separable Gaussian blur in a fragment shader
#pragma kernel CSMain	Shader “Custom/GaussianBlur” {
RWTexture2D Result;	Properties {
Texture2D Source;	_MainTex (“Texture”, 2D) = “white” {}
sampler SourceSampler;	_Radius(“Blur Radius (px)”, Range(1,10)) = 2
float D_dt; // pre‑multiplied D * dt	}
	SubShader {
[numthreads(8,8,1)]	Pass {
void CSMain (uint3 id : SV_DispatchThreadID) {	CGPROGRAM
int2 uv = id.xy;	#pragma vertex vert
float C = Source.Load(int3(uv,0));	#pragma fragment frag
float L = Source.Load(int3(uv + int2(-1,0),0));	#include “UnityCG.cginc”
float R = Source.Load(int3(uv + int2(+1,0),0));
float U = Source.Load(int3(uv + int2(0,+1),0));	sampler2D _MainTex;
float D = Source.Load(int3(uv + int2(0,-1),0));	float _Radius;

// discrete Laplacian	float4 frag(v2f i) : SV_Target {
float lap = (L + R + U + D) - 4 * C;	float2 uv = i.uv;
float updated = C + D_dt * lap;	// horizontal blur
	float4 sum = float4(0,0,0,0);
Result[id.xy] = updated;	float w[5] = {0.204164, 0.304005, 0.193068, 0.072302, 0.016461};
}	for (int x=-2; x<=2; x++)
}	sum += w[abs(x)] * tex2D(_MainTex, uv + float2(x,0)/_ScreenParams.xy);
```
	// vertical blur
Dispatch via	float4 sum2 = float4(0,0,0,0);
```csharp	for (int y=-2; y<=2; y++)
computeShader.SetTexture(0, “Source”, srcRT);	sum2 += w[abs(y)] * sum;
computeShader.SetFloat(“D_dt”, diffusionRate * dt);	return sum2;
computeShader.SetTexture(0, “Result”, dstRT);	}
computeShader.Dispatch(W/8, H/8, 1);	ENDCG
```	}
	}
Key points:	Key points:
- Uses a 5‑point stencil every iteration to push values	- Uses a fixed Gaussian kernel (5 taps) in two passes
up/down toward neighbor average.	- Single‑pass convolution approximates one time‑step of
- You’d dispatch many times to see visible “smoothing” over time.	diffusion with a Gaussian kernel.
- Totally local: only reads immediate N/S/E/W pixels.	- Reads a wider footprint but in a separable way for speed.
- Conserves total “mass” each pass (just redistributes).	- Weights sum to 1, so overall brightness stays constant.

Gaussian blur = diffusion at “infinite time”
Solving the heat equation analytically for time $t$ gives convolution with a Gaussian of variance proportional to $t$.
Discrete diffusion ≈ iterative Gaussian
If you apply the diffusion shader many small steps, you approximate a Gaussian blur whose standard deviation grows as $\sqrt{D,t}$.
Stencil vs. kernel
Diffusion uses a tiny 5‑point stencil iterated in time. Gaussian blur precomputes the combined effect and does it in one (or two separable) passes.

Feel free to tweak the diffusion $\Delta t$ or the Gaussian radius and see how many diffusion steps match a given blur radius!

You’d pick diffusion (iterated Laplacian) over a one‑shot Gaussian blur whenever you need one or more of these qualities:

Progressive control & temporal coupling
- Animated smoothing: If you want the blur to “grow” over time—e.g. an effect that starts sharp and gradually softens—you can simply run more diffusion steps each frame, without re‑computing a bigger Gaussian kernel.
- Feedback systems: In fluid sims or reaction‑diffusion, the field at time t feeds back into your next step. Gaussian is a static filter, diffusion is a dynamic process.
Non‑linear or anisotropic extensions
- Edge‑preserving smoothing: By replacing the constant Laplacian with a nonlinear or anisotropic operator (e.g. Perona–Malik), you diffuse within regions but not across strong edges. You can’t do that with a fixed Gaussian.
- Curvature flows: Mean‑curvature motion, surface fairing, level‑set evolution—these are all PDE‑based and built on the Laplacian kernel.
Custom boundary conditions & geometry
- Dirichlet/Neumann control: You can pin or mirror values at arbitrary boundaries (e.g. obstacles, masks) each diffusion step. A Gaussian convolution layer can’t easily honor per‑pixel boundary rules.
- Irregular domains: On a mesh or unstructured grid you generalize the Laplacian stencil; Gaussian assumes a rectangular grid.
Memory & performance trade‑offs
- Tiny kernel footprint: A 5 point stencil uses minimal cache per pass. If you only need a tiny amount of smoothing (1–3 steps), diffusion can cost less memory bandwidth than a 5×5 or larger separable blur.
- Scalable blur radius: To get a really wide Gaussian (radius > 10 px), you’d need many taps or passes; diffusion gets you arbitrary “radius” simply by more iterations (though at linear time cost).
Physical realism
- True heat‑equation behavior: When simulating physical diffusion (e.g. temperature, chemical concentration, incompressible fluids), you want the field to follow ∂u/∂t = D∇²u exactly. A Gaussian is just the analytic solution at one fixed time—no notion of intermediate “snapshots.”

When Gaussian still wins

One‑off image blur: If you just need a static blur of known radius in a post‑process pass, Gaussian is more efficient (separable convolution, hardware‑accelerated sampling).
Fixed, predictable kernel: Gaussian weights give you direct control over your blur shape in one pass—very handy for UI blur, bloom, depth‑of‑field.

Summary:
Use diffusion if you need time‑dependent smoothing, nonlinear/aniso diffusion, custom boundary rules, or are actually simulating a physical PDE. Use a Gaussian blur for one‑shot, efficient, predictable image filtering when you don’t need the extra flexibility.

Diffusion curves are vector graphic primitives for creating smooth‑shaded images by defining distinct color constraints along curves and propagating these colors across the canvas via a diffusion process solved as a Poisson equation, resulting in compact, resolution‑independent vector representations that support complex gradients and real‑time GPU rendering citeturn0search1turn2view0.

Historical Background

Diffusion curves were introduced by Orzan et al. in 2008 in an ACM Transactions on Graphics paper titled “Diffusion Curves: A Vector Representation for Smooth‑Shaded Images” citeturn0search4turn4search5.

What Are Diffusion Curves?

Diffusion curves partition the 2D graphics space with vector paths, assigning color values to each side of a curve and specifying a blur transition across the curve boundary citeturn0search1.
The final image is obtained by solving a global Poisson equation whose constraints come from the gradients across all diffusion curves, enabling smooth interpolation of colors between curves citeturn2view0.

Core Data Structure

Each diffusion curve is represented as a cubic Bézier spline defined by a set of geometric control points citeturn2view0.
Attributes include two sets of color control points (one for each side of the curve) and blur control points that define the sharpness of the color transition, all linearly interpolated along the spline citeturn2view0.

Rendering Process

Rendering starts by rasterizing color sources: pixels slightly offset from the curve on each side encode color constraints, alongside a gradient field to enforce sharp transitions on the curve citeturn2view0.
A GPU‑based iterative solver then performs diffusion by solving ∇²I = div w per color channel, propagating the constrained colors across the image domain citeturn2view0.
Finally, a spatially varying re‑blur step guided by the curve’s blur attributes smooths the result to the artist’s specification citeturn2view0.

Automatic Curve Extraction

Bitmap images can be vectorized into diffusion curves by detecting edges (e.g., with Canny), fitting Bézier splines via active‑contour snapping, and sampling colors at a fixed normal distance from the curve citeturn2view0.
Outlier color samples are removed using local standard‑deviation checks, and attribute control points are fitted with the Douglas–Peucker polyline‑simplification algorithm to approximate original image gradients citeturn2view0.
Advances in inverse diffusion‑curve optimization further refine curve geometry and color attributes to closely match input imagery citeturn0search11.

Advantages and Applications

Diffusion curves offer a sparse, feature‑based representation that avoids the dense control‑point grids of gradient meshes while maintaining resolution independence and support for keyframe animation citeturn2view0.
Artists gain intuitive freehand workflows: curves need not form closed loops, may intersect arbitrarily, and can carry varying color and blur attributes for rich control citeturn2view0.
Inkscape integrates diffusion curves as a gradient primitive, enabling freeform gradient editing with per‑stroke blur customization citeturn0search12.
Diffusion curves also underpin resolution‑independent texture mapping in graphics pipelines, offering compact, scalable representations of complex shading for real‑time rendering citeturn0search16.

Extensions and Future Directions

Controlled‑diffusion research introduces diffusion strength, anisotropic flow, and barrier curves that block diffusion without emitting color, giving artists greater expressivity citeturn3view0.
Efficient GPU Laplacian solvers have been proposed for real‑time Poisson‑based editing and rendering with diffusion curves, improving performance and scalability citeturn6view0.
Nonlinear and anisotropic diffusion extensions enable edge‑preserving smoothing and flow‑field‑guided color propagation for advanced non‑photorealistic rendering effects citeturn3view0.

Conclusion

By embedding PDE‑based diffusion into vector primitives, diffusion curves bridge raster and vector worlds, combining compact editability with photorealistic shading capabilities citeturn2view0.
Their integration into design tools and ongoing research into controlled diffusion and solver optimizations promise ever more flexible and powerful vector graphics workflows in the future citeturn3view0.

Summary

Diffusion curves solve the Laplace equation, ∇²I = 0, across an image domain with boundary conditions defined by user‑drawn curves, and discretize this using a Jacobi‑like stencil implemented on the GPU citeturn3search10. Unity’s compute shader system (.compute assets with HLSL and ShaderLab) supports this by defining Texture2D and RWTexture2D resources, marking kernels with #pragma kernel, and parallelizing workgroups via the [numthreads(x,y,z)] directive citeturn0search0turn0search5. At each iteration, the compute kernel loads neighbor texels via Load, sums them with a variable stencil radius for faster convergence, and writes results to a ping‑ponged output texture, while constraint pixels are simply re‑written from a boundary color map to preserve curve definitions citeturn4view0. A C# controller script creates two RenderTextures, loads the ComputeShader, sets textures and parameters with SetTexture and SetInt, dispatches multiple iterations via Dispatch, and swaps the textures each frame to converge to the smooth‑shaded result citeturn2search0turn0search1. Typically, 8–16 Jacobi passes with a shrinking stencil size produce visually smooth outputs at interactive rates, outperforming multigrid approaches in both simplicity and temporal stability citeturn4view0.

Compute Shader Implementation (`DiffusionCurves.compute`)

#pragma kernel JacobiIteration
#pragma target 5.0

// Previous iteration’s image
Texture2D<float4>    _PrevTex;
// Output (next iteration)
RWTexture2D<float4>  _ResultTex;
// Boundary colors along curves
Texture2D<float4>    _BoundaryColorTex;
// Mask: 1 for curve constraint pixels, 0 otherwise
Texture2D<uint>      _ConstraintMaskTex;
// Variable stencil radius (shrinks per iteration if desired)
int                  _StencilRadius;

[numthreads(8,8,1)]
void JacobiIteration(uint3 id : SV_DispatchThreadID) {
    int2 uv = id.xy;
    // 1) Preserve curve constraints exactly
    if (_ConstraintMaskTex[uv] == 1) {
        _ResultTex[uv] = _BoundaryColorTex[uv];
    }
    else {
        // 2) Sum N/S/E/W neighbors at distance _StencilRadius
        float4 sum = float4(0,0,0,0);
        int r = _StencilRadius;
        sum += _PrevTex.Load(int3(uv + int2( r, 0), 0));
        sum += _PrevTex.Load(int3(uv + int2(-r, 0), 0));
        sum += _PrevTex.Load(int3(uv + int2(0,  r), 0));
        sum += _PrevTex.Load(int3(uv + int2(0, -r), 0));
        // 3) Average
        _ResultTex[uv] = sum * 0.25;
    }
}

Stencil choice: Using a variable radius (instead of the classic 1‑pixel stencil) accelerates color transport across the image, as shown by Jeschke et al. citeturn4view0.
Thread groups: [numthreads(8,8,1)] lets Unity batch pixels into 8×8 tiles for efficient GPU execution citeturn0search0.

C# Dispatch Script (`DiffusionCurvesController.cs`)

using UnityEngine;

public class DiffusionCurvesController : MonoBehaviour {
    public ComputeShader diffusionShader;
    public Texture2D     boundaryColorTex;   // RGB boundary constraints
    public Texture2D     constraintMaskTex;  // R=1 at curve pixels, else 0
    public int           iterations    = 12;
    public int           stencilRadius = 1;

    RenderTexture texA, texB;
    int kernelID;

    void Start() {
        int w = boundaryColorTex.width;
        int h = boundaryColorTex.height;
        // Create ping‑pong RenderTextures with random write
        texA = MakeRT(w, h);
        texB = MakeRT(w, h);
        // Initialize with boundary colors so curves “seed” the solver
        Graphics.Blit(boundaryColorTex, texA);
        // Grab the kernel index
        kernelID = diffusionShader.FindKernel("JacobiIteration");
    }

    RenderTexture MakeRT(int w, int h) {
        var rt = new RenderTexture(w, h, 0, RenderTextureFormat.ARGBFloat);
        rt.enableRandomWrite = true;
        rt.Create();
        return rt;
    }

    void Update() {
        // Bind static resources
        diffusionShader.SetTexture(kernelID, "_BoundaryColorTex", boundaryColorTex);
        diffusionShader.SetTexture(kernelID, "_ConstraintMaskTex", constraintMaskTex);
        diffusionShader.SetInt   ("_StencilRadius",     stencilRadius);

        // Perform iterative Jacobi passes
        for (int i = 0; i < iterations; i++) {
            diffusionShader.SetTexture(kernelID, "_PrevTex",   texA);
            diffusionShader.SetTexture(kernelID, "_ResultTex", texB);
            // Dispatch size = ceil(width/8) × ceil(height/8)
            int gx = Mathf.CeilToInt(boundaryColorTex.width  / 8f);
            int gy = Mathf.CeilToInt(boundaryColorTex.height / 8f);
            diffusionShader.Dispatch(kernelID, gx, gy, 1);
            // Swap buffers
            var tmp = texA; texA = texB; texB = tmp;
            // Optionally: stencilRadius = ComputeNextRadius(i);
        }

        // Now texA holds the converged diffusion curve image.
        // Assign it to a material or UI element as needed.
    }
}

RenderTexture settings: enableRandomWrite = true allows compute shaders to write into them citeturn0search5.
Ping‑pong: Swapping texA and texB preserves intermediate results without extra allocations citeturn0search1.

Performance & Extensions

Iteration count: Jeschke et al. report that 8–10 variable‑stencil Jacobi iterations achieve a smooth look on 800×800 images at >100 Hz on a mid‑range GPU, with stable results under animation citeturn4view0.
Gradient‑domain variants: For Poisson blending (non‑homogeneous ∇²I = div G), one replaces the uniform 0 Laplacian with a divergence of a guidance field (e.g. image gradients), summing directional weights in the compute kernel citeturn5search0turn5search1.
Seamless cloning in Unity: The open‑source project asus4/unity-poisson-blending demonstrates a full Poisson blend in Unity using a very similar ComputeShader approach, processing hundreds of iterations in milliseconds citeturn2search0.
Alternative GPU methods: Cheind’s poisson-image-editing C++ library and TheArtSquirrel’s Poisson‑blend blog illustrate CPU and GPU reference implementations for seamless cloning, which can inform advanced Unity integrations citeturn5search4turn5search6.

Conclusion

By leveraging Unity’s compute shader pipeline, a diffusion‑curve solver becomes a few dozen lines of HLSL and C#, running entirely on the GPU for real‑time, resolution‑independent vector‑style gradients. Variable stencil sizes dramatically speed convergence, and the same pattern extends naturally to gradient‑domain editing tasks like Poisson blending. This approach both preserves the artistic flexibility of diffusion curves and taps the massive parallelism of modern GPUs for interactive performance.

A practical 2D solver. The inputs to our solver are the Voronoi color image and the distance map obtained in the rasterization step (section 3.2). The former serves as an initial solution guess, and the latter determines the radius of the sampling circle around each pixel, which ensures that no boundaries are crossed. (We reduce all of the radii by 8% to compensate for errors in the distance map, as described in section 3.2.) Averaging all samples from a large circle during diffusion steps is clearly too slow for practical use, so in practice we reduced the number of samples to only four values in and difficult to see even with a side-by-side comparison to the converged solution. Figure 3 shows an example result of SA over a number of iterations. Note the very subtle differences between the solution with 8 steps and the converged solution with 400. The SH strategy is similar to SA, except that the sampling radii are kept at full size for the first half of the diffusion iterations. SH converges towards the correct result faster than SA, but visible banding artifacts may remain for up to 12 or 14 iterations. Figure 5 compares the kinds of errors that result after 8 iterations with different shrinking strategies. At this number of iterations, SH exhibits less than half the error of SA. Even so, the proper shrinking strategy for a given application depends on the desired accuracy and smoothness of the result. Of course, both SA and SH converge to the right solution with an increasing number of steps, simply by the fact that the last iterations are identical to standard diffusion. In addition, in-between solutions are visually appealing because all errors have very low frequency. Figure 5 demonstrates that this type of error is hardly visible for the human eye: the left image with smaller but high frequency error looks less appealing. We also maintain temporal coherence, i.e., flickering is avoided as can be seen in the accompanying video.axis aligned directions. Unfortunately, since each sample now has a great influence on the result, artifacts similar to mach banding can result. Figure 5 (left) shows such horizontal and vertical banding artifacts, mostly visible at the tips of the Siggraph logo. Our solution is to successively shrink the stencil radius during the diffusion. We have identified two different shrinking strategies suitable for different situations, which we will refer to as “shrink always” (SA) and “shrink half” (SH). In SA, we shrink the stencil size linearly at each step. In other words, if the solver will perform n iterations, the radius on iteration i will be scaled by 1 − i n . This lets pixel values align more and more with their local neighborhood at later iterations, producing smoother behavior. The appeal of SA is that it removes visible banding artifacts and color plateaus quickly, almost always within 8 iterations. However, the color distribution is typically not converged at that point, which makes the image appear with slightly higher contrast. This difference is very minor, though,

This response presents a Unity compute shader implementation of the practical 2D diffusion‑curve solver by Jeschke et al., which initializes from a Voronoi color image and a distance map (with radii reduced by 8% to prevent boundary crossing), performs four‑directional sampling per pixel, and incorporates two variable‑stencil shrinking strategies—Shrink Always (SA) and Shrink Half (SH)—to achieve visually smooth convergence in just a few Jacobi‑like iterations citeturn1view0turn0search2.

Overview of the Practical Solver

The solver begins with an initial guess obtained by rasterizing diffusion curves into a Voronoi diagram, coloring each pixel by its nearest curve segment citeturn1view0.
A distance map records, for every pixel, the distance to the closest boundary curve, ensuring sampling circles never cross region boundaries citeturn1view0.

Key Inputs

Voronoi Color Image

The Voronoi color image serves both as the starting solution and as a Dirichlet boundary condition that remains fixed during diffusion iterations citeturn1view0.

Distance Map

The distance map is produced by rendering “slanted” polygons and using depth comparisons; it has an inherent up-to-8 % error due to the fan approximation, which the solver compensates for by scaling radii by 0.92 citeturn1view0.

Sampling Strategy

Rather than computing a full circular average (which is too costly), the solver samples only four axis‑aligned neighbors—up, down, left, and right—per pixel citeturn1view0.
At each iteration, each interior pixel is updated to the average of those four samples, taken at the current stencil radius determined by the distance map and shrinking strategy citeturn1view0.

Shrinking Strategies

Shrink Always (SA)

In SA, the stencil radius at iteration i is linearly scaled by $1 - \tfrac{i}{n}$, where $n$ is the total number of iterations citeturn1view0.
SA typically eliminates banding artifacts within about eight iterations, though full convergence to the exact minimal surface may require many more passes citeturn1view0.

Shrink Half (SH)

SH maintains the full sampling radius for the first $n/2$ iterations and then linearly shrinks it over the remaining steps, accelerating overall convergence but allowing minor banding for up to 12–14 iterations citeturn1view0.
Figure 5 on ResearchGate illustrates that, after eight iterations, SH reduces error by more than half compared to SA citeturn0search4.

Unity Compute Shader Implementation

Compute Shader (HLSL)

The shader declares a read‑only Texture2D<float4> for the current solution and a RWTexture2D<float4> for the next iteration, plus a Texture2D<float> for the distance map and a Texture2D<uint> mask for Dirichlet boundaries citeturn2search0.
A [numthreads(8,8,1)] directive groups threads into 8×8 tiles (64 threads each), balancing parallel throughput and memory access citeturn2search2.
Inside the kernel, if the boundary mask is set, the pixel is reset to its original Voronoi color; otherwise, four neighbor samples are fetched at offsets based on the computed radius for the current iteration and strategy, and their average is written out citeturn2search0.

#pragma kernel CSMain
#pragma target 5.0

Texture2D<float4> _PrevTex;
RWTexture2D<float4> _ResultTex;
Texture2D<float4> _VoronoiTex;
Texture2D<float>   _DistanceMap;
Texture2D<uint>    _BoundaryMask;
int                _Iterations;
int                _Iteration;
int                _Strategy;    // 0 = SA, 1 = SH

[numthreads(8,8,1)]
void CSMain(uint3 id : SV_DispatchThreadID) {
    int2 uv = id.xy;
    if (_BoundaryMask.Load(int3(uv,0)) == 1) {
        _ResultTex[uv] = _VoronoiTex.Load(int3(uv,0));
    } else {
        float d0 = _DistanceMap.Load(int3(uv,0));
        float baseR = d0 * 0.92;  // compensate 8% rasterization error
        float t    = _Iteration / (float)_Iterations;
        float r    = (_Strategy == 0)
                   ? baseR * (1 - t)
                   : (t < 0.5 ? baseR : baseR * (1 - 2*(t - 0.5)));
        int   ir   = max(1, (int)r);

        float4 sum = float4(0,0,0,0);
        sum += _PrevTex.Load(int3(uv + int2( ir,  0), 0));
        sum += _PrevTex.Load(int3(uv + int2(-ir,  0), 0));
        sum += _PrevTex.Load(int3(uv + int2( 0,  ir), 0));
        sum += _PrevTex.Load(int3(uv + int2( 0, -ir), 0));
        _ResultTex[uv] = sum * 0.25;
    }
}

C# Controller Script

A Unity C# script loads the ComputeShader, creates two RenderTexture buffers with enableRandomWrite=true for ping‑ponging, and initializes one with the Voronoi color image for the first pass citeturn2search4.
Each frame, it sets shader parameters (_DistanceMap, _BoundaryMask, _VoronoiTex, _Iterations, _Strategy, _Iteration), dispatches via Dispatch(kernelID, width/8, height/8, 1), then swaps the buffers for the next iteration citeturn3search0.

using UnityEngine;

public class DiffusionSolverController : MonoBehaviour {
    public ComputeShader solver;
    public Texture2D     voronoiTex;
    public Texture2D     distanceMap;
    public Texture2D     boundaryMask;
    public int           iterations = 12;
    public int           strategy   = 0; // 0 = SA, 1 = SH

    RenderTexture texA, texB;
    int kernelID;

    void Start() {
        int w = voronoiTex.width, h = voronoiTex.height;
        texA = MakeRT(w,h); texB = MakeRT(w,h);
        Graphics.Blit(voronoiTex, texA);
        kernelID = solver.FindKernel("CSMain");
    }

    RenderTexture MakeRT(int w,int h) {
        var rt = new RenderTexture(w,h,0,RenderTextureFormat.ARGBFloat);
        rt.enableRandomWrite = true; rt.Create();
        return rt;
    }

    void Update() {
        solver.SetTexture(kernelID, "_PrevTex",   texA);
        solver.SetTexture(kernelID, "_ResultTex", texB);
        solver.SetTexture(kernelID, "_VoronoiTex",voronoiTex);
        solver.SetTexture(kernelID, "_DistanceMap",distanceMap);
        solver.SetTexture(kernelID, "_BoundaryMask",boundaryMask);
        solver.SetInt   ("_Iterations", iterations);
        solver.SetInt   ("_Strategy",   strategy);

        for (int i = 0; i < iterations; i++) {
            solver.SetInt("_Iteration", i);
            int gx = Mathf.CeilToInt(voronoiTex.width  / 8f);
            int gy = Mathf.CeilToInt(voronoiTex.height / 8f);
            solver.Dispatch(kernelID, gx, gy, 1);
            (texA,texB) = (texB,texA);
        }

        // texA now contains the final diffused image
    }
}

Discussion

Performance Considerations

Typically eight iterations suffice to produce visually smooth output, and on a GeForce 8800 GTX a 1024×1024 distance map plus seven diffusion steps took only ~7.6 ms (≈130 Hz) citeturn1view0.

Choosing SA vs SH

Use SA when you need minimal banding artifacts within few iterations, and SH when you prefer faster overall convergence and can tolerate brief low‑frequency banding citeturn1view0.

At its core, diffusion curves rely on specifying Dirichlet constraints (fixed color values) along geometric primitives and then solving a global Poisson (diffusion) equation to smoothly interpolate those colors across the image. You can re-use exactly the same machinery by treating your colored pixels as zero-dimensional “curves” (i.e. point constraints) instead of Beziér splines. In practice this means:

Seed your image: create a mask or texture where only your colored pixels carry their RGB values, and all other pixels are uninitialized (or zero).
Mark constraints: build a binary mask that flags exactly those seed pixels as “boundary” so they remain fixed.
Solve the Poisson equation: run your diffusion-curve/Poisson solver (Jacobi, multigrid, FFT, whatever) over the entire domain; at each iteration, you re-impose the seed colors and update everything else by averaging neighbors (the discrete Laplacian stencil).

All the benefits of diffusion curves—resolution independence, arbitrary seed placement, smooth interpolation, easy GPU implementation—carry over with no curves at all!

1. From Curves to Points: Dirichlet Constraints Anywhere

Diffusion curves (Orzan et al., 2008) define colors on either side of each curve and solve
$ \nabla^2 I = 0 $ with those curve-side values as boundary conditions citeturn1search10. But the mathematics is identical if your “curves” degenerate to single pixels: each seed pixel simply pins a fixed color into the solver. You still solve a homogeneous Poisson problem; you just supply a constraint mask of isolated points instead of splines citeturn1search10.

2. Prior Art: Seeded Diffusion & Poisson Editing

Seeded segmentation via diffusion treats labeled pixels as heat sources and diffuses temperature outwards, using the same 4-neighbour stencil you know from diffusion curves citeturn0search0.
Poisson Image Editing (Perez et al., 2003) fixes pixel values inside a region Ω and solves ∇²f = div v to seamlessly blend colors—a direct analogue to seeding arbitrary pixels citeturn0search1.
Colorization methods in US Patent US 8913074 B2 likewise formulate a cost that penalizes deviation from seed colors vs. neighbor differences, then minimize it with a Poisson-style solver citeturn0search17.
Scribble-based vector colorization extends diffusion curves by letting the user drop colored strokes (points) and diffusing them to create smooth-shaded vector art citeturn1search3.

3. Fast Multipole & Point-Based Extensions

The “Fast Multipole Representation of Diffusion Curves and Points” shows that you can treat both curves and isolated point sources in a unified framework, evaluating the Poisson solution at any pixel via a fast hierarchical summation of Green’s functions citeturn0search12. This bypasses iterative solvers entirely and supports millions of seeds in real time.

4. Implementation Recipe (Unity Compute Shader)

You already have a diffusion‐curve shader. To adapt it:

Textures
- _PrevTex: floating-point RGBA image for the current solution
- _SeedColorTex: same size, zero everywhere except seeds carry their RGB
- _SeedMaskTex: uint mask (1 at seeds, 0 elsewhere)

Compute Kernel

if (_SeedMaskTex[uv] == 1) {
    _ResultTex[uv] = _SeedColorTex[uv];
} else {
    float4 sum = Load(_PrevTex, uv+left)
               + Load(_PrevTex, uv+right)
               + Load(_PrevTex, uv+up)
               + Load(_PrevTex, uv+down);
    _ResultTex[uv] = sum * 0.25;
}

Ping-Pong & Iterate
Dispatch ~8×8 threads per tile, iterate 8–16 passes (more for higher smoothness) citeturn1search10.

Because seeds are only single pixels, this converges to a harmonic interpolation filling the domain with the smoothest blend of your seed colors—exactly what you want when “spreading” a handful of colored pixels citeturn1search10.

5. When to Use This vs. Gaussian or Other Methods

Use seeded diffusion when you need exact color anchors at arbitrary pixels, physical realism (exact Laplace solution), or non-linear/anisotropic extensions (edge-aware diffusion).
Use Gaussian blur when you simply want a uniform blur with no fixed seeds or you need hardware-accelerated separable filters in one pass.

6. Further Reading & Tools

Orzan et al., Diffusion Curves (ACM TOG 2008): original vector primitive description citeturn1search10
Perez et al., Poisson Image Editing: gradient-domain edits & color propagation citeturn0search1
US 8913074 B2, colorization patent: constrained optimization for seed-based color spread citeturn0search17
Hongbo Fu, Scribble-based Colorization: diffusion strokes in vector graphics citeturn1search3
Columbia Fast Multipole: point and curve unified Poisson solver citeturn0search12
IPOL, Poisson Image Editing: open-source reference implementation citeturn0search7
Chris Tralie, Intuitive Poisson Blending: tutorial and code citeturn0search15
Patch Transform (MIT): alternative patch-based editing via PDEs citeturn0search18
PVG, Poisson Vector Graphics: adds point&region primitives to diffusion curves citeturn1search8

With those references and the simple shader tweak above, you can “spread” any set of colored pixels exactly as diffusion curves do—without drawing a single curve.