A Variational Screened Poisson Reconstruction for Whole-Slide Stain Normalization

1. Introduction

1.1. Background

Histopathological analysis is the current gold standard for disease diagnosis, particularly in oncology. With the rapid evolution of digital pathology [1], the slides containing stained tissue sections are routinely digitized into ultra–high-resolution whole-slide images (WSIs) using advanced microscopy scanners. As a result, computational analysis of digital pathology images is increasingly complementing, and in some cases replacing, manual microscopic evaluation by pathologists. Within this paradigm, the quality, consistency, and standardization of digitized pathological images constitute pivotal determinants of diagnostic reproducibility and the performance of subsequent computational analyses.

The WSIs imaging pipeline begins with conventional tissue preparation, followed by whole-slide scanning into standard RGB image formats. Figure 1 provides an overview of the complete acquisition workflow [28].

Clinically, tissue specimens(obtained via surgical resection or biopsy)are fixed in formalin to preserve structural integrity, dehydrated through graded ethanol solutions, and embedded in paraffin wax to form a block. Thin tissue sections (3–6 µm) are cut from this block, mounted onto slides, and stained, most commonly with hematoxylin and eosin (H&E). Finally, a whole-slide scanner digitizes the slide at subcellular resolution to generate a WSI that is stored in a gigapixel images.However, deviations inevitably accumulate due to heterogeneous tissue preparation (e.g., sectioning, reagents, pH) and diverse digitization conditions (e.g., scanner optics, sensor noise).These artifacts introduce substantial domain shifts that impair diagnostic consistency and algorithmic performance. To address this, robust stain normalization is required to harmonize WSIs from heterogeneous sources, thereby standardizing the input for reliable downstream computational analysis.

1.2. Related Works

Research on stain normalization has advanced significantly since the digitalization of pathology. The seminal work of Reinhard et al. [6] introduced statistical color alignment via mean–variance matching in the decorrelated Lab space. However, this global statistics-based normalization inevitably results in the over-smoothing of fine image details.

To improve illumination robustness and chromatic fidelity, subsequent variations modified the transformation domain. Zarella et al. [14] and Magee et al. [15] adopted color projection in HSV/Lab spaces; Roy et al. [17] and Vijh et al. [18] introduced fuzzy-logic constraints; and Nadeem et al. [20] reformulated stain normalization as an optimal-transport problem via multimarginal Wasserstein barycenters. While these methods enhance perceptual smoothness, they remain fundamentally pixel-domain and global in nature, limiting their ability to preserve tissue microstructures.

Macenko et al. [7] shifted attention to structure preservation through stain vector estimation in optical-density (OD) space via singular value decomposition (SVD). Vahadane et al. [8] further incorporated sparse non-negative matrix factorization (SNMF), which jointly optimizes stain basis and concentration matrices. Adaptive deconvolution [19] and IDA decomposition [16] have since been proposed to address spatial heterogeneity while maintaining textural and morphological integrity. Gupta et al. [9] proposed a geometry-inspired, chemical- and tissue-invariant stain-normalization framework (GCTI). GCTI uses SVD/NMF and explicit geometric alignment to correct color variations caused by illumination, stain color vectors, and stain quantity.

Meanwhile, with the advancement of digital pathology, data-driven deep learning models have demonstrated superior performance to conventional stain normalization techniques in recent years. Existing approaches are mainly characterized by two classes: (1) supervised learning models that require paired or carefully curated training data [3,5,13,22,23,24,25], and (2) generative frameworks based on Generative Adversarial Networks (GANs) [4,11,12,21,26,33]. Despite achieving strong quantitative results, such improvements do not necessarily translate into better clinical diagnostic performance. When applied across institutions, sample-driven deep learning models are vulnerable to distribution shifts, which undermine reliability in real-world deployment. Moreover, their lack of explicit physical modeling of the staining process limits interpretability and hinders clinical acceptance.

1.3. Main Contributions

We propose the Screened Poisson Normalization (SPN) model, a variational framework that reformulates stain normalization as a strictly convex inverse problem governed by reaction-diffusion physics. This approach integrates physical interpretability with rigorous mathematical analysis, establishing the well-posedness and stability of the solution to guarantee the preservation of diagnostic morphology against input perturbations. To address the gigapixel scale of clinical pathology, we further derive a domain decomposition algorithm underpinned by the exponential decay of the screened Poisson Green’s function; this theoretical property ensures that interface errors remain exponentially localized, enabling the seamless and parallel reconstruction of Whole-Slide Images without the boundary artifacts typical of patch-based processing.

The paper is organized as follows. Section 2 investigates the mathematical nature of the stain normalization problem and establishes its theoretical foundations. Section 3 then derives Screened Poisson Normalization (SPN) from a reaction–diffusion interpretation of histological staining and formulates stain standardization as a variational inverse problem in the perceptual CIE $Lab$ space. Section 4 presents the rigorous analysis of the model, establishing well-posedness, Lipschitz stability, and the exponential localization property that motivates the non-overlapping subdomain implementation. Section 5 details the numerical realization, including the finite-difference discretization, the DCT-based spectral solver, and practical implementation considerations. Section 6 reports comprehensive experimental validation on multi-center datasets. Finally, Section 7 concludes the paper and discusses limitations and directions for future work.

2. Mathematical Foundations

We define stain normalization as an inverse problem on a bounded Lipschitz domain $\Omega \subset {\mathbb{R}}^2$. Let $I_S,I_T:\Omega \to {\mathbb{R}}^3$ denote the source and reference images, respectively. The task is to recover a normalized field $u:\Omega \to \mathbb{R}$ (processed channel-wise in the decorrelated CIE $Lab$ space [6]) that simultaneously satisfies structural fidelity to $I_S$ and statistical alignment with $I_T$. This section analyzes the ill-posed nature of this reconstruction and establishes its theoretical connection to reaction–diffusion physics.

2.1. Variational Formulation of the Inverse Problem

Mathematically, recovering u is fundamentally under-determined. The statistical constraint operator $\mathcal{C}\left(u\right)\approx \mathcal{C}\left(I_T\right)$ is non-injective, as infinite spatial configurations can share identical global moments. Conversely, relying solely on gradient fidelity $\nabla u\approx \nabla I_S$ determines the solution only up to an additive constant (the null space of the Neumann operator). Consequently, the unconstrained optimization is ill-posed in the sense of Hadamard [38].

