Making a low pass filter for irregular samples
I realize that this is not a filter, but I do not believe any filter will give you a reasonable result given the strange sample locations. This seems to be borne out by several references on the web. This one for example and this PDF one.

I suggest something something along these lines, using projections onto convex sets (POCS).

Algorithm Descroption

Assume that the known red pixel values $R_i$ are at $(x_i,y_i)$ so that your start red image is:
$$
I_{\tt red}(x,y) = \sum_{i=1}^N R_i \delta(x-x_i, y-y_i)
$$

Set $k=0$, and initialize your estimate of the fully sampled red image with $I_{\tt red}(x,y)$:
$$
\hat{I}^k_{\tt red}(x,y) = I_{\tt red}(x,y).
$$

Take the 2D FFT of $\hat{I}^k_{\tt red}$:
$$
F_{\tt red}(n,m) = {\tt FFT}\{\hat{I}^k_{\tt red}(x,y)\}
$$
Because the red channel is slowly varying, remove all frequencies above a cutoff:
$$
F^{lp}_{\tt red}(n,m) = F_{\tt red}(n,m) L(n,m)
$$
where $L$ is a low pass mask.

Form the inverse FFT to get $\hat{I}^{k+1}_{\tt red}(x,y)$
$$
\hat{I}^{k+1}_{\tt red}(x,y) = {\tt FFT}^{-1}\{F^{lp}_{\tt red}(n,m)\}
$$

Set the known pixels of $\hat{I}^{k+1}_{\tt red}(x,y)$ from $I_{\tt red}(x,y)$.

If the difference between $\hat{I}^{k+1}_{\tt red}(x,y)$ and $\hat{I}^{k}_{\tt red}(x,y)$ is small enough, stop, otherwise reiterate with $k=k+1$.

Results

An example run with just a $6 \times 6$ matrix is below. The left image shows the original pixels. The right image shows the POCS version.  The red dots show the original pixel locations in both images.



Different Initial Conditions

I tried setting the initial image to

Ired &lt;- array(100000,c(ImageSize, ImageSize))


and 

Ired &lt;- array(128,c(ImageSize, ImageSize))


and the results are below, so I do not believe the initial values for the unknown pixels matter.

The plot scales the greyscale values, which is why the 100,000 example looks wrong in the left-hand plot.







R Code Below

# 32018
# Non-Bayer 2D interpolation

# Model pixels as:
# P(x,y) = (a_0 + a_1 * x + a_2 * x^2) * (b_0 + b_1 * y + b_2 * y^2)
#        = [a_0 a_1 a_2 b_0 b_1 b_2] * [1 x x^2 1 y y^2]^T 
# And so we need to find the a_i and b_i from specific pixel and coordinate values.
#
# install.packages('ptw')

require(ptw)

pad &lt;- function(x,N)
{
  p1 &lt;- ptw::padzeros(x,N,side="right")
  p2 &lt;- t(ptw::padzeros(t(p1),N,side="right"))

  return(p2)
}

Npoints &lt;- 8
xy &lt;- array(c(c(1,1), c(1,5), c(2,3), c(3,6), c(4,2), c(4,4), c(5,6), c(6,3)), c(2,Npoints))
xy_norm &lt;- (xy - min(xy))/(max(xy) - min(xy))
Pxy &lt;- array(runif(N, 0, 255), c(1,Npoints))

ImageSize &lt;- max(xy)

Ired &lt;- array(0,c(ImageSize, ImageSize))

for (idx in 1:Npoints)
{
  Ired[xy[1,idx], xy[2,idx]] = Pxy[idx]
}

Ihatred &lt;-array(0,c(ImageSize,ImageSize))
Ihatred2 &lt;- Ired
count &lt;- 0

while (sum(abs(Ihatred - Ihatred2)) &gt; 1)
{
  Ihatred &lt;- Ihatred2

  Fred &lt;- fft(pad(Ihatred, 2*ImageSize))

  Fredsize &lt;- dim(Fred)[1]
  cutoff1 &lt;- Fredsize / 8
  cutoff2 &lt;- 7 * Fredsize / 8

  Fred[4:16,] &lt;- 0
  Fred[,4:16] &lt;- 0

  IFred &lt;- fft(Fred, inverse=TRUE) / Fredsize / Fredsize

  Ihatred2 &lt;- Re(IFred[1:ImageSize, 1:ImageSize])

  for (idx in 1:Npoints)
  {
    Ihatred2[xy[1,idx], xy[2,idx]] = Pxy[idx]
  }

  count &lt;- count + 1  
  if (count &gt; 100)
  {
    break
  }  
}


par(mfrow=c(1,2))
image(Ired, col= grey(seq(0, 1, length = 256)))
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}
title('Original with only known red pixels')
image(Ihatred2, col= grey(seq(0, 1, length = 256)))
title('POCS solution')
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}