A General Probability Density Framework for Local Histopolation and Weighted Function Reconstruction from Mesh Line Integrals

Abstract

In this paper, we study the reconstruction of a bivariate function from weighted integrals along the edges of a triangular mesh, a problem of central importance in tomography, computer vision, and numerical approximation. Our approach relies on local histopolation methods defined through unisolvent triples, where the edge weights are induced by suitable probability densities. In particular, we introduce two new two-parameter families of generalized truncated normal distributions, which extend classical exponential-type laws and provide additional flexibility in capturing local features of the target function. These distributions give rise to new quadratic reconstruction operators that generalize the standard linear histopolation scheme, while retaining its simplicity and locality. We establish their theoretical foundations, proving unisolvency and deriving explicit basis functions, and we demonstrate their improved accuracy through extensive numerical tests. Moreover, we design an algorithm for the optimal selection of the distribution parameters, ensuring robustness and adaptivity of the reconstruction. Finally, we show that the proposed framework naturally extends to any bivariate function whose restriction to the edges defines a valid probability density, thus highlighting its generality and broad applicability.

1 Introduction

The classical concept of polynomial interpolation can be extended beyond pointwise evaluations to situations where the available information consists of general linear functionals Rivlin. When these data take the form of integrals over geometric domains, the problem is referred to as histopolation Robidoux, and the corresponding techniques are known as histopolation methods.

Histopolation provides a general framework for reconstructing functions from integral data, a problem of central importance in many areas of science and engineering. In practical applications, direct pointwise samples of the target function are often unavailable; instead, only integral measurements are accessible. This situation arises prominently in computed tomography and medical imaging, where measurements typically consist of line integrals kak2001principles; natterer2001mathematics; palamodov2016reconstruction. Beyond tomography, histopolation has found applications in computer vision, signal and image processing Bosner:2020:AOC, curve and spline approximation Fischer:2005:MPR; Fischer:2007:CSP; Siewer:2008:HIS; Hallik:2017:QLR, fractal functions Barnsley:2023:HFF, and in the conservation of physical quantities in numerical simulations HiemstraJCP. Recently, histopolation methods have received increasing attention BruniErbFekete; bruni2024polynomial, including global polynomial histopolation–regression in one bruni2025polynomial and several dimensions bruno2025bivariate, as well as weighted local polynomial histopolation schemes in the bivariate setting nudo1; nudo2; dell2025reconstructing; dell2025truncated; milov1.

A standard strategy for these reconstruction problems is to partition the computational domain into subdomains (a mesh), approximate locally on each element, and assemble a global approximation from the local contributions. Such schemes are typically conforming or nonconforming, depending on whether the assembled approximation is continuous across element interfaces.

A local approximation method can be rigorously defined through the triple Guessab:2022:SAB

$$\mathcal{M}_{d}=\left(S_{d},\mathbb{F}_{S_{d}},\Sigma_{S_{d}}\right),$$

where

• $S_{d}$ is a polytope in $\mathbb{R}^{d}$, $d\geq 1$;

• $\mathbb{F}_{S_{d}}$ is an $n$-dimensional space of trial functions on $S_{d}$;

• $\Sigma_{S_{d}}=\left\{\mathcal{L}_{j}:j=1,\dots,n\right\}$ is a set of linearly independent linear functionals, called the degrees of freedom.

The triple $\mathcal{M}_{d}$ is said to be unisolvent if, whenever $f\in\mathbb{F}_{S_{d}}$ satisfies

$$\mathcal{L}_{j}(f)=0,\quad j=1,\dots,n,$$

then necessarily $f=0$. A set of functions

$$B=\left\{\varphi_{i}\,:\,i=1,\dots,n\right\}$$

such that

$$\operatorname{span}\left\{\varphi_{i}\,:\,i=1,\dots,n\right\}=S_{d}$$

and

$$\mathcal{L}_{j}\left(\varphi_{i}\right)=\delta_{ij}$$

is said to form the basis functions associated with the triple $\mathcal{M}_{d}$. If the degrees of freedom are integral functionals, the scheme is also referred to as a local histopolation method.

The simplest example is the linear nonconforming histopolation method on triangles, where edge integrals act as degrees of freedom and piecewise linear functions form the trial space. While attractive for its simplicity and efficiency, this approach has important limitations: its piecewise linear structure fails to reproduce oscillatory behaviors or sharp gradients, and higher accuracy can only be achieved through substantial mesh refinement.

In this work we propose an enrichment strategy that enhances local histopolation schemes by introducing quadratic trial spaces together with weighted degrees of freedom. The key idea is to employ probability densities on the edges of the mesh, which act as flexible weights and introduce tunable parameters. In particular, we construct two new two-parameter families of generalized truncated normal distributions, which extend classical exponential-type laws and include special or limiting cases (such as truncated normal and beta-type densities). These families naturally lead to enriched quadratic operators that significantly improve reconstruction accuracy while retaining the locality and simplicity of the original method.

From a theoretical point of view, we prove unisolvency of the enriched triples and derive explicit closed-form basis functions. From a practical perspective, we validate the proposed operators through extensive numerical experiments, which demonstrate their effectiveness in reconstructing both smooth and non-smooth functions. Moreover, as shown in Section 3, the framework is not limited to the two families introduced here: the same reasoning applies to any bivariate function whose restriction to the edges defines a valid probability density, thus highlighting the generality and robustness of the proposed approach.

The paper is organized as follows. In Section 2, we introduce two families of generalized truncated normal distributions, analyze their main properties, and define the corresponding families of quadratic histopolation operators. Moreover, we describe an algorithm for the optimal selection of the distribution parameters, which further enhances the robustness and adaptivity of the reconstruction. In Section 3, we develop the general framework for enriched quadratic histopolation operators and establish their theoretical foundations. In Section 4, we present numerical experiments that confirm the accuracy and effectiveness of the proposed methods. Finally, in Section 5, we conclude with a summary of the contributions and discuss possible directions for future research, including three-dimensional extensions and applications to imaging and approximation problems.

2 Two-parameter families of generalized truncated normal distributions

2.1 Problem statement

The problem we want to address is the following: given only integral measurements of an unknown bivariate function along the edges of a triangular mesh, reconstruct an accurate approximation of the function inside each triangle. In the classical local histopolation method, these measurements are taken as unweighted edge integrals and the reconstruction space is piecewise linear. While simple and efficient, this approach lacks the ability to capture curved or oscillatory features unless the mesh is strongly refined.

Our objective is to overcome this limitation by enriching the standard framework. We introduce probability density functions as edge weights in the integral data, so that each measurement carries additional information about the local behavior of the function. This enrichment is feasible under the assumption that additional integral data are available, beyond the classical unweighted edge averages. In this setting, the reconstruction naturally leads to quadratic operators that preserve locality and simplicity while significantly improving approximation accuracy. The role of this section is to present two flexible two-parameter distribution families that serve as the building blocks of the proposed weighted histopolation scheme.

2.2 Definition of the distribution families

In this section, we introduce two distinct two-parameter families of univariate distributions, each generated by suitable bivariate weight functions defined on a triangle of a given mesh. These families include, as limiting or special cases, several well-known probability laws. They provide the building blocks for a broad class of new piecewise quadratic reconstruction operators, which in turn form the basis of a nonconforming histopolation method locally determined by a unisolvent triple. The goal is to reconstruct a function using information derived from weighted integrals over the edges of a triangular mesh. To this end, we formulate the probability density functions (hereafter, pdfs) associated with the edges of a nondegenerate triangle $T$ with vertices $\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}$ and barycentric coordinates $\lambda_{1},\lambda_{2},\lambda_{3}$. These coordinates are the unique affine functions on $T$, namely

$$\lambda_{i}\left(t\boldsymbol{x}+(1-t)\boldsymbol{y}\right)=t\lambda_{i}(\boldsymbol{x})+(1-t)\lambda_{i}(\boldsymbol{y}),\quad\boldsymbol{x},\boldsymbol{y}\in T,\quad t\in[0,1],$$ (1)

which satisfy the partition of unity property,

$$\sum_{i=1}^{3}\lambda_{i}(\boldsymbol{x})=1,\quad\boldsymbol{x}\in T,$$

and the Kronecker delta property at the vertices,

$$\lambda_{i}\left(\boldsymbol{v}_{j}\right)=\delta_{ij}=\begin{cases}1,&i=j,\\[4.0pt] 0,&i\neq j,\end{cases}\quad i,j=1,2,3.$$

As a consequence, we have

$$\lambda_{i}(\boldsymbol{x})=0,\quad\boldsymbol{x}\in s_{i},\quad i=1,2,3.$$ (2)

For the forthcoming analysis, we recall the lower incomplete gamma function, which will play a central role. It is defined for $s>0$ and $z>0$ by

$$\gamma(s,z)=\int_{0}^{z}t^{s-1}e^{-t}dt.$$ (3)

Among its numerous properties, one that we shall exploit is the limiting relation

$$\lim_{z\to 0}\frac{\gamma(s,z)}{z^{s}}=\frac{1}{s}.$$ (4)

Indeed, since

$$e^{-z}\leq e^{-t}\leq 1,\quad 0\leq t\leq z,$$

we obtain

$$t^{s-1}e^{-z}\leq t^{s-1}e^{-t}\leq t^{s-1}.$$

Integrating over $[0,z]$ with respect to $t$ and dividing by $z^{s}$, we obtain the inequalities

$$\frac{e^{-z}}{s}\leq\frac{\gamma(s,z)}{z^{s}}\leq\frac{1}{s}.$$

Taking the limit as $z\to 0$ establishes (4). This relation will be crucial for analyzing the convergence of our densities when one of their parameters tends to infinity. For convenience, we introduce the modified incomplete gamma function

$$\gamma^{\mathrm{mod}}(s,z):=\frac{\gamma(s,z)}{z^{s}}=\frac{1}{z^{s}}\int_{0}^{z}t^{s-1}e^{-t}dt,$$ (5)

so that (4) takes the simpler form

$$\lim_{z\to 0}\gamma^{\mathrm{mod}}(s,z)=\frac{1}{s}.$$ (6)

If a numerical routine is available to evaluate $\gamma(s,z)$, then $\gamma^{\mathrm{mod}}(s,z)$ can be computed directly by scaling. For example, in Matlab one may use gammainc(z,s) gamma(s), while in Mathematica the corresponding function is Gamma[s,0,z].

2.3 First family of distributions

The first family is defined through the following bivariate function with shape parameter $\mu\geq 1$ and scale parameter $\sigma>0$ by

$$K_{\sigma,\mu}(\boldsymbol{x})=a_{\sigma,\mu}\left(H^{2}(\boldsymbol{x})\right)^{\mu-1}e^{-\frac{1}{2}\left(\frac{H(\boldsymbol{x})}{\sigma^{2}}\right)^{\mu}},\quad\boldsymbol{x}\in T,$$ (7)

where $H$ is expressed in terms of barycentric coordinates as

$$H(\boldsymbol{x})=2\sum_{i=1}^{3}\lambda_{i}^{2}(\boldsymbol{x})-1,\quad\boldsymbol{x}\in T,$$ (8)

and $a_{\sigma,\mu}$ is a normalization factor ensuring that integrals involving the associated generalized truncated normal distribution are properly scaled. It is worth noting that $K_{\sigma,\mu}$ is a weight function, as it is nonnegative throughout the triangle $T$. Furthermore, when restricted to any edge of $T$, $H$ reduces to a quadratic form, and the induced univariate functions not only serve as weights but also define valid probability density functions on the edges. An advantage of the general form (7) is that it allows us to tune the free parameters $\sigma$ and $\mu$ when reconstructing a function with specific features. In particular, different behaviors can be produced by varying the exponential term in (7).

In the following, we parametrize the edge $s_{j}$ as

$$\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2},\quad t\in[-1,1].$$ (9)

Then, by (1) and (2), together with (8), we obtain

$$H\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)=t^{2}.$$ (10)

We now define the first class of univariate pdfs associated with the weight function $K_{\sigma,\mu}$. It follows directly from the definition of $K_{\sigma,\mu}$ that on the edge $s_{i}$ of $T$ one has

$$k_{\sigma,\mu}(t):=K_{\sigma,\mu}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)=a_{\sigma,\mu}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}},\quad t\in\left[-1,1\right].$$ (11)

An important observation is that $k_{\sigma,\mu}(t)$ is independent of the choice of edge of the triangle $T$. We shall henceforth work with the normalized probability density functions (pdfs) with unit zeroth moment, which we denote with a tilde. Thus, in order to obtain a normalized pdf, the coefficient $a_{\sigma,\mu}$ in the bivariate weight function $K_{\sigma,\mu}(\boldsymbol{x})$ defined in (7) is introduced so that

$$\int_{-1}^{1}k_{\sigma,\mu}(t)dt=1.$$

As shown in Lemma 2, for any $\mu\geq 1$ and $\sigma>0$, the normalization constant $a_{\sigma,\mu}$ is given by

$$a_{\sigma,\mu}=\frac{\mu}{2^{\frac{4\mu-3}{2\mu}}\sigma^{4\mu-3}\gamma\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)},$$

which can be written more compactly in terms of the modified incomplete gamma function as

$$a_{\sigma,\mu}=\frac{\mu}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}.$$ (12)

Hence, for any $\sigma>0$ and $\mu\geq 1$, substituting (12) into (11), it follows that the normalized pdfs on $[-1,1]$ are given by

$$\widetilde{k}_{\sigma,\mu}(t)=\frac{\mu}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}},\quad t\in\left[-1,1\right].$$ (13)

As a first observation, this family extends the class of exponential-type distributions. Indeed, when $\mu=1$, equation (13) reduces to the well-known doubly truncated normal distribution on the interval $[-1,1]$, here written in terms of $\gamma^{\mathrm{mod}}$; see Jawitz2004. Indeed, by straightforward calculation, and using the notation (6), this distribution can be written as

$$\widetilde{k}_{\sigma,1}(t)=\frac{1}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{2}}\right)}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)},\quad t\in\left[-1,1\right].$$

Recall that this latter distribution has no shape parameter, and thus cannot model all phenomena. Nevertheless, it has been extensively studied and plays a central role in probability theory, see, e.g., Xinyi2012functions and the references therein. Therefore, our pdfs constitute a new two-parameter family of exponential type. In Figure 1, we plot several instances of $\widetilde{k}_{\sigma,\mu}$, illustrating the strong sensitivity of the distribution to both the shape parameter $\mu$ and the scale parameter $\sigma$. This highlights the increased flexibility and importance of the generalized family.

Figure 1: Normalized probability density functions $\widetilde{k}_{\sigma,\mu}$ on $[-1,1]$. Left: effect of the scale parameter, with fixed $\mu=2$ and different values of $\sigma$. Right: effect of the shape parameter, with fixed $\sigma=0.50$ and varying values of $\mu$, illustrating the increased sharpness and concentration around $t=0$. These plots emphasize the flexibility provided by the two parameters.

Remark 1

The pdfs $\widetilde{k}_{\sigma,\mu}$ also include a beta-type distribution as a limiting case. Indeed, using (6) with

$$s=\frac{4\mu-3}{2\mu},\quad z=\frac{1}{2\sigma^{2\mu}},$$

we obtain

$$\lim_{\sigma\to\infty}\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)=\frac{2\mu}{4\mu-3}.$$ (14)

Therefore, for any $t\in[-1,1]$, the limiting distribution is

$$\widetilde{k}_{\infty,\mu}(t):=\lim_{\sigma\to\infty}\widetilde{k}_{\sigma,\mu}(t)=\frac{4\mu-3}{2}\left(t^{2}\right)^{2\mu-2}.$$ (15)

Thus, $\widetilde{k}_{\sigma,\mu}$ converges to a beta-type distribution as $\sigma\to\infty$.

In order to present our local reconstruction operator, which is based on the pdfs $\widetilde{k}_{\sigma,\mu}$, we now introduce some additional notations and recall the basic setting for clarity of exposition. The operator is built locally on each nondegenerate triangle $T$ of a given mesh. For notational convenience, we set

$$\boldsymbol{v}_{4}=\boldsymbol{v}_{1},\quad\boldsymbol{v}_{5}=\boldsymbol{v}_{2},$$

as well as

$$\lambda_{4}=\lambda_{1},\quad\lambda_{5}=\lambda_{2}.$$

Then, using the parametrization of the side $s_{j}$ given in (9), we introduce the linear functionals

$$\mathcal{I}_{j}^{\mathrm{CH}}(f)=\frac{1}{\left\lvert s_{j}\right\rvert}\int_{s_{j}}f(\boldsymbol{x})d\boldsymbol{x}=\frac{1}{2}\int_{-1}^{1}f\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)d\boldsymbol{x},\quad j=1,2,3,$$ (16)

and denote

$$\Sigma^{\mathrm{CH}}=\left\{\mathcal{I}_{j}^{\mathrm{CH}}\,:\,j=1,2,3\right\}.$$

Thus, the linear nonconforming histopolation method can be locally represented by the triple

$$\mathcal{CH}=\left(T,\mathbb{P}_{1}(T),\Sigma^{\mathrm{CH}}\right),$$ (17)

where $\mathbb{P}_{1}(T)$ is the space of linear polynomials on $T$, and its associated basis functions are expressed in barycentric coordinates as

$$\varphi_{i}^{\mathrm{CH}}=1-2\lambda_{i},\quad i=1,2,3,$$

as can be directly verified. The reconstruction operator associated to the triple (17) is

$$\begin{array}[]{rcl}{\pi}^{{\mathrm{CH}}}:C(T)&\rightarrow&{\mathbb{P}}_{1}(T)\\ f&\mapsto&\displaystyle{\sum_{j=1}^{3}}\mathcal{I}^{\mathrm{CH}}_{j}\left(f\right)\varphi^{\mathrm{CH}}_{j}.\end{array}$$ (18)

To describe our first new two-parameter family of nonconforming histopolation operators, we introduce a class of enriched linear functionals, defined by

$${\mathcal{I}}_{j}^{\mathrm{enr}}(f) = \int_{-1}^{1}f\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{k}_{\sigma,\mu}(t)dt,\quad j=1,2,3,$$ (19)

$${\mathcal{L}}_{j}^{\mathrm{enr}}(f) = \int_{-1}^{1}{p}_{\sigma,\mu}(t)f\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{k}_{\sigma,\mu}(t)dt,\quad j=1,2,3,$$ (20)

where

$${p}_{\sigma,\mu}(t)=t^{2}-\frac{\gamma^{\mathrm{mod}}\left(\frac{4\mu-1}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}.$$ (21)

Since $\widetilde{k}_{\sigma,\mu}$ is a probability density, the functionals $\mathcal{I}_{j}^{\mathrm{enr}}$ and $\mathcal{L}_{j}^{\mathrm{enr}}$ can be interpreted as weighted averages of $f$ and $p_{\sigma,\mu}f$ along the side $s_{j}$ of $T$. The rationale behind the choice of $p_{\sigma,\mu}$, which will become clear in Section 3, lies in its central role in the subsequent analysis. Using (10) and the parametrization (9), this polynomial can equivalently be expressed in terms of the function $H$ defined in (8) as

$${p}_{\sigma,\mu}(t) = \widetilde{H}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)$$

$$:= H\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)-\frac{\gamma^{\mathrm{mod}}\left(\frac{4\mu-1}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}.$$

Remark 2

We remark that, using (11), the product $p_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)$ appearing in the definition of $\mathcal{L}_{j}^{\mathrm{enr}}$ can be further written as

$$p_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)=\left(\widetilde{H}K_{\sigma,\mu}\right)\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right).$$

Moreover, we note that when $\mu=1$, equation (15) implies

$$\widetilde{k}_{\infty,1}(t)=\frac{1}{2},$$

and in this particular case the enriched functionals in (19) reduce to the classical unweighted histopolation functionals of (16).

We define the associated set of enriched linear functionals as

$$\Sigma^{\mathrm{enr}}_{1,\sigma,\mu}(T)=\left\{{\mathcal{I}}_{j}^{\mathrm{enr}},{\mathcal{L}}^{\mathrm{enr}}_{j}\,:\,j=1,2,3\right\}$$

and introduce the enriched triple

$$\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}=\left(T,\mathbb{P}_{2}(T),\Sigma^{\mathrm{enr}}_{1,\sigma,\mu}(T)\right),$$ (22)

where ${\mathbb{P}}_{2}(T)$ is the set of all quadratic polynomials defined on $T.$ This triple can be regarded as a generalization of the linear (unweighted) histopolation scheme. We employ the proposed enriched quadratic operators $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$ as building blocks for the function reconstruction method. Their quadratic nature allows them to represent curved features and ensures greater accuracy compared to the linear case.

The first step in our analysis is to show that, for any $\sigma>0$ and $\mu\geq 1$, the system $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$ defines a well-posed unisolvent triple for $\mathbb{P}_{2}(T)$. The second step is to determine the associated basis functions. These functions are then employed to define the first two-parameter family of local function reconstruction operators.

We now summarize our main findings together with some immediate consequences. Detailed proofs will be provided at the end of this subsection.

Theorem 2.1

The triple $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$ is unisolvent for any $\sigma>0$ and $\mu\geq 1$.

We next derive closed-form expressions for the basis functions associated with the enriched triple $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$. Specifically, we determine a basis

$${B}_{2,\sigma,\mu}=\left\{\varphi_{i},\psi_{i}\,:\,i=1,2,3\right\}$$

of the vector space $\mathbb{P}_{2}(T)$ satisfying the conditions

$${{\mathcal{I}}}_{j}^{\mathrm{enr}}\left(\varphi_{i}\right) =\delta_{ij},$$ (23)

$${{\mathcal{L}}}^{\mathrm{enr}}_{j}\left(\varphi_{i}\right) =0,$$ (24)

$${{\mathcal{I}}}_{j}^{\mathrm{enr}}\left(\psi_{i}\right) =0,$$ (25)

$${{\mathcal{L}}}^{\mathrm{enr}}_{j}\left(\psi_{i}\right) =\delta_{ij},$$ (26)

for any $i,j=1,2,3$.

Theorem 2.2

The basis functions associated with the enriched triple $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$ are given by

$$\varphi_{i} =1-2\lambda_{i},\quad i=1,2,3,$$ (27)

$$\psi_{i} =-A_{\sigma,\mu}\varphi_{i}+\frac{2}{\left\|{p}_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}\left(-\lambda_{i}^{2}+\lambda^{2}_{i+1}+\lambda^{2}_{i+2}\right),\quad i=1,2,3,$$ (28)

where

$$\left\|{p}_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}=\int_{-1}^{1}{p}^{2}_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)dt,$$

$$A_{\sigma,\mu}=\frac{1+t_{2,\sigma,\mu}}{\left\|{p}_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}},$$ (29)

and

$$t_{2,\sigma,\mu}=\int_{-1}^{1}t^{2}\widetilde{k}_{\sigma,\mu}(t)dt,$$

denotes the second moment of $\ \widetilde{k}_{\sigma,\mu}.$

We now have all the necessary ingredients to define the two-parameter local reconstruction operator associated with the pdfs $\widetilde{k}_{\sigma,\mu}$. This operator is constructed using the basis functions of Theorem 3.2 and is defined by

$$\begin{array}[]{rcl}{\pi}_{1,\sigma,\mu}^{{\mathrm{enr}}}:C(T)&\rightarrow&{\mathbb{P}}_{2}(T)\\ f&\mapsto&\displaystyle{\sum_{j=1}^{3}{\mathcal{I}}^{\mathrm{enr}}_{j}\left(f\right){\varphi}_{j}+\sum_{j=1}^{3}{\mathcal{L}}^{\mathrm{enr}}_{j}\left(f\right)}{\psi}_{j}.\end{array}$$ (30)

The proofs are based on certain properties of the pdfs $\widetilde{k}_{\sigma,\mu}$ and the polynomial $p_{\sigma,\mu}$. To establish Theorems 2.1 and 2.2, we first present several auxiliary lemmas and recall some concepts that will be useful in the sequel. In the discussion that follows, particular emphasis will be placed on the first three lemmas, which will be invoked repeatedly.

We begin with a lemma concerning the incomplete gamma function (3).

Lemma 1

For any $\rho>0$ and $s>-1$, the following identity holds

$$\int_{0}^{1}t^{s}e^{-(\rho t)^{2}}dt=\frac{1}{2}\gamma^{\mathrm{mod}}\left(\frac{s+1}{2},\rho^{2}\right),$$

where $\gamma^{\mathrm{mod}}$ denotes the modified incomplete gamma function defined in (5).

Proof

By using the change of variables $u=\rho^{2}t^{2}$, we obtain

$$\int_{0}^{1}t^{s}e^{-(\rho t)^{2}}dt=\frac{1}{2\rho^{s+1}}\int_{0}^{\rho^{2}}u^{\frac{s+1}{2}-1}e^{-u}du.$$

The result then follows directly from the definition (5) with parameters $\frac{s+1}{2}$ and $\rho^{2}$.

Exploiting the symmetry of the probability density $\widetilde{k}_{\sigma,\mu}$ defined in (13), we first state a number of its fundamental properties. The next lemma shows that $\widetilde{k}_{\sigma,\mu}$ is indeed a probability density function, that is, it has unit zeroth moment.

Lemma 2

For any $\mu\geq 1$ and $\sigma>0$, the function $\widetilde{k}_{\sigma,\mu}$ is a probability density function.

Proof

By symmetry and with the substitution $u=t^{\mu}$, we obtain

$$\int_{-1}^{1}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}}dt = 2\int_{0}^{1}t^{4\mu-4}e^{-\left(\frac{t^{\mu}}{\sqrt{2}\sigma^{\mu}}\right)^{2}}dt$$ (31)

$$= \frac{2}{\mu}\int_{0}^{1}u^{\frac{3\mu-3}{\mu}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{\mu}}\right)^{2}}du.$$

Applying Lemma 1 with

$$s=\frac{3\mu-3}{\mu},\quad\rho=\frac{1}{\sqrt{2}\sigma^{\mu}},$$

we can write

$$\frac{2}{\mu}\int_{0}^{1}u^{\frac{3\mu-3}{\mu}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{\mu}}\right)^{2}}du=\frac{1}{\mu}\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right).$$ (32)

Then, combining (31) and (32), it follows that

$$\int_{-1}^{1}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}}dt=\frac{1}{\mu}\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right).$$ (33)

Therefore, in order for $\widetilde{k}_{\sigma,\mu}$ to be a probability density function, the normalization constant $a_{\sigma,\mu}$ in (12) must be the reciprocal of the right-hand side of (33), so that

$$\int_{-1}^{1}\widetilde{k}_{\sigma,\mu}(t)dt=1.$$

The next lemma allows us to derive a compact formula for the moments of $\widetilde{k}_{\sigma,\mu}$ in terms of the modified incomplete gamma functions defined in (5).

Lemma 3

Let $\left\{t_{m,\sigma,\mu}\right\}_{m=0}^{\infty}$ be the sequence of moments of $\widetilde{k}_{\sigma,\mu}$. Then

$$t_{m,\sigma,\mu}=\int_{-1}^{1}t^{m}\widetilde{k}_{\sigma,\mu}(t)dt=\begin{cases}0,&\text{if $m=2k+1$ }\\ \dfrac{\gamma^{\mathrm{mod}}\left(\dfrac{2k+4\mu-3}{2\mu},\dfrac{1}{2\sigma^{2\mu}}\right)}{\gamma^{\mathrm{mod}}\left(\dfrac{4\mu-3}{2\mu},\dfrac{1}{2\sigma^{2\mu}}\right)},&\text{if $m=2k$}.\end{cases}$$ (34)

Proof

By symmetry, $\widetilde{k}_{\sigma,\mu}$ is even, hence for all $k\geq 0$, it follows that

$$t_{2k+1,\sigma,\mu}=0.$$

For the even moments, performing the change of variables $u=t^{\mu}$ gives

$$t_{2k,\sigma,\mu} = \int_{-1}^{1}t^{2k}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}}dt$$ (35)

$$= 2\int_{0}^{1}t^{2k+4\mu-4}e^{-\left(\frac{t^{\mu}}{\sqrt{2}\sigma^{\mu}}\right)^{2}}dt$$

$$= \frac{2}{\mu}\int_{0}^{1}u^{\frac{2k+3\mu-3}{\mu}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{\mu}}\right)^{2}}du.$$

Applying Lemma 1 with

$$s=\frac{2k+3\mu-3}{\mu},\quad\rho=\frac{1}{\sqrt{2}\sigma^{\mu}},$$

we can write

$$\frac{2}{\mu}\int_{0}^{1}u^{\frac{2k+3\mu-3}{\mu}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{\mu}}\right)^{2}}du=\frac{1}{\mu}\gamma^{\mathrm{mod}}\left(\frac{2k+4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right).$$ (36)

Then, combining (35) and (36), it follows that

$$\int_{-1}^{1}t^{2k}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}}dt=\frac{1}{\mu}\gamma^{\mathrm{mod}}\left(\frac{2k+4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right).$$ (37)

Then, by recalling the explicit expression of the pdfs $\widetilde{k}_{\sigma,\mu}$ in (13) and using (37), we obtain

$$t_{2k,\sigma,\mu} = \int_{-1}^{1}t^{2k}\widetilde{k}_{\sigma,\mu}(t)dt$$

$$= \frac{\mu}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}\int_{-1}^{1}t^{2k}\left(t^{2}\right)^{2\mu-2}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{\mu}}dt$$

$$= \frac{\gamma^{\mathrm{mod}}\left(\frac{2k+4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}.$$

Remark 3

By setting $k=0$ in (34), one immediately recovers the normalization condition, thus confirming that $\widetilde{k}_{\sigma,\mu}(t)$ is indeed a probability density function.

Remark 4

The polynomial ${p}_{\sigma,\mu}$ defined in (21) can be expressed in terms of the second moment as

$$p_{\sigma,\mu}(t)=t^{2}-t_{2,\sigma,\mu}.$$ (38)

Indeed, by taking $m=2$ in identity (34), one has

$$t_{2,\sigma,\mu}=\frac{\gamma^{\mathrm{mod}}\left(\frac{4\mu-1}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{4\mu-3}{2\mu},\frac{1}{2\sigma^{2\mu}}\right)},$$

which directly yields (38).

The limiting moments of $\widetilde{k}_{\sigma,\mu}(t)$ as $\sigma\to\infty$ are easily obtained from the following remark.

Remark 5

We observe that

$$\lim_{\sigma\to\infty}\int_{-1}^{1}t^{m}\widetilde{k}_{\sigma,\mu}(t)dt=\int_{-1}^{1}t^{m}\widetilde{k}_{\infty,\mu}(t)dt.$$

Indeed, from the definition (15) we have

$$\int_{-1}^{1}t^{m}\widetilde{k}_{\infty,\mu}(t)dt=\frac{4\mu-3}{2}\int_{-1}^{1}t^{m}\left(t^{2}\right)^{2\mu-2}dt.$$

If $m$ is odd, this integral vanishes by symmetry. If instead $m=2k$, a straightforward calculation yields

$$\int_{-1}^{1}t^{2k}\widetilde{k}_{\infty,\mu}(t)dt=\frac{4\mu-3}{2}\int_{-1}^{1}t^{2k}\left(t^{2}\right)^{2\mu-2}dt=\frac{4\mu-3}{2k+4\mu-3}.$$

On the other hand, using (34) together with (6), we find

$$t_{2k,\infty,\mu}:=\lim_{\sigma\to\infty}t_{2k,\sigma,\mu}=\frac{4\mu-3}{2k+4\mu-3},$$

which are exactly the moments of order $2k$ associated with the pdfs $\widetilde{k}_{\infty,\mu}(t)$.

The following result shows that $p_{\sigma,\mu}$, defined in (21), is indeed an orthogonal polynomial.

Lemma 4

The polynomial $p_{\sigma,\mu}(t)$ is orthogonal to all linear polynomials on $[-1,1]$ with respect to the pdfs $\widetilde{k}_{\sigma,\mu}$.

Proof

Since $p_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)$ is an even function, it suffices to show that

$$\int_{-1}^{1}p_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)dt=0.$$

Using the expression of $p_{\sigma,\mu}(t)$ from (38), we have

$$\int_{-1}^{1}p_{\sigma,\mu}(t)\widetilde{k}_{\sigma,\mu}(t)dt = \int_{-1}^{1}\left(t^{2}-t_{2,\sigma,\mu}\right)\widetilde{k}_{\sigma,\mu}(t)dt$$

$$= \int_{-1}^{1}t^{2}\widetilde{k}_{\sigma,\mu}(t)dt-t_{2,\sigma,\mu}\int_{-1}^{1}\widetilde{k}_{\sigma,\mu}(t)dt$$

$$= 0.$$

In the last equality we used the facts that

$$\int_{-1}^{1}\widetilde{k}_{\sigma,\mu}(t)dt=1,$$

since $\widetilde{k}_{\sigma,\mu}$ is a probability density function, and

$$\int_{-1}^{1}t^{2}\widetilde{k}_{\sigma,\mu}(t)dt=t_{2,\sigma,\mu},$$

which is precisely the second moment of $\widetilde{k}_{\sigma,\mu}$.

We will need the average values of the barycentric coordinates on the edges of $T$, weighted by the probability density $\widetilde{k}_{\sigma,\mu}(t)$.

Lemma 5

For any $i,j=1,2,3$, it holds

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}\right)=\frac{1}{2}\left(1-\delta_{ij}\right).$$ (39)

Proof

From (19), we have

$$\mathcal{I}_{j}^{\mathrm{enr}}\left(\lambda_{i}\right)=\int_{-1}^{1}\lambda_{i}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{k}_{\sigma,\mu}(t)dt.$$

If $i=j$, then by (2) it follows immediately that

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{j}\right)=0,\quad j=1,2,3.$$

If instead $i\neq j$, then one of the values $\lambda_{i}(\boldsymbol{v}_{j+1})$ or $\lambda_{i}(\boldsymbol{v}_{j+2})$ equals $1$, while the other equals $0$. Since $\widetilde{k}_{\sigma,\mu}$ is a probability density function and barycentric coordinates are affine functions, we get

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}\right) =\int_{-1}^{1}\left(\frac{1+t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)+\frac{1-t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)\right)\widetilde{k}_{\sigma,\mu}(t)dt$$

$$=\frac{1}{2}\int_{-1}^{1}\left((1+t)\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)+(1-t)\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)\right)\widetilde{k}_{\sigma,\mu}(t)dt$$

$$=\frac{1}{2}\int_{-1}^{1}\widetilde{k}_{\sigma,\mu}(t)dt=\frac{1}{2},$$

where in the last equality we have used the fact that $\widetilde{k}_{\sigma,\mu}$ is a probability density. Then the lemma is proved.

We also need explicit expressions for the average values of the square of the barycentric coordinates on each edge of $T$, weighted by $\widetilde{k}_{\sigma,\mu}(t)$.

Lemma 6

For any $i,j=1,2,3$, it holds

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\frac{1+t_{2,\sigma,\mu}}{4}\left(1-\delta_{ij}\right),$$ (40)

where $t_{2,\sigma,\mu}$ is the second moment of $\widetilde{k}_{\sigma,\mu}$.

Proof

From (19), we have

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\int_{-1}^{1}\lambda_{i}^{2}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{k}_{\sigma,\mu}(t)dt.$$

If $i=j$, then by (2), it follows immediately that

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{j}^{2}\right)=0,\quad j=1,2,3.$$

If instead $i\neq j$, we obtain

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right) =\int_{-1}^{1}\left(\frac{1+t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)+\frac{1-t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)\right)^{2}\widetilde{k}_{\sigma,\mu}(t)dt$$

$$=\frac{1}{4}\int_{-1}^{1}\left((1+t)\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)+(1-t)\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)\right)^{2}\widetilde{k}_{\sigma,\mu}(t)dt.$$ (41)

Proceeding as in the proof of Lemma 5, the integral in (41) decomposes as

$${\mathcal{I}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\frac{1}{4}\int_{-1}^{1}\widetilde{k}_{\sigma,\mu}(t)dt+\frac{1}{4}\int_{-1}^{1}t^{2}\widetilde{k}_{\sigma,\mu}(t)dt.$$

Since $\widetilde{k}_{\sigma,\mu}$ is a probability density, the first integral equals $1$, while the second equals $t_{2,\sigma,\mu}$. This proves (40).

The following result is an immediate consequence of Lemma 4.

Lemma 7

For any $i,j=1,2,3$, the following identity holds

$${\mathcal{L}}_{j}^{\mathrm{enr}}\left(\lambda_{i}\right)=0.$$ (42)

Proof

The result follows directly from the definition of $\mathcal{L}_{j}^{\mathrm{enr}}$ given in (20), together with the orthogonality property of the polynomial $p_{\sigma,\mu}(t)$ established in Lemma 4.

Remark 6

From the definition of $\mathcal{L}_{j}^{\mathrm{enr}}$ and the orthogonality of the polynomial $p_{\sigma,\mu}$ established in Lemma 4, it follows that

$$\mathcal{L}_{j}^{\mathrm{enr}}\left(p_{1}\right)=0,\quad j=1,2,3,$$

for every linear polynomial $p_{1}\in\mathbb{P}_{1}(T)$.

Lemma 8

For any $i,j=1,2,3$, the following identity holds

$${\mathcal{L}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\frac{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}{4}\left(1-\delta_{ij}\right).$$ (43)

Proof

From equation (20) and the fact that barycentric coordinates are affine functions, we obtain

$${\mathcal{L}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right) =\int_{-1}^{1}p_{\sigma,\mu}(t)\lambda_{i}^{2}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{k}_{\sigma,\mu}(t)dt$$

$$=\int_{-1}^{1}p_{\sigma,\mu}(t)\left(\frac{1+t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)+\frac{1-t}{2}\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)\right)^{2}\widetilde{k}_{\sigma,\mu}(t)dt.$$

If $i=j$, then by (2) it follows that

$${\mathcal{L}}_{j}^{\mathrm{enr}}\left(\lambda_{j}^{2}\right)=0,\quad j=1,2,3.$$

If instead $i\neq j$, then $\lambda_{i}\left(\boldsymbol{v}_{j+1}\right)=1$ or $\lambda_{i}\left(\boldsymbol{v}_{j+2}\right)=1$, and one of them must be zero. Using this fact, and exploiting the orthogonality property of $p_{\sigma,\mu}$ established in Lemma 4, together with the identity (38) expressing $t^{2}$ in terms of $p_{\sigma,\mu}$, we obtain

$${\mathcal{L}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right) =\frac{1}{4}\int_{-1}^{1}p_{\sigma,\mu}(t)t^{2}\widetilde{k}_{\sigma,\mu}(t)dt$$

$$=\frac{1}{4}\int_{-1}^{1}p_{\sigma,\mu}(t)\left(p_{\sigma,\mu}(t)+t_{2,\sigma,\mu}\right)\widetilde{k}_{\sigma,\mu}(t)dt.$$

Using again the orthogonality property of $p_{\sigma,\mu}$ the desired result follows.

We are now ready for the proof that the triple $\mathcal{H}^{\mathrm{enr}}_{1,\sigma,\mu}$ is unisolvent. To this aim, we first introduce some notations that will be used throughout the remainder of the paper. In particular, we consider the matrix

$$A:=\begin{bmatrix}0&1&1\\ 1&0&1\\ 1&1&0\end{bmatrix}.$$ (44)

As $\det(A)=2\neq 0$, the matrix $A$ is invertible, with inverse given by

$$A^{-1}=\frac{1}{2}\begin{bmatrix}-1&1&1\\ 1&-1&1\\ 1&1&-1\end{bmatrix}.$$ (45)

Proof of Theorem 2.1. It is sufficient to prove that the only polynomial $p_{2}\in\mathbb{P}_{2}(T)$ satisfying

$$\mathcal{I}_{j}^{\mathrm{enr}}\!\left(p_{2}\right) =0,\quad j=1,2,3,$$ (46)

$$\mathcal{L}_{j}^{\mathrm{enr}}\!\left(p_{2}\right) =0,\quad j=1,2,3,$$ (47)

is the zero polynomial. To this aim, let $p_{2}\in\mathbb{P}_{2}(T)$, which can be expressed in barycentric form as

$$p_{2}=\sum_{i=1}^{3}a_{i}\lambda_{i}+\sum_{i=1}^{3}b_{i}\lambda_{i}^{2},$$ (48)

with coefficients $a_{i},b_{i}\in\mathbb{R}$, $i=1,2,3$. By Lemmas 7 and 8, condition (47) reduces to

$$0={\mathcal{L}}^{\mathrm{enr}}_{j}\left(p_{2}\right) =\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}b_{i}\mathcal{L}^{\mathrm{enr}}_{j}\left(\lambda^{2}_{i}\right)=\frac{\left\|p_{\sigma,\mu}\right\|^{2}_{{{2,\sigma,\mu}}}}{4}\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}b_{i}\left(1-\delta_{ij}\right).$$

This can equivalently be written in matrix form as

$$\dfrac{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}{4}A{\boldsymbol{b}}={\boldsymbol{0}},$$

where $A$ is the matrix defined in (44), and

$${\boldsymbol{b}}^{T}=\left(b_{1},b_{2},b_{3}\right),\quad{\boldsymbol{0}}^{T}=(0,0,0).$$

Since $A$ is non-singular, we conclude that

$$b_{1}=b_{2}=b_{3}=0.$$

Consequently, by (48), the polynomial $p_{2}$ reduces to

$$p_{2}=p_{1}=a_{1}\lambda_{1}+a_{2}\lambda_{2}+a_{3}\lambda_{3},\quad a_{1},a_{2},a_{3}\in\mathbb{R}.$$

Applying Lemma 5, condition (46) becomes

$$0 ={\mathcal{I}}_{j}^{\mathrm{enr}}\left(p_{2}\right)=\frac{1}{2}\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}a_{i}\left(1-\delta_{ij}\right).$$

This can equivalently be written in matrix form as

$$\frac{1}{2}A{\boldsymbol{a}}={\boldsymbol{0}},$$ (49)

where

$${\boldsymbol{a}}^{T}=\left(a_{1},a_{2},a_{3}\right).$$

Since $A$ is non-singular, system (49) admits only the trivial solution

$$a_{1}=a_{2}=a_{3}=0,$$

which implies $p_{2}=0$ and thus completes the proof. ∎

The proof of Theorem 3.2 relies on the derivation of the explicit analytic expressions of the basis functions $\varphi_{i}$ and $\psi_{i}$. We illustrate the procedure by computing $\varphi_{1}$ and $\psi_{1}$, as the remaining cases follow in the same way.

Proof of Theorem 2.2.
We begin by deriving the analytic expression of the basis function $\varphi_{1}$. According to its definition, $\varphi_{1}$ can be expressed as

$$\varphi_{1}=\sum_{i=1}^{3}a_{i}\lambda_{i}+\sum_{i=1}^{3}b_{i}\lambda_{i}^{2},$$ (50)

for some coefficients $a_{i},b_{i}\in\mathbb{R}$, $i=1,2,3$. To determine these coefficients, we use condition (24) together with Lemmas 7 and 8, which yield

$$0=\mathcal{L}^{\mathrm{enr}}_{j}\left(\varphi_{1}\right)=\sum_{i=1}^{3}b_{i}\mathcal{L}^{\mathrm{enr}}_{j}\left(\lambda^{2}_{i}\right)=\frac{\|p_{\sigma,\mu}\|^{2}_{2,{\sigma,\mu}}}{4}\sum_{i=1}^{3}b_{i}\left(1-\delta_{ij}\right).$$

As in the previous proof, this system admits only the trivial solution

$$b_{1}=b_{2}=b_{3}=0.$$

Substituting these values into (50), $\varphi_{1}$ reduces to

$$\varphi_{1}=\sum_{i=1}^{3}a_{i}\lambda_{i}.$$

Applying Lemma 5, condition (23) becomes

$$\delta_{1j} ={\mathcal{I}}^{\mathrm{enr}}_{j}\left(\varphi_{1}\right)=\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}a_{i}{\mathcal{I}}^{\mathrm{enr}}_{j}\left(\lambda_{i}\right)=\frac{1}{2}\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}a_{i}\left(1-\delta_{ij}\right).$$

This is equivalently written in matrix form as

$$\left(1,0,0\right)^{T}=\frac{1}{2}A{\boldsymbol{a}}$$

where

$${\boldsymbol{a}}^{T}=\left(a_{1},a_{2},a_{3}\right).$$

A straightforward computation yields

$$a_{1}=-1,\quad a_{2}=a_{3}=1,$$

and hence

$$\varphi_{1}=1-2\lambda_{1}.$$

It remains to compute the expression of the basis function $\psi_{1}$. By definition, it can be expressed as

$$\psi_{1}=\sum_{i=1}^{3}c_{i}\lambda_{i}+\sum_{i=1}^{3}d_{i}\lambda_{i}^{2},$$ (51)

for some coefficients $c_{i},d_{i}\in\mathbb{R}$, $i=1,2,3$. To determine these coefficients, we use condition (26) together with Lemmas 7 and 8, which yield

$$\delta_{1j} ={\mathcal{L}}^{\mathrm{enr}}_{j}\left(\psi_{1}\right)=\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}d_{i}{\mathcal{L}}^{\mathrm{enr}}_{j}\left(\lambda_{i}^{2}\right)=\frac{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}{4}\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}d_{i}\left(1-\delta_{ij}\right).$$

This is equivalently written in matrix form as

$$\left(1,0,0\right)^{T}=\frac{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}{4}A{\boldsymbol{d}}$$

where

$${\boldsymbol{d}}^{T}=\left(d_{1},d_{2},d_{3}\right).$$

Solving this system, we have

$$d_{1}=-\frac{2}{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}},\quad d_{2}=d_{3}=\frac{2}{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}.$$

Substituting these values into (51), $\psi_{1}$ reduces to

$$\psi_{1}=\sum_{i=1}^{3}c_{i}\lambda_{i}+\frac{2}{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}\left(-\lambda_{1}^{2}+\lambda^{2}_{2}+\lambda^{2}_{3}\right).$$

Applying Lemmas 5 and 6, condition (25) yields

$$0 ={\mathcal{I}}^{\mathrm{enr}}_{j}\left(\psi_{1}\right)=\frac{1}{2}\sum_{\begin{subarray}{c}i=1\end{subarray}}^{3}c_{i}\left(1-\delta_{ij}\right)+A_{\sigma,\mu}\left(1-\delta_{ij}\right),$$

where $A_{\sigma,\mu}$ is defined in (29). This system admits the unique solution

$$c_{1}=A_{\sigma,\mu},\quad c_{2}=c_{3}=-A_{\sigma,\mu}.$$

Hence $\psi_{1}$ takes the form

$$\psi_{1}=-A_{\sigma,\mu}\varphi_{1}+\frac{2}{\left\|p_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}\left(-\lambda_{1}^{2}+\lambda^{2}_{2}+\lambda^{2}_{3}\right),$$

which concludes the proof. ∎

2.4 Second family of distributions

We now present a second two-parameter family of exponential-type distributions, which generalizes the truncated normal distribution, and is given by the following bivariate function with parameters $\mu\geq 1$ and $\sigma>0$

$$G_{\sigma,\mu}(\boldsymbol{x})=\frac{(2\mu-1)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}H^{\mu-1}(\boldsymbol{x})e^{-\frac{1}{2}\left(\frac{H(\boldsymbol{x})}{\sigma^{2}}\right)^{2\mu-1}},$$ (52)

where $H$ is given by (8). Unlike the function of the first family, $K_{\sigma,\mu}$ in (7), the function $G_{\sigma,\mu}$ does not, in general, define a weight function for all values of $\mu$, because $H$ changes sign on $T$. In fact, at the barycenter

$$\bar{\boldsymbol{v}}:=\frac{1}{3}\left(\boldsymbol{v}_{1}+\boldsymbol{v}_{2}+\boldsymbol{v}_{3}\right)$$

one has $H(\bar{\boldsymbol{v}})=-\frac{1}{3}$, hence

$$H^{\mu-1}(\bar{\boldsymbol{v}})=\left(-\frac{1}{3}\right)^{\mu-1},$$

which is negative whenever $\mu$ is an integer and $\mu-1$ is odd.

With the normalization in (52), the restriction of $G_{\sigma,\mu}$ to any edge of $T$ yields a probability density function on $[-1,1]$ (see Lemma 9). Indeed, using (1) and (2) we obtain

$$\widetilde{g}_{\sigma,\mu}(t) =G_{\sigma,\mu}\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)$$

$$=\frac{2\mu-1}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}\left(t^{2}\right)^{\mu-1}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{2\mu-1}},\quad t\in[-1,1].$$ (53)

Remark 7

When $\mu=1$, the edge-restricted densities $\widetilde{g}_{\sigma,1}(t)$ and $\widetilde{k}_{\sigma,1}(t)$ coincide for $t\in[-1,1]$. In particular, both reduce to the standard truncated normal distribution on $[-1,1]$.

Figure 2: Normalized probability density functions $\widetilde{g}_{\sigma,\mu}$ on $[-1,1]$. Left: effect of the scale parameter, with fixed $\mu=2$ and different values of $\sigma$. Right: effect of the shape parameter, with fixed $\sigma=0.50$ and varying values of $\mu$, illustrating the increased sharpness and concentration around $t=0$. These plots emphasize the flexibility provided by the two parameters.

Figure 2 shows that $\widetilde{g}_{\sigma,\mu}$ depends strongly on both the shape parameter $\mu$ and the scale parameter $\sigma$, thereby motivating the two-parameter generalization.

To define our second two-parameter family of nonconforming histopolation operators, we introduce the following enriched linear functionals

$${\mathcal{D}}_{j}^{\mathrm{enr}}(f) =\int_{-1}^{1}f\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{g}_{\sigma,\mu}(t)dt,\quad j=1,2,3,$$ (54)

$${\mathcal{S}}_{j}^{\mathrm{enr}}\left(f\right) =\int_{-1}^{1}q_{\sigma,\mu}(t)f\left(\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2}\right)\widetilde{g}_{\sigma,\mu}(t)dt,\quad j=1,2,3.$$ (55)

where

$${q}_{\sigma,\mu}(t)=t^{2}-\frac{\gamma^{\mathrm{mod}}\left(\frac{2\mu+1}{2(2\mu-1)},\frac{1}{2\sigma^{4\mu-2}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}.$$ (56)

We then define the following triple

$$\mathcal{H}^{\mathrm{enr}}_{2,\sigma,\mu}=\left(T,\mathbb{P}_{2}(T),\Sigma_{2,\sigma,\mu}^{\mathrm{enr}}(T)\right),$$

where

$$\Sigma^{\mathrm{enr}}_{2,\sigma,\mu}(T)=\left\{{\mathcal{D}}_{j}^{\mathrm{enr}},{{\mathcal{S}}}^{\mathrm{enr}}_{j}\,:\,j=1,2,3\right\}.$$

Remark 8

The distribution $\widetilde{g}_{\sigma,\mu}$ admits a beta-type (power) density as a limiting case. Indeed, using (4) with

$$s=\frac{1}{2},\quad z=\frac{1}{2\sigma^{4\mu-2}}$$

it follows that

$$\lim_{z\rightarrow 0^{+}}\gamma^{\mathrm{mod}}\left(\frac{1}{2},z\right)=\lim_{z\rightarrow 0^{+}}\frac{\gamma\left(\frac{1}{2},z\right)}{\sqrt{z}}=2.$$

Therefore, by (53), for all $t\in[-1,1]$ we obtain the pointwise limit

$$\widetilde{g}_{\infty,\mu}(t):=\lim_{\sigma\rightarrow\infty}\widetilde{g}_{\sigma,\mu}(t)=\frac{2\mu-1}{2}\left(t^{2}\right)^{\mu-1}.$$ (57)

Hence, $\widetilde{g}_{\infty,\mu}$ is the symmetric beta-type (power) density on $[-1,1]$, highlighting the limiting behaviour of the family as $\sigma\to\infty$.

Our first goal is to prove that $\mathcal{H}_{2,\sigma,\mu}^{\mathrm{enr}}$ is unisolvent. In a second step, we characterize the basis functions associated with this triple. Proofs are deferred to later subsections. The main results and some immediate consequences are summarized below.

Theorem 2.3

The triple $\mathcal{H}^{\mathrm{enr}}_{2,\sigma,\mu}$ is unisolvent for any $\sigma>0$ and $\mu\geq 1$.

We next derive closed-form expressions for the basis functions associated with the triple $\mathcal{H}^{\mathrm{enr}}_{2,\sigma,\mu}$. Specifically, we determine a basis

$${B}^{\prime}_{2,\sigma,\mu}=\left\{\eta_{i},\rho_{i}\,:\,i=1,2,3\right\}$$

of the vector space $\mathbb{P}_{2}(T)$ satisfying the conditions

$${{\mathcal{D}}}_{j}^{\mathrm{enr}}\left(\eta_{i}\right) =\delta_{ij},$$

$${{\mathcal{S}}}^{\mathrm{enr}}_{j}\left(\eta_{i}\right) =0,$$

$${\mathcal{D}}_{j}^{\mathrm{enr}}\left(\rho_{i}\right) =0,$$

$${{\mathcal{S}}}^{\mathrm{enr}}_{j}\left(\rho_{i}\right) =\delta_{ij},$$

for any $i,j=1,2,3$.

Theorem 2.4

The basis functions associated with the enriched triple $\mathcal{H}^{\mathrm{enr}}_{2,\sigma,\mu}$ are given by

$$\eta_{i} =1-2\lambda_{i},\quad i=1,2,3,$$

$$\rho_{i} =-C_{\sigma,\mu}\eta_{i}+\frac{2}{\left\|{q}_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}\left(-\lambda_{i}^{2}+\lambda^{2}_{i+1}+\lambda^{2}_{i+2}\right),\quad i=1,2,3,$$

where

$$\left\|{q}_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}=\int_{-1}^{1}{q}^{2}_{\sigma,\mu}(t)\widetilde{g}_{\sigma,\mu}(t)dt,$$

$$C_{\sigma,\mu}=\frac{1+s_{2,\sigma,\mu}}{\left\|q_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}},$$