To restore well-posedness, the problem must be cast within a Tikhonov regularization framework [39]. We seek a solution $u\in H^1(\Omega )$ that minimizes a strictly convex energy functional of the form:

$$\mathcal{J}\left(u\right):={\displaystyle \frac{1}{2}}{\parallel \nabla u-\mathbf{g}\parallel }_{L^2(\Omega )}^2+{\displaystyle \frac{\lambda }{2}}{\parallel u-T\parallel }_{L^2(\Omega )}^2,\lambda >0.$$

(1)

Here, the first term enforces structural consistency with a guidance field $\mathbf{g}$ derived from the source, while the second term anchors the solution to a chromatic prior T. The regularization parameter $\lambda$ balances these competing objectives and ensures the coercivity of the bilinear form, guaranteeing the existence of a unique, stable weak solution.

2.2. Physical Isomorphism: Reaction–Diffusion Kinetics

The variational structure of (1) is not merely a mathematical regularization; it is isomorphic to the physical transport laws [30] governing histological staining. The accumulation of dye within porous tissue follows reaction–diffusion kinetics [27,29,31]. For a dye species with concentration $C(\mathbf{x},t)$, the governing equation is:

$${\partial }_{t}C=\nabla ·(D\nabla C)+R(C,\mathbf{x}),$$

(2)

where D is the diffusion coefficient and R describes the binding reaction. Near chemical equilibrium, the binding process approximates linear kinetics $R\left(C\right)\approx -k(C-C_{\mathrm{eq}})$, where k is the reaction rate and $C_{\mathrm{eq}}$ is the local equilibrium concentration determined by tissue affinity.

At the thermodynamic steady state (${\partial }_{t}C=0$), assuming isotropic diffusion, (2) reduces to the screened Poisson equation:

$$-\Delta C+\kappa C=\kappa C_{\mathrm{eq}},\mathrm{with}\kappa :=k/D.$$

(3)

This physical equilibrium (3) shares the exact Euler-Lagrange form of the variational problem (1). The regularization parameter $\lambda$ thus corresponds to the physical ratio of reaction to diffusion rates ($\kappa$), while the anchor T corresponds to the chemical equilibrium potential. This theoretical equivalence provides the physical basis for modeling stain normalization as a screened Poisson process.

3. The Screened Poisson Normalization (SPN) Model

This section details the Screened Poisson Normalization (SPN) model. Operating in the CIE $Lab$ color space [6] and using optical density (OD) to represent dye concentration [33,34], we cast stain normalization as a variational reconstruction problem.

For each channel $c\in \{L^{*},a^{*},b^{*}\}$, we define the normalized field $u_c\in H^1(\Omega )$ as the unique minimizer of the strictly convex energy functional:

$${\mathcal{J}}_c\left(u_c\right):={\int }_{\Omega }\left({\displaystyle \frac{1}{2}}{\parallel \nabla u_c-{\mathbf{g}}_c\parallel }^2+{\displaystyle \frac{{\lambda }_c}{2}}{(u_c-T_c)}^2\right)d\mathbf{x}.$$

(4)

The first term ensures structural fidelity by penalizing gradient deviations from a guidance field ${\mathbf{g}}_c$, while the second acts as a Tikhonov regularizer, anchoring the solution to a chromatic prior $T_c$. The parameter ${\lambda }_c>0$ governs the trade-off between these objectives and ensures the problem is well-posed.

Minimizing (4) yields the associated Euler–Lagrange equation, which takes the form of a steady-state screened Poisson equation:

$$-\Delta u_c+{\lambda }_c{u}_c={\lambda }_c{T}_c-\nabla ·g_c\mathrm{in}\Omega ,$$

(5)

subject to the adiabatic boundary condition $(\nabla u_c-{\mathbf{g}}_c)·\mathbf{n}=0$ on $\partial \Omega$. Physically, this equation describes a diffusion process driven by a source term $f_c:={\lambda }_c{T}_c-\nabla ·{\mathbf{g}}_c$, which balances the injection of high-frequency morphological details with low-frequency chromatic information.

To close the system, we construct the chromatic anchor $T_c$ and guidance field ${\mathbf{g}}_c$ from the source $I_S$ and reference $I_T$. The anchor $T_c$ is defined by matching the global mean $\mu$ and standard deviation $\sigma$ of the reference:

$$T_c\left(\mathbf{x}\right)={\displaystyle \frac{{\sigma }_T^{\left(c\right)}}{{\sigma }_S^{\left(c\right)}}}\left({\left(I_S\right)}_c\left(\mathbf{x}\right)-{\mu }_S^{\left(c\right)}\right)+{\mu }_T^{\left(c\right)}.$$

(6)

This linear mapping provides a stable baseline for chromaticity that is robust to local variations. To preserve morphology, we parameterize the guidance field ${\mathbf{g}}_c$ as a spatially varying blend of source and anchor gradients:

$${\mathbf{g}}_c\left(\mathbf{x}\right)=(1-w\left(\mathbf{x}\right)){\alpha }_c\nabla {\left(I_S\right)}_c\left(\mathbf{x}\right)+w\left(\mathbf{x}\right)\nabla T_c\left(\mathbf{x}\right),$$

(7)

where ${\alpha }_c$ is a robust contrast factor defined by the ratio of median gradient magnitudes, i.e., ${\alpha }_c=median\parallel \nabla {\left(I_T\right)}_c\parallel /median\parallel \nabla {\left(I_S\right)}_c\parallel$. The structural weight $w\left(\mathbf{x}\right)\in (0,1]$ is derived from the source luminance channel $I_{S,L^{*}}$ via:

$$w\left(\mathbf{x}\right)=exp\left(-{\displaystyle \frac{\parallel \nabla I_{S,L^{*}}{\left(\mathbf{x}\right)\parallel }^2}{{\tau }^2}}\right),\tau >0.$$

(8)

In structurally salient regions ($w\to 0$), the guidance field strictly adheres to the source morphology $\nabla {\left(I_S\right)}_c$, whereas in homogeneous regions, the model relaxes toward the smooth chromatic anchor. The resulting solution corresponds to the steady state of a system with screening length ${\ell }_c\sim {\lambda }_c^{-1/2}$, preventing local artifacts from propagating globally.

4. Theoretical Analysis

This section analyzes the mathematical properties of the Screened Poisson Normalization (SPN) model. We treat the variational formulation as a linear inverse problem in the Sobolev space $H^1(\Omega )$, establishing the well-posedness of the solution, its stability with respect to input data, and the convergence of the domain decomposition scheme used for gigapixel processing.

The variational problem derived in Section 3 is formally equivalent to finding a weak solution $u\in H^1(\Omega )$ such that

$$B(u,v)=F\left(v\right),\forall v\in H^1(\Omega ),$$

(9)

where the bilinear form $B:H^1(\Omega )\times H^1(\Omega )\to \mathbb{R}$ and the linear functional $F:H^1(\Omega )\to \mathbb{R}$ are defined by:

$$B(u,v):={\int }_{\Omega }(\nabla u·\nabla v+\lambda uv)d\mathbf{x},F\left(v\right):={\int }_{\Omega }(\mathbf{g}·\nabla v+\lambda Tv)d\mathbf{x}.$$

We assume $\Omega \subset {\mathbb{R}}^2$ is a bounded Lipschitz domain, $\lambda >0$, $\mathbf{g}\in {\left[L^2(\Omega )\right]}^2$, and $T\in L^2(\Omega )$.

Theorem 1

(Existence and Uniqueness). For any screening parameter $\lambda >0$, guidance field $\mathbf{g}\in {\left[L^2(\Omega )\right]}^2$, and chromatic anchor $T\in L^2(\Omega )$, the variational problem (9) admits a unique solution $u^{*}\in H^1(\Omega )$.

Proof. 

The result follows from the Lax–Milgram theorem. First, the bilinear form B is continuous on $H^1(\Omega )\times H^1(\Omega )$ by the Cauchy–Schwarz inequality:

$$\left|B(u,v)\right|\le {\parallel \nabla u\parallel }_{L^2}{\parallel \nabla v\parallel }_{L^2}+{\lambda \parallel u\parallel }_{L^2}{\parallel v\parallel }_{L^2}\le {max(1,\lambda )\parallel u\parallel }_{H^1}{\parallel v\parallel }_{H^1}.$$

Second, B is strictly coercive due to the strictly positive screening parameter $\lambda$. Testing with $v=u$:

$$B(u,u)={\parallel \nabla u\parallel }_{L^2}^2+{\lambda \parallel u\parallel }_{L^2}^2\ge min(1,\lambda ){\parallel u\parallel }_{H^1}^2.$$

Since F defines a bounded linear functional on $H^1(\Omega )$, the existence of a unique minimizer is guaranteed.    □

With the solution uniquely defined, we quantify its stability with respect to perturbations in the input data $(\mathbf{g},T)$.

Theorem 2

(Lipschitz Stability). Let $u_1,u_2\in H^1(\Omega )$ be the unique weak solutions corresponding to data $({\mathbf{g}}_1,T_1)$ and $({\mathbf{g}}_2,T_2)$, respectively. The solution depends continuously on the data:

$$\parallel u_1-u_2{\parallel }_{H^1(\Omega )}\le {\displaystyle \frac{1}{min(1,\lambda )}}\left(\parallel {\mathbf{g}}_1-{\mathbf{g}}_2{\parallel }_{L^2(\Omega )}+\lambda {\parallel T_1-T_2\parallel }_{L^2(\Omega )}\right).$$

(10)

Proof. 

Let $w=u_1-u_2$. Subtracting the respective weak formulations yields

$$B(w,v)={\int }_{\Omega }\left(({\mathbf{g}}_1-{\mathbf{g}}_2)·\nabla v+\lambda (T_1-T_2)v\right)d\mathbf{x},\forall v\in H^1(\Omega ).$$

Choosing the test function $v=w$ and applying the coercivity of B alongside the Cauchy–Schwarz inequality gives

$${min(1,\lambda )\parallel w\parallel }_{H^1}^2\le B(w,w)\le \left(\parallel {\mathbf{g}}_1-{\mathbf{g}}_2{\parallel }_{L^2}+\lambda {\parallel T_1-T_2\parallel }_{L^2}\right){\parallel w\parallel }_{H^1}.$$

Dividing by ${\parallel w\parallel }_{H^1}$ completes the proof.    □

We next quantify the trade-off between morphological preservation and color normalization. The following estimate ensures that structural modifications are strictly bounded by the magnitude of the chromatic transfer.

Lemma 1

(Energy Comparison Principle). Let $u^{*}$ be the global minimizer of $\mathcal{J}\left(u\right)$. For any comparison function $v\in H^1(\Omega )$,

$$\parallel \nabla u^{*}{-\mathbf{g}\parallel }_{L^2(\Omega )}^2+\lambda \parallel u^{*}{-T\parallel }_{L^2(\Omega )}^2\le {\parallel \nabla v-\mathbf{g}\parallel }_{L^2(\Omega )}^2+\lambda {\parallel v-T\parallel }_{L^2(\Omega )}^2.$$

(11)

Proof. 

This is an immediate consequence of the optimality condition $\mathcal{J}\left(u^{*}\right)\le \mathcal{J}\left(v\right)$.    □

Theorem 3

(Structural Fidelity Bound). Assume $\mathbf{g}\approx \nabla I_S$. The $L^2$-deviation of the normalized gradients $\nabla u^{*}$ from the guidance field satisfies:

$$\parallel \nabla u^{*}{-\mathbf{g}\parallel }_{L^2(\Omega )}\le \sqrt{\lambda }\parallel I_S{-T\parallel }_{L^2(\Omega )}+{\parallel \nabla I_S-\mathbf{g}\parallel }_{L^2(\Omega )}.$$

(12)

Proof. 

Set $v=I_S$ in Lemma 1. Neglecting the non-negative term $\lambda \parallel u^{*}{-T\parallel }_{L^2}^2$ on the left-hand side yields

$$\parallel \nabla u^{*}{-\mathbf{g}\parallel }_{L^2}^2\le \parallel \nabla I_S{-\mathbf{g}\parallel }_{L^2}^2+\lambda {\parallel I_S-T\parallel }_{L^2}^2.$$

The result follows from the subadditivity of the square root function, $\sqrt{a^2+b^2}\le a+b$.    □

Finally, we address the convergence of the tiled reconstruction required for gigapixel Whole Slide Images. Let $\Omega$ be partitioned into non-overlapping subdomains $\left\{{\Omega }_k\right\}$. The local solution $u_{\mathrm{tile}}$ on ${\Omega }_k$ satisfies the natural boundary condition $(\nabla u_{\mathrm{tile}}-\mathbf{g})·\mathbf{n}=0$ on $\partial {\Omega }_k$. The discrepancy $e:=u_{\mathrm{global}}-u_{\mathrm{tile}}$ arises solely from the flux mismatch $\psi :=(\nabla u_{\mathrm{global}}-\mathbf{g})·\mathbf{n}$ at the interface $\Gamma :=\partial {\Omega }_k\setminus \partial \Omega$.

Lemma 2

(Green’s Function Asymptotics). The fundamental solution $G_{\lambda }\left(\mathbf{x}\right)$ of the operator $-\Delta +\lambda \mathcal{I}$ in ${\mathbb{R}}^2$ satisfies the decay estimate for $\left|\mathbf{x}\right|>0$:

$$|G_{\lambda }\left(\mathbf{x}\right)|\le C{\displaystyle \frac{e^{-\sqrt{\lambda }\left|\mathbf{x}\right|}}{\sqrt{\left|\mathbf{x}\right|}}}.$$

(13)

Theorem 4

(Exponential Localization of Interface Error). For any $\mathbf{x}\in {\Omega }_k$, the error $e\left(\mathbf{x}\right)$ decays exponentially with the distance to the interface Γ:

$$\left|e\left(\mathbf{x}\right)\right|\le Kexp\left(-{\displaystyle \frac{d(\mathbf{x},\Gamma )}{{\ell }_c}}\right),{\ell }_c={\lambda }^{-1/2}.$$

(14)

Proof. 

The error function satisfies $-\Delta e+\lambda e=0$ in ${\Omega }_k$ with Neumann data ${\partial }_{n}e=\psi$ on $\Gamma$. Using the boundary layer potential representation

$$e\left(\mathbf{x}\right)={\int }_{\Gamma }{G}_{\lambda }(\mathbf{x}-\mathbf{y})\psi \left(\mathbf{y}\right)d\sigma \left(\mathbf{y}\right),$$

and noting that $|\mathbf{x}-\mathbf{y}|\ge d(\mathbf{x},\Gamma )$ for all $\mathbf{y}\in \Gamma$, we invoke the monotonicity of the Green’s kernel:

$$\left|e\left(\mathbf{x}\right)\right|\le \left(\underset{r\ge d(\mathbf{x},\Gamma )}{sup}\left|G_{\lambda }\left(r\right)\right|\right){\parallel \psi \parallel }_{L^1(\Gamma )}.$$

Combining this uniform bound with the exponential decay of $G_{\lambda }$ yields the stated estimate.    □

Remark 1.

Theorem 4 implies that boundary artifacts are confined to a layer of characteristic width $\mathcal{O}\left({\ell }_c\right)$. For sufficiently large tiles, the interior solution converges exponentially to the global reconstruction, mathematically justifying the use of non-overlapping domain decomposition.

5. Optimization and Implementation

For each channel $c\in \{L^{*},a^{*},b^{*}\}$, the SPN model seeks a scalar field $U_c:\Omega \to \mathbb{R}$ that solves the screened Poisson equation (5),where $U_c$ denotes the discrete (grid-based) version of $u_c$.

5.1. Discrete and Non-Overlapping Domain Decomposition

The gigapixel resolution of WSI images renders global algorithms computationally prohibitive. Leveraging the localization properties of the SPN operator established in Section 4, we adopt a patch-wise implementation strategy.

We discretize the continuous domain $\Omega$ on a regular $H\times W$ lattice:

$${\Omega }_h=\{(x,y)\in {\mathbb{Z}}^2\mid 0\le x<H,0\le y<W\},$$

with unit grid spacing. The scalar field $U_c$ is represented as a matrix ${\mathbf{U}}_c\in {\mathbb{R}}^{H\times W}$, and the vector field ${\mathbf{g}}_c$ as a tensor in ${\mathbb{R}}^{H\times W\times 2}$. To avoid directional bias, we employ centered finite differences for the interior gradients:

$${\left({\nabla }_h{U}_c\right)}_{x,y}={\left[\begin{array}{c}{D}_x{U}_c\\ D_y{U}_c\end{array}\right]}_{x,y}=\left[\begin{array}{c}(U_{c,x+1,y}-U_{c,x-1,y})/2\\ (U_{c,x,y+1}-U_{c,x,y-1})/2\end{array}\right],1\le x\le H-2,1\le y\le W-2.$$

(15)

At the image boundaries, we utilize forward/backward differences consistent with the homogeneous Neumann condition ${\partial }_n{U}_c=0$.

The discrete divergence ${\nabla }_h·{\mathbf{g}}_c$ is defined as the negative adjoint of ${\nabla }_h$, ensuring that the discrete operator retains the symmetric positive-definite structure of the continuous formulation. Combining ${\nabla }_h$ and ${\nabla }_h·$ yields the standard five-point Laplacian stencil:

$${\left({\Delta }_h{U}_c\right)}_{x,y}=U_{c,x+1,y}+U_{c,x-1,y}+U_{c,x,y+1}+U_{c,x,y-1}-4U_{c,x,y},$$

(16)

with appropriate modifications at the boundaries. On the vectorized domain, this defines a sparse block-Toeplitz with Toeplitz-blocks (BTTB) matrix $\mathbf{L}\in {\mathbb{R}}^{N\times N}$ (where $N=HW$). Solving the linear system

$$(\mathbf{L}-{\lambda }_c\mathbf{I}){\mathbf{u}}_c={\mathbf{r}}_c$$

directly via sparse factorization is infeasible at the WSI scale. Instead, we exploit the spectral properties of $\mathbf{L}$ for efficient inversion.

5.2. Spectral Diagonalization and DCT Solver

Under homogeneous Neumann boundary conditions, the eigenfunctions of the discrete Laplacian (16) coincide with the basis functions of the 2D Type-II Discrete Cosine Transform (DCT-II). Let $\mathcal{F}$ denote the orthonormal forward 2D DCT-II operator and ${\mathcal{F}}^{-1}$ its inverse (DCT-III) [42]. For an input image $\mathbf{F}\in {\mathbb{R}}^{H\times W}$, the spectral coefficients $\widehat{\mathbf{F}}=\mathcal{F}\left\{\mathbf{F}\right\}$ are given by:

