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

Question

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 )