Iterative gaussian convolution and decimation when the decimation/convolution ratios are not divisible by two?

I have a bit of a mathematics problem.

I have an image $$I_1$$ that I have previously blurred (from an image $$I_0$$) with a Gaussian of standard deviation $$\sigma_1$$.

I would like to blur image $$I_1$$ with a Gaussian of standard deviation $$\sigma_X$$ to produce a new image, image $$I_2$$. The result, image $$I_2$$, should be as if $$I_0$$ has been blurred with a Gaussian of standard deviation $$\sigma_2$$.

This problem is not difficult, we know that the convolution of two Gaussian filters is the third filter. I can compute the $$\sigma_X$$ by: $$\sigma_X = \sqrt{ \sigma_2^2 - \sigma_1^2 }$$ (this works because application of a Gaussian of $$\sigma_1$$ followed by application of Gaussian of $$\sigma_2$$, is as if a single Gaussian of $$\sigma_3$$ had been applied, where $$\sigma_3^2 = \sigma_1^2 + \sigma_2^2$$).

I can directly apply a Gaussian with $$\sigma_X$$ to image $$I_1$$ to achieve my desired image $$I_2$$.

Furthermore, I would like to decimate the image $$I_2$$ by a factor of $$N$$ (e.g., $$0.5$$ is decimation by half, $$0.25$$ is by fourth) after the blur. Keep in mind that neither of $$1/N$$ nor $$\log_2(\sigma_2/\sigma_1)$$ must be whole numbers, though they will always be positive.

Here is the problem. In the case where $$\sigma_X$$ is very large, the process becomes computationally expensive. Since I will be decimating the image anyways, I would like to apply the method often used in (dyadic) Gaussian pyramids, where one applies a blur followed by a decimation step, and then repeats this. However, since $$\log_2(\sigma_2/\sigma_1)$$ and $$1/N$$ are not guaranteed to be whole numbers, it is not straightforward to break the transform down into smaller blur-decimate steps so that I never have to apply a very large Gaussian convolution.

I have a feeling that this is possible.

Currently, I have:

    double sigmaratio = sigma2 / sigma1;
double xsigma = sqrt( sigma2*sigma2 - sigma1*sigma1);
double N = decimation_ratio //(e.g. 0.5 means decimate by half, 0.25 is by fourth, etc.)
double sig_halves = log2( sigmaratio );
double N_halves = -1 * log2(N); //will be =1 if we halve it once, =2 if twice, etc.

//Taking floor to get the whole part...
int intsighalves = (int) sig_halves;
int intNhalves = (int) N_halves;
int num = intsighalves;
if( intNhalves < num )
{
num = intNhalves;
}
if( num >= 1 )
{
double Ndiff = N_halves - (double)num;
double sigdiff = sig_halves - (double)num;
//I need to apply num dyadic steps with some sigmaXprime at each step (getting 2x as large each time) and decimation by half at each step.
}


My question is: how do I choose the sigmaXprime at each dyadic step for the whole parts, and when/how do I choose the partial sigmaXprimepartial and the partial decimation in order to get an exact implementation of the image $$I_1$$ convolved with Gaussian $$\sigma_X$$ followed by decimation by $$N$$?

Edit: Just a note, this may not be possible due to loss of information at each decimation step, but I think it should be possible since in the case of dyadic Gaussian pyramids, application of a single Gaussian / decimation step with the correct parameters is exactly equal to application of iterated Gaussian / decimation steps (see my answer to https://dsp.stackexchange.com/questions/667/image-pyramid-without-decimation/55654#55654 )

• It is not clear what you want. Could you please try explaining it clearer? – Royi Mar 2 at 18:14
• Thank you for the feedback, I will try. I need to figure out how to write latex here or else it will never be clear I think. – rveale Mar 2 at 18:21
• For Math just write as you'd write in LaTeX: For In Line Math: $<In Line Math Code>$, Display Math: $$<Display Math Code>$$. – Royi Mar 2 at 20:58
• Use FFT for the convolution. Decimation is by an integer factor. Resampling is changing the sample frequency, this is achieved by a convolution with $\frac{\sin(n)}{n}$ similar to that for the upsampling. Since we want non-circular convolution we need zero padding, and $\frac{\sin x}{x}$ is replaced by $\frac{\sin x}{x} w(x)$ for some window $w$. – reuns Mar 4 at 20:26