$${\widehat{F}}_{u,v}={\alpha }_u{\alpha }_v\sum _{x=0}^{H-1}\sum _{y=0}^{W-1}{F}_{x,y}cos\left[{\displaystyle \frac{\pi u(2x+1)}{2H}}\right]cos\left[{\displaystyle \frac{\pi v(2y+1)}{2W}}\right],$$

(17)

where $0\le u<H$, $0\le v<W$, and ${\alpha }_u,{\alpha }_v$ are normalization constants. In this basis, the Laplacian is diagonal:

$$\mathcal{F}{\left\{{\Delta }_h{U}_c\right\}}_{u,v}={\Lambda }_{u,v}{\widehat{U}}_{c,u,v},$$

with eigenvalues given analytically by:

$${\Lambda }_{u,v}=2\left(cos{\displaystyle \frac{\pi u}{H}}+cos{\displaystyle \frac{\pi v}{W}}\right)-4,{\Lambda }_{u,v}\in [-8,0].$$

(18)

The mode $(u,v)=(0,0)$ corresponds to the DC component with ${\Lambda }_{0,0}=0$.

Let ${\mathbf{R}}_c={\nabla }_h·{\mathbf{g}}_c-{\lambda }_c{\mathbf{T}}_c$ denote the discrete right-hand side. Applying $\mathcal{F}$ to the discrete screened Poisson equation decouples the system:

$$({\Lambda }_{u,v}-{\lambda }_c){\widehat{U}}_{c,u,v}={\widehat{R}}_{c,u,v},\forall (u,v).$$

(19)

The solution is obtained via a simple spectral filter:

$${\widehat{U}}_{c,u,v}={\displaystyle \frac{{\widehat{R}}_{c,u,v}}{{\Lambda }_{u,v}-{\lambda }_c}},{\mathbf{U}}_c={\mathcal{F}}^{-1}\left\{{\widehat{\mathbf{U}}}_c\right\}.$$

(20)

Stability and Regularization.

Since ${\lambda }_c>0$ and ${\Lambda }_{u,v}\le 0$, the denominator ${\mathcal{D}}_{u,v}:={\Lambda }_{u,v}-{\lambda }_c$ satisfies ${sup}_{u,v}{\mathcal{D}}_{u,v}=-{\lambda }_c<0$. In particular, the DC mode $(0,0)$ is well-defined, eliminating the global offset ambiguity typical of pure Poisson reconstruction. From an optimization perspective, the screening term ${\lambda }_c{\parallel U_c-T_c\parallel }_2^2$ acts as a Tikhonov regularizer, stabilizing the inversion of low-frequency modes and anchoring the solution to $T_c$.

5.3. Tiled SPN Implementation for Gigapixel WSIs

Direct application of the spectral solver to a full-resolution WSI is impractical due to memory constraints. We therefore adopt a non-overlapping domain decomposition [43] where the slide domain $\Omega$ is partitioned into tiles ${\left\{{\Omega }_k\right\}}_{k=1}^K$. The key requirement is that each local problem must approximate the restriction of the single global SPN functional.

Figure 2. Tiled SPN pipeline. In the first phase, we compute slide-level global prior statistics. In the second phase, we restrict these priors to each tile and solve the screened Poisson equation using a DCT-based spectral solver, assembling the normalized tiles into a coherent WSI.

Figure 2. Tiled SPN pipeline. In the first phase, we compute slide-level global prior statistics. In the second phase, we restrict these priors to each tile and solve the screened Poisson equation using a DCT-based spectral solver, assembling the normalized tiles into a coherent WSI.

Our implementation proceeds in two phases:

1. Global Prior Estimation. We first compute the global statistics of the source ($I_S$) and reference ($I_T$) slides. For each channel c, we calculate the global moments $({\mu }_S^{\mathrm{glob}},{\sigma }_S^{\mathrm{glob}})$ and $({\mu }_T,{\sigma }_T)$, define the chromatic anchor $T_c$ via the affine map (Eq. 6), and estimate the global contrast factor ${\alpha }_c$ (Eq. 7). These quantities are computed once at the slide level and provide coherent priors for all local reconstructions.

2. Local Spectral Reconstruction. We solve the screened Poisson equation independently on each tile. For a tile ${\Omega }_k$, we restrict the global priors ($T_{c,k}=T_c{|}_{{\Omega }_k}$ and $g_{c,k}=g_c^{\mathrm{glob}}{|}_{{\Omega }_k}$), assemble the local right-hand side ${\rho }_{c,k}={\nabla }_h·g_{c,k}-{\lambda }_c{T}_{c,k}$, and solve:

$$({\Delta }_h-{\lambda }_c)U_{c,k}={\rho }_{c,k}\mathrm{on}{\Omega }_k.$$

(21)

The homogeneous Neumann boundary conditions are implicitly enforced by the DCT-II/III basis. The normalized WSI is obtained by concatenating the tile solutions, ${\left({\tilde{I}}_S\right)}_c{|}_{{\Omega }_k}=U_{c,k}$. By Theorem 4, interface discrepancies are exponentially localized near the boundaries with characteristic width $\mathcal{O}\left({\lambda }_c^{-1/2}\right)$.

Algorithm 1 Tiled Screened Poisson Normalization (SPN)

Require:  Source $I_S$, Reference $I_T$, parameters $\tau ,{\lambda }_c$, tiling ${\left\{{\Omega }_k\right\}}_{k=1}^K$ Ensure:  Normalized WSI ${\tilde{I}}_S$ 1: Phase 1: Global Statistics 2: for each channel $c\in \{L^{*},a^{*},b^{*}\}$ do 3:     Compute global source moments on $I_S$ and target moments on $I_T$ 4:     Define global chromatic anchor $T_c$ via Eq. (6) 5:     Compute global contrast factor ${\alpha }_c$ via Eq. (7) 6: end for 7: Phase 2: Tile-wise Spectral Reconstruction 8: for each tile ${\Omega }_k$ in parallel do 9:     Extract source patch $P_k\leftarrow I_S{|}_{{\Omega }_k}$ 10:     Compute structure weight $w_k\left(\mathbf{x}\right)\leftarrow exp(-\parallel \nabla {\left(P_k\right)}_{L^{*}}\left(\mathbf{x}\right){\parallel }^2/{\tau }^2)$ 11:     for each channel $c\in \{L^{*},a^{*},b^{*}\}$ do 12:         Restrict anchor: $T_{c,k}\leftarrow T_c{|}_{{\Omega }_k}$ 13:         Form guidance field: 14:         ${\mathbf{g}}_{c,k}\leftarrow (1-w_k){\alpha }_c\nabla {\left(P_k\right)}_c+w_k\nabla T_{c,k}$ 15:         Assemble RHS: ${\rho }_{c,k}\leftarrow {\nabla }_h·{\mathbf{g}}_{c,k}-{\lambda }_c{T}_{c,k}$ 16:         Solve $({\Delta }_h-{\lambda }_c)U_{c,k}={\rho }_{c,k}$ via DCT spectral filter (20) 17:         Update ${\tilde{I}}_S{|}_{{\Omega }_k}\leftarrow U_{c,k}$ 18:     end for 19: end for 20: return${\tilde{I}}_S$

The boundary behavior of this tiled scheme is governed by the screened Poisson operator analyzed in Section 4. The screening term ${\lambda }_c>0$ introduces a characteristic length scale ${\ell }_c\approx {\lambda }_c^{-1/2}$, ensuring exponential decay of flux mismatches away from tile interfaces. Consequently, the assembled reconstruction is structurally and chromatically consistent with the global solution, with seams reduced to negligible boundary effects.

5.4. Computational Complexity

The computational cost is dominated by the 2D DCT spectral solves. For a single tile with $N_k=H_k{W}_k$ pixels:

• Assembly: Computing gradients, weights, and the right-hand side ${\rho }_{c,k}$ scales linearly, i.e., $\mathcal{O}\left(N_k\right)$.
• Spectral Solve: The 2D DCT is implemented via separable 1D fast cosine transforms, with a cost of $\mathcal{O}(N_{k}log{N}_k)$.

The total cost for a WSI with $N={\sum }_k{N}_k$ pixels is $\mathcal{O}(NlogN)$, with memory complexity $\mathcal{O}\left(N_k\right)$ per tile. This approach allows for trivial parallelization over tiles. In comparison, sparse direct solvers typically scale as $\mathcal{O}\left(N^{1.5}\right)$ for 2D problems, while multigrid methods, though $\mathcal{O}\left(N\right)$, require complex data structures. The tiled DCT implementation thus provides an effective balance between asymptotic efficiency and implementation simplicity for gigapixel images.

6. Experiments and Results

6.1. Experimental Setup

Datasets 

To evaluate the effectiveness and robustness of the proposed algorithm, we conducted comparative experiments on multiple publicly available and clinical datasets exhibiting diverse staining conditions and imaging variations. In total, the evaluation encompassed 88 patients and more than 300 whole-slide histopathological images (WSIs).

• Dataset 1: MITOS-ATYPIA-14 [35], released as part of the MITOS-ATYPIA Grand Challenge, was annotated by the Pathology Department of Pitié-Salpêtrière Hospital in Paris. It consists of 14 pairs of breast carcinoma slides, each scanned using two distinct digital pathology scanners—Aperio ScanScope XT and Hamamatsu NanoZoomer—allowing controlled analysis of inter-scanner color variation.
• Dataset 2: Multi-Scanner Squamous Cell Carcinoma (SCC) dataset [36] contains 44 samples of canine cutaneous squamous cell carcinoma digitized using 5 different scanning systems, producing a total of 220 WSIs.
• Dataset 3: HE-Staining Variation (HEV) dataset [37] provided by Heidelberg University. The HEV dataset comprises follicular thyroid carcinoma slides prepared under nine distinct staining protocols: standard H&E, over-stained H&E (longHE), under-stained H&E (shortHE), hematoxylin-only (onlyH), eosin-only (onlyE), over-stained hematoxylin (longH), over-stained eosin (longE), under-stained hematoxylin (shortH), and under-stained eosin (shortE). These variations facilitate the evaluation of color normalization under chemically induced staining inconsistencies.
• Dataset 4: SUSY-BF-10 clinical dataset (local collection). To address the temporal color degradation that public datasets do not capture, we constructed a clinical dataset from Sun Yat-sen Memorial Hospital, Sun Yat-sen University. The SUSY-BF-10 dataset includes 25 long-term archived breast carcinoma slides, each preserved for over ten years and rescanned to analyze fading and chromatic shifts over time.

A summary of dataset characteristics and representative visualizations is provided in Table 1 and Figure 3.

Implementation Details 

In surveying the literature on stain normalization, we observe that most existing methods are evaluated primarily on pre-extracted WSI tiles and are often positioned as a preprocessing step for downstream tasks such as classification or segmentation. In contrast, the ultimate goal of digital pathology in the clinical setting is to provide pathologists with visually consistent slides, especially in cross-site, multi-center collaborations and remote diagnostic workflows. Under these scenarios, a robust and interpretable color normalization scheme at the whole-slide level is of central importance.

Motivated by this observation, we design our experiments to assess performance from two complementary perspectives:

1.

Full-resolution WSI normalization. We evaluate the behavior of each method when applied directly to entire, full-size histopathology slides.

2.

Tile-based reconstruction. We further consider the practically common pipeline in which a large WSI is first partitioned into tiles, normalized tile-wise, and then reassembled into a full-resolution slide. This setting is particularly relevant for memory-constrained or streaming-based deployments.

In addition to the four real-world datasets used in our study, we also construct a synthetic test image (synthetic_image) with controlled color and structural patterns. This synthetic benchmark provides an intuitive and stress-test-style assessment of robustness and generalization across a wide range of stain variations.

All experiments are conducted on a server equipped with an Intel(R) Xeon(R) Platinum 8352V CPU @ 2.10 GHz, an NVIDIA RTX 4090 GPU with 24 GB of VRAM, and 512 GB of system memory. The software stack is based on Python 3.8 and PyTorch 2.5.

Baselines 

We compare the proposed method against three representative baselines:

