When computing double integrals with change of variables, one needs to apply the scaling factor called the Jacobian.
${\displaystyle \left|J\right| = \left|\frac{\partial \left(x,y\right)}{\partial \left(u,v\right)}\right| = \left|\frac{1}{\frac{\partial \left(u,v\right)}{\partial \left(x,y\right)}}\right| }$
But I am not sure which one should be used. When I look at examples on this it seems that one would apply the second reciprocal form of the Jacobian when your initial space has a region that is non-rectangular and in the new space, after your change of variables you end up with a rectangular region.
Hope someone will clarify this.
Here is an example image:
Here is where I feel there is confusion, based on this image I have here:
For what follows below, one can surely weaken several regularity assumptions, but I just want to convey the formula. Suppose $f:\Bbb{R}^2\to\Bbb{R}$ is given and you're given a smooth diffeomorphism $\phi:\Bbb{R}^2\to\Bbb{R}^2$. To set the notation, let us denote the coordinates on the domain of $\phi$ by $(u,v)$ and the coordinates on the target of $\phi$ by $(x,y)$. So, $\phi$ is like a mapping of "$(u,v)$ coordinates to $(x,y)$ coordinates". Then, given a set $A\subset\Bbb{R}^2$ in the "$(u,v)$ coordinates", we have \begin{align} \int_{\phi(A)}f(x,y)\,dx\,dy&=\int_{A}f(\phi(u,v))\cdot\left|\det D\phi_{(u,v)}\right|\,du\,dv\\ &\equiv \int_{A}f(\phi(u,v))\cdot\left|\det \frac{\partial (x,y)}{\partial(u,v)}\right|\,du\,dv, \end{align} So, if you want to really stretch the limits of logic and notation, it's like $dx\,dy =\left|\det\frac{\partial(x,y)}{\partial(u,v)}\right|\,\,du\,dv$, so the $(u,v)$ appears once "on top" in the form of $du\,dv$ and once "on the bottom" as $\partial(u,v)$, so they "cancel".
Also, just to clarify, in this theorem, there's no such thing as "$(u,v)$ are curved coordinates" while "$(x,y)$ are cartesian coordinates". We're just using certain symbols to conveniently express a given formula.
For example, in $\Bbb{R}^n$, if I have a function $f:\Bbb{R}^n\to\Bbb{R}$ and a diffeomorphism $\phi:\Bbb{R}^n\to\Bbb{R}^n$, I can label the coordinates however I wish. So, I can say I would like to label coordinates on the domain of $\phi$ as $\xi=(\xi^1,\dots, \xi^n)$ and on the target space of $\phi$ as $\zeta=(\zeta^1,\dots, \zeta^n)$. Then, for $A\subset\Bbb{R}^n$ the formula would look like \begin{align} \int_{\phi(A)}f(\zeta^1,\dots \zeta^n)\,d\zeta^1\cdots\,d\zeta^n&=\int_Af(\phi(\xi))\cdot\left|\det \frac{\partial (\zeta^1,\cdots, \zeta^n)}{\partial (\xi^1,\cdots, \xi^n)}\right|\,d\xi^1\cdots\,d\xi^n \end{align} Of course, the most concise way of writing the formula is to not introduce arbitrary letters like $\xi,\zeta,x,y,u,v$ etc. All we need are the function $f$, the change of variable mapping $\phi$ and a set $A\subset\text{domain}(\phi)$. Then, the formula reads \begin{align} \int_{\phi(A)}f&=\int_{A}f\circ \phi \cdot |\det D\phi|. \end{align}
Edit:
As a general rule, I would suggest you get out of the habit of thinking in terms of symbols like $x,y$ or $u,v$ or $r,\theta$ etc. The choice of letters does NOT dictate mathematics. For example, I could say that $(\mu,@)\mapsto (\mu\cos@,\mu\sin@)$ is a polar coordinate mapping, and if I really wanted to stretch notation, I can even say that $(x,y)\mapsto (x\cos y, x\sin y)$ defines a polar coordinate mapping, where now $x$ is radial distance and $y$ denotes the angle. So, you should be flexible with your letters, and focus on what the actual statement being made is.
When attempting to solve questions involving the change of variables for integration, you always have to identify three things and then pattern match to get the right formula
Again, don't get hung up on the letters I used $f,\phi,A$ (I just happen to like these letters. If you look in another place, definitely some different notation is used).
Anyway, in your example, we consider the function $f:\Bbb{R}^2\to\Bbb{R}$ defined as $f(x,y)=x^2+y^2$. Next, we consider the mapping $\phi:\Bbb{R}^2\to\Bbb{R}^2$ defined as $\phi(x,y)=(x^2-y^2,2xy)$. Lastly, we consider the rectangular region $R'=[1,9]\times [4,8]$. Then, the region $R$ in your question is nothing but the preimage set $\phi^{-1}(R')$. So, by pattern matching the definitions and applying the change-of-variables theorem, we get \begin{align} \int_{R}(x^2+y^2)\,dx\,dy&=\int_{\phi^{-1}(R')}f(x,y)\,dx\,dy\\ &=\int_{R'}(f\circ \phi^{-1})(u,v)\cdot |\det D(\phi^{-1})_{(u,v)}|\,du\,dv\tag{$*$} \end{align} So, we have to evaluate the integral on the RHS. As it is currently written, it looks very daunting, because of the appearance of the inverse function $\phi^{-1}$. We have an explicit formula for $\phi$, but as you know, inverting things is always a nightmare, and right now we have to invert things and then also calculate the derivative of the inverse. What we can do to simplify our calculation is invoke the inverse function theorem which tells us how to calculate derivatives of inverse functions. We have \begin{align} \int_{R'}(f\circ \phi^{-1})(u,v)\cdot |\det D(\phi^{-1})_{(u,v)}|\,du\,dv&= \int_{R'}f(\phi^{-1}(u,v))\cdot \frac{1}{|\det D\phi_{\phi^{-1}(u,v)}|}\,du\,dv\\ &=\int_{R'}\left(f\cdot \frac{1}{|\det D\phi|}\right)(\phi^{-1}(u,v))\,du\,dv \end{align} Notice now that we are evaluating the product $f\cdot \frac{1}{|\det D\phi|}$ at the same point $\phi^{-1}(u,v)$. Surely you can see that for any $(x,y)\in\Bbb{R}^2$ we have $f(x,y)=x^2+y^2$ and \begin{align} \frac{1}{|\det D\phi_{(x,y)}|}=\frac{1}{\left|\det \begin{pmatrix} 2x&-2y\\2y & 2x\end{pmatrix}\right|} = \frac{1}{4(x^2+y^2)} \end{align} Therefore, the product $f\cdot \frac{1}{|\det D\phi|}=\frac{1}{4}$ is a constant function, meaning no matter at what point you plug in, the output is always $\frac{1}{4}$. Therefore, by continuing our calculation, \begin{align} \int_{R}(x^2+y^2)\,dx\,dy&=\int_{R'}\frac{1}{4}\,du\,dv=\frac{\text{area}(R')}{4}. \end{align}
I just want to highlight once again that the reason why we invoked the formula $|\det D(\phi^{-1})_{(u,v)}|=\frac{1}{|\det D\phi_{\phi^{-1}(u,v)}|}$ is because in formula $(*)$, the function $f$ is being evaluated at the point $\phi^{-1}(u,v)$ whereas $|det D(\phi^{-1})_{(u,v)}|$ is evaluated at $(u,v)$. This difference in the points of evaluation makes it difficult to explcicitly calculate stuff. So, by changing to the reciprocal, everything is now evaluated at the same point $\phi^{-1}(u,v)$, so the calculation is easier.
Let me emphasize: I did not have to do this reciprocal business. If I was a masochist, I could in principle literally go ahead an try to invert the relationship and find an explicit formula for $\phi^{-1}$, compose with $f$ and so on. This would have been unnecessarily tedious.