First, the horizontal and vertical summation passes are separated. Then for each pass, it's using a DDA  to generate the source points for each output point
(yes, it's mapping the points in reverse direction: for each output position, find the fractional value of the source point that would be rescaled to this location).
Then you use the fractional part of the position to scale the input sample values. It's indeed an integral (/summed) rescaling because of that. Could use smarter
interpolation (similar to trapezoid method for integrals), but it's not really improving the result much. Midpoint is fine, and computation is easier (just one mult).
Hope the following schema helps understanding: