Is it possible to combine decimation and low pass filtering in one step? Not necessarily only for images but also for general signals.
Yes, that's what people usually do when they implement downsampling: since of the output of the anti-aliasing filter, you throw away N-1 samples, why even calculate these?
The trick is to decompose your filter into polyphase components, which enables you to run the resulting filter operation only once per output of the downsampling, instead of once per input. There's plenty of reference implementations - from GNU Radio's decimating FIR filters, to rescalers in image processing hardware.
Think of it this way:
The trick is to take your original filter $[h_0, h_1, h_2, h_3, \ldots, h_N, h_{N+1}, h_{N+2},\ldots,h_{2N}, h_{2N+1}, \ldots]$ and just split it up into filters where there's only one non-zero entry every $N$ coefficients. Choose the non-zero-value positions so that the first polyphase component filter gets $h_0, h_N, h_{2N},\ldots $, the second gets $h_1, h_{N+1}, h_{2N+1},\ldots$ and so on.
Add up the result of these filters, when you feed in the same input, to "undo" the splitting. This doesn't change anything, it's the same filter, just split into $N$ filters with many zeros in them, but with the non-zero elements in different positions.
After the addition, you decimate by $N$. Ok, you can do that before the addition, so now you have one input stream, fed into $N$ subfilters, each with a lot of zeros in them, each followed by a decimation by $N$.
Now you have a special kind of filter that only had every Nth filter tap occupied, so the first subfilter's coefficient vector is $[h_0, 0, \ldots, 0, h_N, 0, \ldots, 0, h_{2N}, 0 \ldots]$, and you'd decimate by $N$ afterwards, you could as well just swap decimation and filter, and just use the filter $[h_0,h_N,h_{2N},\ldots]$. The two things are identical in effect; this is called Noble Identity.
So, we can "pull the decimation up front" for that filter. You can, in fact, do that to all subfilters (you will have to add delay so that it works out mathematically for the non-zero-phase polyphase components, but the idea doesn't change. You have one input stream, going into $N$ different delays, decimate-by-$N$ decimators, subfilters, and a summation.
As it happens, this means that only one "branch" at a time actually gets input per input cycle.