and

$$s_{2,\sigma,\mu}=\int_{-1}^{1}t^{2}\widetilde{g}_{\sigma,\mu}(t)dt,$$

denotes the second moment of $\ \widetilde{g}_{\sigma,\mu}.$

We now have all the necessary ingredients to define the two-parameter local reconstruction operator associated with the pdfs $\widetilde{g}_{\sigma,\mu}$. This operator is constructed using the basis functions of Theorem 2.4 and is defined by

$$\begin{array}[]{rcl}{\pi}_{2,\sigma,\mu}^{{\mathrm{enr}}}:C(T)&\rightarrow&{\mathbb{P}}_{2}(T)\\ f&\mapsto&\displaystyle{\sum_{j=1}^{3}{\mathcal{D}}^{\mathrm{enr}}_{j}\left(f\right){\eta}_{j}+\sum_{j=1}^{3}{\mathcal{S}}^{\mathrm{enr}}_{j}\left(f\right)}{\rho}_{j}.\end{array}$$

As in the previous case, the proofs rely on specific properties of the density $\widetilde{g}_{\sigma,\mu}$ and of the polynomial $q_{\sigma,\mu}$. To prove Theorems 2.3 and 2.4, we first collect a number of auxiliary lemmas and recall some standard notions needed in the sequel. Exploiting the symmetry of the edge density $\widetilde{g}_{\sigma,\mu}$ defined in (53), we begin by recording some of its basic properties. The next lemma shows that $\widetilde{g}_{\sigma,\mu}$ is a probability density on $[-1,1]$.

Lemma 9

For any $\mu\geq 1$ and $\sigma>0$, the function $\widetilde{g}_{\sigma,\mu}$ is a probability density function.

Proof

By symmetry and with the substitution $u=t^{2\mu-1}$, we obtain

$$\int_{-1}^{1}\left(t^{2}\right)^{\mu-1}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{2\mu-1}}dt = 2\int_{0}^{1}t^{2\mu-2}e^{-\left(\frac{t^{2\mu-1}}{\sqrt{2}\sigma^{2\mu-1}}\right)^{2}}dt$$ (58)

$$= \frac{2}{2\mu-1}\int_{0}^{1}e^{-\left(\frac{u}{\sqrt{2}\sigma^{2\mu-1}}\right)^{2}}du.$$

Applying Lemma 1 with

$$s=0,\quad\rho=\frac{1}{\sqrt{2}\sigma^{2\mu-1}},$$

we can write

$$\frac{2}{2\mu-1}\int_{0}^{1}e^{-\left(\frac{u}{\rho}\right)^{2}}du=\frac{1}{2\mu-1}\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right).$$ (59)

Then, combining (58) and (59), it follows that

$$\int_{-1}^{1}\left(t^{2}\right)^{\mu-1}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{2\mu-1}}dt=\frac{1}{2\mu-1}\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right).$$ (60)

Therefore, by the definition of $\widetilde{g}_{\sigma,\mu}$ in (53), whose normalization constant equals

$$\frac{2\mu-1}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)},$$

we get

$$\int_{-1}^{1}\widetilde{g}_{\sigma,\mu}(t)dt=1.$$

This proves the claim.

The next lemma allows us to derive a compact formula for the moments of $\widetilde{g}_{\sigma,\mu}$ in terms of the modified incomplete gamma functions defined in (5).

Lemma 10

Let $\{s_{m,\sigma,\mu}\}_{m=0}^{\infty}$ be the sequence of normalized moments of $\widetilde{g}_{\sigma,\mu}$. Then

$$s_{m,\sigma,\mu}=\int_{-1}^{1}t^{m}\widetilde{g}_{\sigma,\mu}(t)dt=\begin{cases}0,&\text{if }m=2k+1,\\[2.84526pt] \dfrac{\gamma^{\mathrm{mod}}\left(\dfrac{2k+2\mu-1}{2(2\mu-1)},\dfrac{1}{2\sigma^{4\mu-2}}\right)}{\gamma^{\mathrm{mod}}\left(\dfrac{1}{2},\dfrac{1}{2\sigma^{4\mu-2}}\right)},&\text{if }m=2k.\end{cases}$$ (61)

Proof

By symmetry, $\widetilde{g}_{\sigma,\mu}$ is even, hence for all $k\geq 0$, it follows that

$$s_{2k+1,\sigma,\mu}=0.$$

For the even moments, performing the change of variables $u=t^{2\mu-1}$ gives

$$s_{2k,\sigma,\mu} = \frac{2\mu-1}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}\int_{-1}^{1}t^{2k}\left(t^{2}\right)^{\mu-1}e^{-\frac{1}{2}\left(\frac{t^{2}}{\sigma^{2}}\right)^{2\mu-1}}dt$$

$$= \frac{2(2\mu-1)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}\int_{0}^{1}t^{2k+2\mu-2}e^{-\left(\frac{t^{2\mu-1}}{\sqrt{2}\sigma^{2\mu-1}}\right)^{2}}dt$$

$$= \frac{2}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}\int_{0}^{1}u^{\frac{2k}{2\mu-1}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{2\mu-1}}\right)^{2}}du.$$

Applying Lemma 1 with

$$s=\frac{2k}{2\mu-1},\quad\rho=\frac{1}{\sqrt{2}\sigma^{2\mu-1}}$$

we can write

