I have a discrete image of size $2^N \times 2^N$ which I would like to iteratively downsample to produce a pyramid with image sizes $2^k\times 2^k, \, k=0,\ldots,N$. That is, each subsequent image will have half the resolution along each axis. So I would like to better understand how this is modeled mathematically. My current understanding of the ideal downsampling of a continuous image $f$ is that one wants to compute the discrete image $\boldsymbol{g} \in \mathbb{R}^{2^{N-1}\times 2^{N-1}}$ such that:
$$g_{ij}=\int_{x_{2i}}^{x_{2(i+1)}}\int_{y_{2j}}^{y_{2(j+1)}}f(x,y)\,dy\,dx,$$
where pixel $(i,j)$ on the coarser level corresponds to some rectangle $[x_{2i},x_{2(i+1)}]\times[y_{2j},y_{2(j+1)}]$. The main trouble is that I am not given the continuous function, but instead some discrete data $\boldsymbol{h} \in \mathbb{R}^{2^N \times 2^N}$. To reconstruct a continuous function using a linear model I would have something of the form:
$$f(\boldsymbol{z}) = \sum_{i=1}^{2^N}\sum_{j=1}^{2^N}c_{ij}\phi_{ij}(\boldsymbol{z}),$$
where $c_{ij}$ are some coefficients, and $\phi_{ij}$ are some predetermined basis functions. A standard constraint used to derive $c_{ij}$ is the interpolation requirement: $h_{ij} = f(\boldsymbol{z}_{ij})$, where for example the samples are taken in the middle of each pixel: $\boldsymbol{z}_{ij} = \frac{1}{2}(x_i+x_{i+1}, y_j+y_{j+1})$. This is a linear system, and provided that the basis functions are linearly independent, then there is a unique set of coefficient $\boldsymbol{c} \in \mathbb{R}^{2^N\times 2^N}$ corresponding to $\boldsymbol{h}\in \mathbb{R}^{2^N\times 2^N}$.
Correct me if I am wrong, but I believe that if $\phi_{ij}$ are box functions, then computing the integral for $g$ is equivalent to using bilinear basis functions interpolation and then sampling at the midpoint. I believe a similar thing holds for higher degrees, where if $f$ is assumed to be formed by polynomials of degree $n$ then evaluating the integral would be equivalent to sampling a $n+1$ degree interpolated version of $\boldsymbol{h}$. The box filter version seems to be special in that it requires the basis functions to be Dirac deltas with a scaling factor of $\frac{1}{4}$ and $\boldsymbol{c} = \boldsymbol{h}$.
And this is the point where I am not certain about things anymore. I have seen the $sinc$ function described as an ideal interpolator in the context of the Nyquist-Shannon sampling theorem, and I have also seen it discussed in the context of downsampling. How does such a $sinc$ downsampling relate to the above mathematical modelling of the problem? I am not quite sure how this downsampling is achieved in the first place or what it means when a discrete set of samples are considered. I understand that $\mathcal{F}[f * sinc] = \mathcal{F}[f] \cdot box$ can be used to cut off frequencies above the Nyquist frequency when $f$ is a continuous function. On the other hand I have seen people just put zeroes as DFT coefficients (notably this is the setting of a discrete signal). Is the idea that they consider $f$ a sum of Dirac deltas weighted by $h_{ij}$, and then consider the convolution with an appropriate $sinc$ to cut-off frequencies that are too high? This would essentially result in considering $sinc$ as basis functions. Where is the integration step that I have at the very beginning however? Am I misunderstanding how $sinc$ downsampling is supposed to work? Also there seem to be some concerns with ringing artifacts arising from the above? What is the mathematical reason for this?