• Reinhard. The classical global Reinhard color transfer algorithm [6] serves as our primary reference for global-statistics–based normalization. It matches low-order statistics between source and target images in a decorrelated color space and remains widely used in digital pathology.
• GCTI. We include a geometry-aware method with local alignment in the stain vector space (GCTI) [9], which explicitly exploits the underlying geometric structure of color distributions to compensate for scanner- and protocol-induced variations in microscopy images within a unified framework.
• CycleGAN. As a deep learning–based baseline, we employ a CycleGAN-style stain transfer model [10]. This method learns a non-linear mapping between source and target stain domains from unpaired data and represents a strong data-driven alternative to hand-crafted normalization schemes.

These baselines jointly cover global-statistics approaches, geometry-aware model-based methods, and modern generative adversarial networks, providing a comprehensive context for evaluating the proposed gradient-domain diffusion normalization.

6.2. Evaluation Metrics

In practice, there is still no widely accepted, unified standard for quantitative evaluation in the stain normalization task. Existing work typically adopts a mixture of structural similarity measures, color-space or gamut-based criteria, and performance gains on downstream tasks (e.g., classification or segmentation). From a clinical perspective, however, the ideal outcome is clear: the normalized image should (i) preserve the structural and morphological content of the source image and (ii) match the color space and distribution of the target image as closely as possible. Achieving this would allow pathology images acquired at different hospitals, regions, and from heterogeneous devices to be mapped to a common color standard, thereby reducing color-induced variability, improving diagnostic efficiency for pathologists, and ultimately promoting more standardized diagnostic criteria.

To better assess the performance of different algorithms under this perspective, we evaluate: (1) the structural similarity between the source image and the normalized result, and (2) the chromatic alignment between the normalized result and the target image. Concretely, we use the Structural Similarity Index (SSIM) [44] as a morphology-oriented metric between the source $I_S$ and the normalized image ${\tilde{I}}_S$, and the 1st-Wasserstein distance [45] between the color distributions of ${\tilde{I}}_S$ and the target $I_T$. Considering the entanglement of structure and color in the RGB space, all metrics are computed in the CIELab color space. The specific definitions are as follows.

Structural Similarity Index (SSIM). To quantify the preservation of tissue morphology and local texture, we employ the SSIM metric. Let $\mathbf{x}$ and $\mathbf{y}$ denote local window patches extracted from the luminance channels of the source image $I_S$ and the normalized image ${\tilde{I}}_S$, respectively. The SSIM is defined as:

$$\mathrm{SSIM}(\mathbf{x},\mathbf{y})={\displaystyle \frac{(2{\mu }_x{\mu }_y+C_1)(2{\sigma }_{xy}+C_2)}{({\mu }_x^2+{\mu }_y^2+C_1)({\sigma }_x^2+{\sigma }_y^2+C_2)}},$$

(22)

where ${\mu }_x,{\mu }_y$ represent the local means, ${\sigma }_x^2,{\sigma }_y^2$ denote the local variances, and ${\sigma }_{xy}$ is the cross-covariance between the two patches. The constants $C_1={\left(k_{1}L\right)}^2$ and $C_2={\left(k_{2}L\right)}^2$ are included to ensure numerical stability, where L is the dynamic range of pixel values (e.g., $L=255$ for 8-bit images). The final SSIM score is computed as the mean over all sliding windows in the image domain.

Wasserstein Distance (W-D). To evaluate chromatic alignment, we calculate the 1st-Wasserstein distance (also known as Earth Mover’s Distance) between the color distributions of the target image $I_T$ and the result ${\tilde{I}}_S$. For a specific color channel $c\in \{a^{*},b^{*}\}$, let $P_c$ and $Q_c$ represent the probability distributions of pixel intensities in $I_T$ and ${\tilde{I}}_S$, respectively. The Wasserstein distance $W_1(P_c,Q_c)$ is defined as the minimal cost to transform one distribution into the other:

$$W_1(P_c,Q_c)=\underset{\gamma \in \Pi (P_c,Q_c)}{inf}{\mathbb{E}}_{(u,v)\sim \gamma }\left[\right|u-v\left|\right],$$

(23)

where $\Pi (P_c,Q_c)$ denotes the set of all joint distributions $\gamma (u,v)$ whose marginals are $P_c$ and $Q_c$. Since we operate on 1D marginal distributions of color channels, Eq. (23) has a closed-form solution expressed via the cumulative distribution functions (CDFs), $F_T\left(z\right)$ and $F_R\left(z\right)$:

$$W_1(P_c,Q_c)={\int }_{-\infty }^{+\infty }|F_T\left(z\right)-F_R\left(z\right)|dz.$$

(24)

A lower $W_1$ score indicates superior alignment of the statistical color profile, without assuming Gaussianity of the underlying distributions.

6.3. Results and Visual Analysis

We complement the quantitative evaluation with a detailed visual comparison across four representative datasets and a synthetic benchmark. This analysis focuses on both full-size whole-slide views and tiled reconstructions.

Full size baseline test

As illustrated in the first two rows of Figure 4, both Reinhard and SPN achieve faithful chromatic adaptation. This performance is consistent with their explicit statistical matching in the CIELab space, where the $a^{*}$ and $b^{*}$ channels isolate chromaticity from luminance. In contrast, the third and fourth rows reveal the limited generalization of CycleGAN on out-of-distribution data: when applied to non-MITOS samples, the model introduces a characteristic “deep red” bias inherent to its training set rather than correctly adapting to the target domain. Furthermore, compared to the global Reinhard method, SPN and GCTI preserve high-frequency morphological details, exhibiting significantly sharper cellular boundaries at high magnification.

SCC dataset: full-size vs. tiled normalization.

The SCC dataset contains slides acquired from five different scanners with codes cs2, gt450, nz20, nz210, and p100. We use cs2 as the source device and treat the remaining four devices as distinct targets. Figure 5 compares color normalization on both full-size images and tiled images. The full-size resolution is $2000\times 1600$, and the tiled setting uniformly partitions each image into a $5\times 4$ grid of non-overlapping patches.

As shown in Figure 5, tile-wise normalization via Reinhard and GCTI introduces distinct boundary discontinuities arising from inconsistent local statistics. While SPN effectively eliminates these seams through gradient preservation, CycleGAN displays a diminished checkerboard effect relative to the statistical methods. This partial suppression is attributed to the shared generator, which enforces a consistent color transform across tiles. Yet, this same mechanism limits adaptability, causing CycleGAN to hallucinate the MITOS style onto the SCC samples rather than respecting the target domain’s specific appearance.

