Pure averaging isn't the best way to go, that's basically a box filter. You can use the interpolation kernels, but you need to widen them proportionally to the amount you're downscaling. Instead of your bicubic working with a 4x4 area, for 3x downscaling it needs 12x12. Ideally for ½ scale you should be using 8x8 already.
A box filter for scaling down cubic and linear is much better than using impulse filters with the interpolated values at the floating point coordinates in question. If someone wants to improve the code in GEGL (used for rotations, scaling, perspective transforms, live on image warping as well as other transforms in the development version of GIMP) the relevant place to put it in is here:
For the linear case a "tent-filter" instead of a box filter would be really correct, a similar box filter function exist in cubic. Both of these receive a reverse jacobian matrix indicating the shape and extent of surrounding area in source image to sample from.
To eliminate moire artifacts you need to remove all frequencies above the Nyquist limit. In resizing applications there are two Nyquist limits, one for the input and one for the output; you need to filter for whichever is smallest. When upsizing the input limit is always smallest, so the filters can be constant. When downsizing the output limit is smallest so that's the one you need to tune your filter for. That's why I suggest widening the interpolation formula when downsizing.
I've been meaning for years to make a blog post on this subject. I don't think many people realize that an interpolation formula is also a filter formula, and that it can be manipulated and analyzed as such.
I have my own implementation of the classic filters that I use for my resizing tasks. It works in linear gamma space as suggested in the article. I've implemented lots of different algorithms, and I've settled on Lanczos-5 as the best overall compromise. One interesting observation, the Catmull-Rom bicubic interpolation is nearly indistinguishable from Lanczos-2.
