SBProfile::draw : Flux/imscale scaling questionable when full jacobian is provided. about galsim HOT 4 CLOSED

Yes, you're right about this. It should be flux_ratio *= std::abs(invdet * imscale) to get an extra factor of imscale. Thanks for catching this!

Fortunately, it's a pretty severe edge case. The normal paths of drawing in GalSim usually end up at this point with imscale=1, so the error doesn't cause any problems (which is probably why no one ever noticed it until now).

I think the only way to hit this case is to:

1. Use obj.drawReal() directly rather than the usual obj.drawImage().
2. The object being drawn needs to be transformed with a non-diagonal transformation (either shear or rotation).
3. The image being drawn onto needs to have a non-unit pixel scale (but note that it's not allowed to have a non-PixelScale wcs).

So, it's definitely worth fixing of course, but it's notable that there is no code in our current test suite that hits this situation, never mind has a test that should have noticed the flux being wrong. So I suspect it most likely hasn't ever caused an error in anyone's images. Thank goodness. :)

I'll make a PR with the fix shortly.

from galsim.

I'm glad to have found something valuable.
If I can add a question, here is one: does in the code

            double mA = jac[0];
            double mB = jac[1];
            double mC = jac[2];
            double mD = jac[3];

means that the jacobian matrix is set as (convention)

\begin{matrix}
mA & mB \\
mC & mD
\end{matrix}

or

\begin{matrix}
mA & mC \\
mB & mD
\end{matrix}

(it is a matter of transpose or not)

In fact I'm a bit trouble by this in GalSim/SBProfile.cpp/draw()

            double mA = jac[0];
            double mB = jac[1];
            double mC = jac[2];
            double mD = jac[3];
            double det = mA*mD - mB*mC;
            double invdet = 1./det;
            double new_x0 = (mD*x0 - mB*y0) * invdet;
            double new_y0 = (-mC*x0 + mA*y0) * invdet;

which would make sense if

J = \begin{matrix}
mA & mB \\
mC & mD
\end{matrix} 

=>

J^{-1} = (mA mD - mB mC)^{-1} \begin{matrix}
mD & -mB \\
-mC & mA
\end{matrix} 

as then

\begin{matrix}
new_x0 \\ new_y0
\end{matrix}    
=  J^{-1} \times  \begin{matrix}
x0 \\ y0
\end{matrix}   

But looking at drawK method I read

            double mA = jac[0];
            double mB = jac[1];
            double mC = jac[2];
            double mD = jac[3];
            double new_x0 = mA*x0 + mC*y0;
            double new_y0 = mB*x0 + mD*y0;

which would mean that jac is this routine is like the transpose of the previous one.
Have I missed something?
Thanks

from galsim.

The effect of the Jacobian is transposed in Fourier space. So that's correct. J = [[ a, b], [c, d]]. But in Fourier space, you apply J^T.

from galsim.

Ok. +1

from galsim.