$$\frac{2}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}\int_{0}^{1}u^{\frac{2k}{2\mu-1}}e^{-\left(\frac{u}{\sqrt{2}\sigma^{2\mu-1}}\right)^{2}}du=\frac{\gamma^{\mathrm{mod}}\left(\frac{2k+2\mu-1}{2(2\mu-1)},\frac{1}{2\sigma^{4\mu-2}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}.$$

Therefore

$$s_{2k,\sigma,\mu}=\frac{\gamma^{\mathrm{mod}}\left(\frac{2k+2\mu-1}{2(2\mu-1)},\frac{1}{2\sigma^{4\mu-2}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)}.$$

As an immediate application of the integral representation of the moments in Lemma 10, we obtain an alternative expression for the polynomial in (56).

Remark 9

The polynomial ${q}_{\sigma,\mu}$ defined in (56) can be expressed in terms of the second moment as

$$q_{\sigma,\mu}(t)=t^{2}-s_{2,\sigma,\mu}.$$ (62)

Indeed, by taking $m=2$ in identity (61), one has

$$s_{2,\sigma,\mu}=\frac{\gamma^{\mathrm{mod}}\left(\frac{2\mu+1}{2(2\mu-1)},\frac{1}{2\sigma^{4\mu-2}}\right)}{\gamma^{\mathrm{mod}}\left(\frac{1}{2},\frac{1}{2\sigma^{4\mu-2}}\right)},$$

which directly yields (62).

Remark 10

Setting $m=0$ in (61) gives $s_{0,\sigma,\mu}=1$, confirming that $\widetilde{g}_{\sigma,\mu}$ is a probability density.

Remark 11

We observe that

$$\lim_{\sigma\to\infty}\int_{-1}^{1}t^{m}\widetilde{g}_{\sigma,\mu}(t)dt=\int_{-1}^{1}t^{m}\widetilde{g}_{\infty,\mu}(t)dt.$$

Indeed, from the definition (57) we have

$$\int_{-1}^{1}t^{m}\widetilde{g}_{\infty,\mu}(t)dt=\frac{2\mu-1}{2}\int_{-1}^{1}t^{m}\left(t^{2}\right)^{\mu-1}dt.$$

If $m$ is odd, this integral vanishes by symmetry. If instead $m=2k$, a straightforward calculation yields

$$\int_{-1}^{1}t^{2k}\widetilde{g}_{\infty,\mu}(t)dt=\frac{2\mu-1}{2}\int_{-1}^{1}t^{2k}\left(t^{2}\right)^{\mu-1}dt=\frac{2\mu-1}{2k+2\mu-1}.$$

On the other hand, using (61) together with (6), we find

$$s_{2k,\infty,\mu}:=\lim_{\sigma\to\infty}s_{2k,\sigma,\mu}=\frac{2\mu-1}{2k+2\mu-1},$$

which are exactly the moments of order $2k$ associated with the pdfs $\widetilde{g}_{\infty,\mu}(t)$.

The next result establishes that $q_{\sigma,\mu}$, defined in (56), is the degree-two orthogonal polynomial associated with the weight $\widetilde{g}_{\sigma,\mu}$ on $[-1,1]$.

Lemma 11

The polynomial $q_{\sigma,\mu}(t)$ is orthogonal to all linear polynomials on $[-1,1]$ w.r.t. the weight function $\widetilde{g}_{\sigma,\mu}$.

Proof

Since $q_{\sigma,\mu}(t)\widetilde{g}_{\sigma,\mu}(t)$ is an even function, it suffices to show that

$$\int_{-1}^{1}q_{\sigma,\mu}(t)\widetilde{g}_{\sigma,\mu}(t)dt=0.$$ (63)

Using the expression of $q_{\sigma,\mu}$ from (62), we have

$$\int_{-1}^{1}q_{\sigma,\mu}(t)\widetilde{g}_{\sigma,\mu}(t)dt = \int_{-1}^{1}\left(t^{2}-s_{2,\sigma,\mu}\right)\widetilde{g}_{\sigma,\mu}(t)dt$$ (64)

$$= \int_{-1}^{1}t^{2}\widetilde{g}_{\sigma,\mu}(t)dt-s_{2,\sigma,\mu}\int_{-1}^{1}\widetilde{g}_{\sigma,\mu}(t)dt$$ (65)

$$= 0.$$ (66)

In the last equality we used the facts that

$$\int_{-1}^{1}\widetilde{g}_{\sigma,\mu}(t)dt=1$$

since $\widetilde{g}_{\sigma,\mu}(t)$ is a probability density function, and

$$\int_{-1}^{1}t^{2}\widetilde{g}_{\sigma,\mu}(t)dt=s_{2,\sigma,\mu}$$

which is precisely the second moment of $\widetilde{g}_{\sigma,\mu}(t)$.

In what follows, we prove several crucial lemmas needed in the sequel.

Lemma 12

For any $i,j=1,2,3$, it holds

$${\mathcal{D}}_{j}^{\mathrm{enr}}\left(\lambda_{i}\right)=\frac{1}{2}\left(1-\delta_{ij}\right).$$ (67)

Proof

The proof is similar to that of Lemma 5 and is therefore omitted.

Lemma 13

For any $i,j=1,2,3$, it holds

$${\mathcal{D}}_{j}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\frac{1+s_{2,\sigma,\mu}}{4}\left(1-\delta_{ij}\right),$$ (68)

where $s_{2,\sigma,\mu}$ is the second moment of $\widetilde{g}_{\sigma,\mu}$.

Proof

The proof is similar to that of Lemma 6 and is therefore omitted.

Lemma 14

For any $i,j=1,2,3$, it holds

$${\mathcal{S}}_{{j}}^{\mathrm{enr}}\left(\lambda_{i}\right)=0.$$ (69)

Proof

The result follows directly from (55) and Lemma 11.

Remark 12

From the definition of $\mathcal{S}_{j}^{\mathrm{enr}}$ and the orthogonality of the polynomial $q_{\sigma,\mu}$ established in Lemma 11, it follows that

$$\mathcal{S}_{j}^{\mathrm{enr}}\left(p_{1}\right)=0,\quad j=1,2,3,$$ (70)

for every linear polynomial $p_{1}\in\mathbb{P}_{1}(T)$.

Lemma 15

For any $i,j=1,2,3$, it holds

$${\mathcal{S}}_{{j}}^{\mathrm{enr}}\left(\lambda_{i}^{2}\right)=\frac{\left\|q_{\sigma,\mu}\right\|^{2}_{2,{{\sigma,\mu}}}}{4}\left(1-\delta_{ij}\right).$$ (71)

Proof

The proof is similar to that of Lemma 8 and is therefore omitted.

Now, we are ready to prove that the triple $\mathcal{H}_{2,\sigma,\mu}^{\mathrm{enr}}$ is unisolvent.

Proof of Theorem 2.3. Using Lemmas 12, 13, 14, and 15, the proof follows along the same lines as that of Theorem 2.1 and is therefore omitted. ∎

Proof of Theorem 2.4
Using Lemmas 12, 13, 14, and 15, the proof follows along the same lines as that of Theorem 2.2 and is therefore omitted.

∎

2.5 Parameter selection by global mesh-level tuning

We discuss a complementary strategy to select the parameters of the edge densities (e.g., $\mu,\sigma$ for the families introduced in Section 2). Given a validation set of test functions

$$\left\{f^{(r)}\right\}_{r=1}^{R}$$

and a sequence of meshes

$$\left\{\mathcal{T}_{n}\right\}_{n=0}^{N},$$

we determine the optimal parameters by minimizing the global $L^{1}$ error over a candidate grid $\mathcal{P}$ of parameter pairs:

$$\left(\mu^{\star},\sigma^{\star}\right)\in\arg\min_{(\mu,\sigma)\in\mathcal{P}}\sum_{r}\sum_{n}\left\|f^{(r)}-{\pi}_{1,\sigma,\mu}^{{\mathrm{enr}}}\left[f^{(r)}\right]\right\|_{L^{1}(\Omega;\mathcal{T}_{n})},$$

where ${\pi}_{1,\sigma,\mu}^{{\mathrm{enr}}}$ is defined in (30). This procedure provides a robust, once-for-all choice of the parameters.

Algorithm 1 Global parameter tuning by grid search

1:Validation functions $\left\{f^{(r)}\right\}_{r=1}^{R}$; meshes $\{\mathcal{T}_{n}\}_{n=0}^{N}$; candidate grid $\mathcal{P}$ of parameter pairs $(\mu,\sigma)$; reconstruction operator $\pi^{\mathrm{enr}}_{1,\mu,\sigma}$

2:Optimal parameters $\left(\mu^{\star},\sigma^{\star}\right)$

3:$E_{\min}\leftarrow+\infty$;  $\left(\mu^{\star},\sigma^{\star}\right)\leftarrow\text{undefined}$

4:for all $(\mu,\sigma)\in\mathcal{P}$ do

5:  $E(\mu,\sigma)\leftarrow 0$

6:  for $r=1$ to $R$ do

7:   for $n=0$ to $N$ do

8:     $u\leftarrow\pi^{\mathrm{enr}}_{1,\mu,\sigma}[f^{(r)};\mathcal{T}_{n}]$

9:     $E(\mu,\sigma)\leftarrow E(\mu,\sigma)+\|f^{(r)}-u\|_{L^{1}(\Omega;\mathcal{T}_{n})}$

10:   end for

11:  end for

12:  if $E(\mu,\sigma)<E_{\min}$ then

13:   $E_{\min}\leftarrow E(\mu,\sigma)$;  $\left(\mu^{\star},\sigma^{\star}\right)\leftarrow(\mu,\sigma)$

14:  end if

15:end for

16:return $\left(\mu^{\star},\sigma^{\star}\right)$

This grid–search strategy is simple to implement and guarantees the existence of a well–defined pair of optimal parameters $\left(\mu^{\star},\sigma^{\star}\right)$. Once selected, these values can be applied uniformly across all mesh elements, thus ensuring that the reconstruction procedure is completely specified and reproducible. Although the search involves multiple reconstructions on the validation set, it is performed only once as an offline preprocessing step. The modest computational cost is largely compensated by a substantial gain in robustness and accuracy in the proposed weighted histopolation framework. In practice, the method behaves like a standard hyperparameter tuning stage in machine learning, providing an automatic and principled way to calibrate the edge densities before applying the reconstruction algorithm to new data.

3 Generalization to Arbitrary Edge Probability Densities

In the previous sections, we focused on two specific two-parameter families of generalized truncated normal distributions, each giving rise to valid probability densities on the edges of a triangular element. A key observation is that, regardless of the chosen family, the resulting edge-restricted densities share the same structural properties: they are normalized probability densities, symmetric with respect to the edge parameter, and admit finite moments of all orders. These properties alone are sufficient to establish unisolvence of the enriched triple and to construct explicit quadratic reconstruction operators.

This motivates us to consider a more general setting. Let $T\subset\mathbb{R}^{2}$ be a nondegenerate triangle with edges $s_{j}=\left[\boldsymbol{v}_{j+1},\boldsymbol{v}_{j+2}\right]$, and denote by

$$\gamma_{j}(t)=\frac{1+t}{2}\boldsymbol{v}_{j+1}+\frac{1-t}{2}\boldsymbol{v}_{j+2},\quad t\in[-1,1],$$

the affine parametrization of the edge $s_{j}$. Let $\omega\in L^{1}([-1,1])$ be a general pdf on $[-1,1]$. We assume further that $\omega$ admits finite moments

$$\mu_{m,\omega}=\int_{-1}^{1}t^{m}\omega(t)dt,\quad m\in\mathbb{N}.$$ (72)

3.1 Edge functionals and orthogonal polynomials

Given such a density $\omega$, we define enriched edge functionals in analogy with the previous constructions by

$$\mathcal{I}^{\omega}_{j}(f) =\int_{-1}^{1}f(\gamma_{j}(t))\omega(t)dt,$$ (73)

$$\mathcal{L}^{\omega}_{j}(f) =\int_{-1}^{1}q(t)f(\gamma_{j}(t))\omega(t)dt,$$ (74)

where $q\in\mathbb{P}_{r}([-1,1])$, $r\geq 2$, is a polynomial satisfying the orthogonality conditions

$$\int_{-1}^{1}q(t)\omega(t)dt=0,\quad\int_{-1}^{1}tq(t)\omega(t)dt=0,$$ (75)

together with the non-degeneracy condition on the second moment

$$\mu_{2,\omega}:=\int_{-1}^{1}t^{2}q(t)\omega(t)dt\neq 0.$$ (76)

These assumptions ensure that the $q$-weighted edge functionals (74) vanish on affine polynomials while detecting quadratic contributions.

3.2 Unisolvency and reconstruction

We then define the enriched triple

$$\mathcal{H}^{\omega}=\left(T,\mathbb{P}_{2}(T),\Sigma^{\omega}(T)\right),$$

where

$$\Sigma^{\omega}(T)=\left\{\mathcal{I}^{\omega}_{j},\mathcal{L}^{\omega}_{j}\,:\,j=1,2,3\right\}.$$

The following result shows that the structural properties observed for the specific families extend to any choice of probability density $\omega$.

Theorem 3.1

The triple $\mathcal{H}^{\omega}$ is unisolvent.

Proof(Sketch of proof)

The argument follows the same lines as in the proofs of Theorems 2.3 and 2.19. Indeed, let $p_{2}\in\mathbb{P}_{2}(T)$ such that

$$\mathcal{I}_{j}^{\omega}\left(p_{2}\right)=\mathcal{L}_{j}^{\omega}\left(p_{2}\right)=0,\quad j=1,2,3.$$

Along each edge $s_{j}$, the trace of $p_{2}$ is quadratic. By construction, $\mathcal{L}^{\omega}_{j}$ annihilates constant and linear terms, while detecting the quadratic coefficient. Hence, vanishing of $\mathcal{L}^{\omega}_{j}\left(p_{2}\right)$ implies that all edge traces of $p_{2}$ are affine. This forces $p_{2}\in\mathbb{P}_{1}(T)$. The remaining conditions $\mathcal{I}^{\omega}_{j}\left(p_{2}\right)=0$ then yield a homogeneous system for the linear coefficients, which admits only the trivial solution. This completes the proof.

This generalization shows that the enriched quadratic reconstruction framework is not restricted to the two-parameter families considered earlier, but in fact applies to any choice of probability density on the edges. The specific distributions $\widetilde{k}_{\sigma,\mu}$ and $\widetilde{g}_{\sigma,\mu}$ therefore appear as special cases within a broad and flexible class, highlighting the robustness and universality of the proposed approach.

Theorem 3.2

If the probability density function $\omega$ is even, then the basis functions associated with the enriched triple $\mathcal{H}^{\omega}$ take the form

$$\varphi_{\omega,i} =1-2\lambda_{i},\quad i=1,2,3,$$

$$\psi_{\omega,i} =-A_{\omega}\varphi_{\omega,i}+\frac{2}{\left\|{q}\right\|^{2}_{2,{\omega}}}\left(-\lambda_{i}^{2}+\lambda^{2}_{i+1}+\lambda^{2}_{i+2}\right),\quad i=1,2,3,$$

where

$$\left\|{q}\right\|^{2}_{2,{\omega}}=\int_{-1}^{1}q^{2}(t)\omega(t)dt,$$

$$A_{\omega}=\frac{1+\mu_{2,\omega}}{\left\|{q}\right\|^{2}_{2,{\omega}}}.$$

3.3 Practical constructions of $q$

When $q$ is chosen to be a quadratic polynomial, it is possible to construct it explicitly from the moments $\mu_{m,\omega}$. We outline below three practical approaches.

1. Canonical quadratic form. Seek $q(t)=t^{2}-a-bt$ subject to the orthogonality conditions (75). Using the notation (72), the orthogonality conditions yield the linear system

   $$\begin{cases}\mu_{2,\omega}-a-b\mu_{1,\omega}=0,\\ \mu_{3,\omega}-a\mu_{1,\omega}-b\mu_{2,\omega}=0,\end{cases}\quad\Longrightarrow\quad b=\frac{\mu_{1,\omega}\mu_{2,\omega}-\mu_{3,\omega}}{\mu_{1,\omega}^{2}-\mu_{2,\omega}},\quad a=\mu_{2,\omega}-b\mu_{1,\omega},$$

   provided $\mu_{2,\omega}\neq\mu_{1,\omega}^{2}$ (nonzero variance). Optionally normalize $q$ by

   $$\widehat{q}=\frac{q}{\nu}$$

   where

   $$\nu=\mu_{4,\omega}-a\mu_{2,\omega}-b\mu_{3,\omega}.$$

2. Orthogonal polynomial. Construct the family of orthogonal polynomials $\left\{\pi_{k}\right\}_{k\geq 0}$ with respect to the inner product

   $$\left\langle f,g\right\rangle_{\omega}=\int_{-1}^{1}f(t)g(t)\omega(t)dt.$$

   Taking $q=\pi_{2}$ (or its normalized version) ensures that conditions (75)–(76) are satisfied by construction milovanovic1997orthogonal; mastroianni2008interpolation.

3. Higher degree. Any polynomial $q$ of degree $\geq 2$ fulfilling (75)–(76) is admissible. The quadratic choice is however preferred, since it requires only low-order moments and ensures stability.

Remark 13

If the chosen density $\omega$ is symmetric with respect to the origin, then all odd moments vanish, i.e.,

$$\mu_{2k+1,\omega}=0,\quad k\in\mathbb{N}.$$

In this case, the expressions in (A) simplify considerably: one obtains $b=0$ and hence

$$q(t)=t^{2}-\mu_{2,\omega}.$$

This situation arises in particular for both the first and the second family of densities introduced in Section 2, which are symmetric and therefore lead to especially simple quadratic forms.

3.4 Numerical procedure

Before presenting the numerical results, we summarize the computational workflow adopted in our experiments. The procedure consists of two phases:

• Offline parameter tuning: optimal density parameters $\left(\mu^{\star},\sigma^{\star}\right)$ are selected by Algorithm 1, applied to a small validation set of functions and a sequence of meshes.

• Reconstruction and error evaluation: once the optimal parameters are fixed, they are used uniformly across all tests to compare the classical histopolation scheme $\mathcal{CH}$ with the enriched scheme $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$.

The pseudocode of the complete workflow is reported in Algorithm 2.

Algorithm 2 Numerical testing workflow

1:Test functions $\left\{f_{r}\right\}_{r=1}^{R}$; mesh sequence $\left\{\mathcal{T}_{n}\right\}_{n=0}^{N}$; parameter grid $\mathcal{P}$

2:$L^{1}$ errors for $\mathcal{CH}$ and $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$

3:$\left(\mu^{\star},\sigma^{\star}\right)\leftarrow$ parameter tuning via Algorithm 1

4:for each test function $f_{j}$ do

5:  for each mesh $\mathcal{T}_{n}$ do

6:   $u^{\mathrm{CH}}\leftarrow\pi^{\mathrm{CH}}\left[f_{j};\mathcal{T}_{n}\right]$

7:   $u^{\mathrm{enr}}\leftarrow\pi_{1,\sigma,\mu}^{\mathrm{enr}}\left[f_{j};\mathcal{T}_{n}\right]$

8:   Compute $E^{\mathrm{CH}}_{j,n}=\left\|f_{j}-u^{\mathrm{CH}}\right\|_{L^{1}(\Omega;\mathcal{T}_{n})}$

9:   Compute $E^{\mathrm{enr}}_{j,n}=\left\|f_{j}-u^{\mathrm{enr}}\right\|_{L^{1}(\Omega;\mathcal{T}_{n})}$

10:  end for

11:end for

12:return Error arrays $\left\{E^{\mathrm{CH}}_{j,n},E^{\mathrm{enr}}_{j,n}\right\}$

This explicit workflow highlights how the parameter tuning step integrates seamlessly into the testing phase: after a one-time offline calibration, the enriched scheme is applied to new test functions and meshes in a fully specified and reproducible way.

4 Numerical tests

In this section, we present a series of numerical experiments to assess the performance and accuracy of the proposed enriched histopolation methods. All tests are carried out on the square domain $\Omega=[-1,1]^{2}$, using a collection of benchmark functions that exhibit a variety of analytic behaviors, including smooth, oscillatory, and localized features. The test functions are defined as follows

$$f_{1}(x,y) = \sqrt{x^{2}+y^{2}},\ \ f_{2}(x,y)=e^{-4\left(x^{2}+y^{2}\right)}\sin\left(\pi(x+y)\right),$$

$$f_{3}(x,y) = \sin(2\pi x)\sin(2\pi y),\ \,f_{4}(x,y)=\sin\left(4\pi(x+y)\right),\ \,f_{5}(x,y)=\frac{1}{25(x^{2}+y^{2})+1},$$

$$f_{6}(x,y) = 0.75\exp\biggl(-\frac{(9(x+1)/2-2)^{2}}{4}-\frac{(9(y+1)/2-2)^{2}}{4}\biggr)$$

$$+0.75\exp\biggl(-\frac{(9(x+1)/2+1)^{2}}{49}-\frac{(9(y+1)/2+1)}{10}\biggr)$$

$$+0.5\exp\biggl(-\frac{(9(x+1)/2-7)^{2}}{4}-\frac{(9(y+1)/2-3)^{2}}{4}\biggr)$$

$$-0.2\exp\biggl(-(9(x+1)/2-4)^{2}-(9(y+1)/2-7)^{2}\biggr).$$

The function $f_{6}$ is the well-known Franke function, widely used as a benchmark for approximation methods franke1982scattered. For the spatial discretization, we consider a sequence of regular Friedrichs–Keller triangulations Knabner, denoted by

$$\mathcal{T}_{n}=\left\{T_{i}:i=1,\dots,2(n+1)^{2}\right\},$$

where each $\mathcal{T}_{n}$ consists of $2(n+1)^{2}$ triangles; see Fig. 3.

Figure 3: Successive regular Friedrichs–Keller triangulations $\mathcal{T}_{n}$, with $n=20$ (top left), $n=30$ (top right), $n=40$ (bottom left), and $n=50$ (bottom right).

To evaluate the quality of the reconstruction, we compute the error in the $L^{1}$ norm produced by the classical linear histopolation method $\mathcal{CH}$ defined in (17) and compare it with the error obtained using the enriched histopolation method $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$ introduced in (22), with the optimal parameters $\mu^{\star}=2$ and $\sigma^{\star}=1$. These parameters have been determined by means of the parameter tuning Algorithm 1. All experiments have been carried out according to Algorithm 2, implemented in Matlab. The integrals needed for assembling the reconstruction operators and for evaluating the $L^{1}$ errors were computed on each triangular element using high-order quadrature rules, with a sufficiently dense set of nodes to ensure accurate numerical integration and to rule out any artifacts due to quadrature errors. The $L^{1}$ norm is chosen as it handles discontinuities effectively without requiring limiting procedures (cf. dell2025truncated). The results of this comparison are reported in Figures 4, 5 and 6. A clear improvement of the enriched method over the standard scheme is observed across all test functions and mesh refinements. In particular, $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$ achieves significantly smaller $L^{1}$ errors, demonstrating its superior capability in reproducing local oscillations, sharp gradients, and singular features. Moreover, the enriched approach demonstrates an accelerated error decay under mesh refinement, thereby confirming that the additional degrees of freedom provided by the probability density weights yield a genuine enhancement in approximation power.

Figure 4: Semi-log plot of the $L^{1}$ approximation error for $f_{1}$ (left) and $f_{2}$ (right) obtained with the classical histopolation method $\mathcal{CH}$ (blue) and the weighted quadratic enriched method $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$ (red), as the number of triangles in the Friedrichs–Keller triangulations increases.

Figure 5: Semi-log plot of the $L^{1}$ approximation error for $f_{3}$ (left) and $f_{4}$ (right) obtained with the classical histopolation method $\mathcal{CH}$ (blue) and the weighted quadratic enriched method $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$ (red), as the number of triangles in the Friedrichs–Keller triangulations increases.

Figure 6: Semi-log plot of the $L^{1}$ approximation error for $f_{5}$ (left) and $f_{6}$ (right) obtained with the classical histopolation method $\mathcal{CH}$ (blue) and the weighted quadratic enriched method $\mathcal{H}^{\mathrm{enr}}_{1,\sigma^{\star},\mu^{\star}}$ (red), as the number of triangles in the Friedrichs–Keller triangulations increases.

5 Conclusions and Future works

In this work we have introduced two new two-parameter families of generalized truncated normal distributions and exploited them to design enriched local histopolation schemes based on unisolvent triples. These distributions serve as flexible edge weights and provide additional degrees of freedom, which translate into improved accuracy and adaptability of the reconstruction operators. The resulting quadratic operators preserve the simplicity of the classical linear histopolation method while significantly enhancing their approximation capabilities, especially for functions with oscillations, steep gradients, or localized singularities.

From a theoretical point of view, we have established unisolvency and derived explicit closed-form basis functions, thus ensuring that the proposed schemes are mathematically well-posed and computationally viable. We have also proposed an algorithm for the optimal selection of the distribution parameters, which further increases the robustness and adaptivity of the approach. Numerical experiments have confirmed the theoretical findings, showing systematic improvements over the linear nonconforming histopolation approach. Moreover, as shown in Section 3, the same reasoning extends to any bivariate weight function whose restriction to the edges of the triangle defines a valid probability density. This highlights the generality and robustness of the proposed framework, which is not limited to the two specific families presented in this paper.

Several directions for future research naturally arise from this study. First, extending the framework to three-dimensional meshes would broaden its applicability to volumetric data reconstruction and tomographic imaging. Second, adaptive strategies for selecting the distribution parameters could be developed, so as to optimize accuracy depending on local features of the target function. Third, connections with other classes of orthogonal polynomials and probability densities could be investigated, possibly leading to new families of enriched operators. Finally, applications to real-world problems in imaging, computer vision, and numerical approximation remain a promising avenue for further exploration.

Conflict of interest

Not Applicable.