HEV dataset: Multi-stained pathology images.

The HEV dataset encompasses nine distinct staining configurations, offering the richest color diversity among the datasets analyzed. We utilize standard HE as the source domain, treating the remaining eight protocols as targets. For brevity, and because full-slide comparisons show limited variation, Figure 6 focuses exclusively on the tiled results.

The CycleGAN results in Figure 6 mirror those in Figure 5: blocking artifacts are diminished relative to Reinhard and GCTI due to the global constraints of the generator. More significantly, this experiment highlights the robustness limitations of the GCTI framework. GCTI employs a local wedge-finding algorithm that projects pixel intensities onto a subspace defined by SVD, inferring stain vectors via angular percentiles (e.g., $0.1$th and $99.9$th). This mechanism is susceptible to failure in extreme scenarios: (1) Illumination correction singularities, where near-zero background estimates lead to division-by-zero errors and intensity saturation; and (2) Thresholding instability in uniform regions, where the lack of color variance renders percentile-based thresholds arbitrary. Figure 6 clearly demonstrates the former, where instability in the background estimation step produces extensive artifacts characterized by signal saturation and loss of texture.

Evaluation on Synthetic Image

To rigorously assess algorithmic stability under limiting conditions, we introduce a synthetic image characterized by a smooth background gradient and a central geometric object (Figure 7).

The first row of Figure 7 presents the full-resolution normalization results.SPN and Reinhard achieve the most faithful reconstruction of the background gradient, demonstrating robust structural consistency across both low-frequency trends and high-frequency edges. In contrast, GCTI exhibits a catastrophic failure mode driven by algorithmic instability: the central large circle is misclassified as background, causing the denominator in the illumination correction step to approach zero. This singularity leads to numerical explosion, resulting in a saturated, near-uniform output in the background region.

Additionally, the magnified regions in the second row reveal that CycleGAN produces spurious ringing artifacts around the circular structures, a phenomenon inconsistent with the behavior of the variational and statistical methods. This artifact stems from the generative network’s reliance on learned upsampling filters (transposed convolutions) which, unlike rigorous interpolation schemes, can hallucinate high-frequency noise to satisfy the discriminator. Such artifacts underscore the critical deficit of deep generative models regarding interpretability: unlike the proposed framework, which is constrained by a physical diffusion process, the GAN introduces opaque, architectural artifacts that compromise diagnostic confidence. The tiled results (Row 3) further expose the stability limits of the baseline methods, particularly regarding well-posedness on disjoint subdomains. First, statistical methods break down on homogeneous patches. Reinhard normalization relies on matching first- and second-order moments; when a target tile consists of a pure background (zero variance), the scaling factors become undefined, rendering the transfer function ill-posed. Similarly, GCTI fails in these homogeneous regions because the lack of spectral diversity renders the SVD-based stain estimation rank-deficient, leading to numerical breakdown.

Conversely, SPN maintains both global gradient continuity and local geometric fidelity across all tiles. This robustness confirms that the proposed variational framework effectively regularizes the normalization problem, preventing the degenerate behavior observed in purely statistical approaches when local data is sparse.

6.4. Quantitative Results and Discussion

Table 2 presents a comprehensive quantitative comparison of our proposed SPN against three state-of-the-art methods: CycleGAN, Reinhard, and GCTI. We utilize the Structural Similarity Index Measure (SSIM) to evaluate structure preservation (↑) and the Wasserstein Distance (W-D) to assess stain distribution matching (↓).

It is worth noting that the classic Reinhard and GCTI algorithms remain strong baselines, particularly on relatively simple and highly structured datasets such as SYSU-BF-10 and MITOS-14. In these cases, Reinhard achieves competitive or even state-of-the-art results (e.g., SSIM of 0.981 on SYSU-BF-10). This is likely because these datasets exhibit unimodal color distributions where simple global statistical alignment (mean and standard deviation matching) is sufficient. However, its performance degrades significantly on complex datasets like HEV (W-D often $>4.0$), where stain variations are non-linear and multi-modal.

As previously discussed, CycleGAN utilizes fixed weights trained on MITOS. While it yields consistent structural outputs (omitted entries in Table 2), its inability to adapt to the specific target reference leads to large domain discrepancies, evidenced by high W-D scores (e.g., 25.624 on SCC).

Our method demonstrates robust performance across diverse scenarios, with its advantages becoming most pronounced on the synthetic_image. This dataset serves as a rigorous benchmark for “general stain normalization” as it provides paired ground truth. On this set, SPN achieves the best results in both metrics (SSIM: 0.950, W-D: 3.706). This confirms that our gradient-guided diffusion mechanism effectively learns the intrinsic mapping between domains without overfitting to specific dataset artifacts, demonstrating superior generalization compared to baselines.

7. Conclusions

We presented Screened Poisson Normalization (SPN), a physics-informed framework that reformulates stain normalization as a screened Poisson reconstruction problem. By modeling the staining process as a reaction–diffusion system, SPN explicitly bridges statistical moment matching with physically interpretable gradient-domain regularization.

Our theoretical and experimental analysis yields three key conclusions:

1.

Structural Fidelity and Scalability: Operating in the gradient domain allows SPN to preserve diagnostic morphological details (e.g., nuclear boundaries) often degraded by pure statistical alignment. Furthermore, the exponential decay of the screened Green’s function enables rigorous non-overlapping domain decomposition, permitting seamless gigapixel processing without tiling artifacts.

2.

Determinism vs. Data-Dependence: In contrast to deep learning approaches (e.g., CycleGAN) which are prone to domain bias and generative hallucinations on out-of-distribution data, SPN is training-free and deterministic. This ensures that all output structures are causally linked to the input, providing the explainability and reproducibility essential for clinical diagnostics.

3.

Robustness to Heterogeneity: SPN demonstrates consistent performance across diverse scanners, staining protocols, and archived slides. The global anchoring (screening) term effectively stabilizes the solution in tissue-sparse or background-dominated regions, overcoming the numerical brittleness observed in local geometry-based methods.

In summary, SPN offers a transparent, mathematically well-posed alternative to "black-box" style transfer, establishing a stable foundation for computational pathology in heterogeneous clinical settings. Future work will focus on parallel optimization for high-throughput deployment and extending the gradient-domain formulation to multi-stain separation and volumetric reconstruction.