How can we convert a probability density according to solid angle to a density according to area? In Chapter 14.5.3 of Physically Based Rendering - From Theory to Implementation (3rd edition), the authors claim that a probability density $p_\omega$ according to solid angle can be converted to a density $p_A$ according to area by $$p_A=p_\omega\frac{\left|\cos\theta_i\right|}{\left|p_i-p_{i+1}\right|^2},\tag1$$ where $\theta_i$ is the angle between the surface normal at $p_{i+1}$ and the vector to $p_i$.


  The problem is: $(1)$ doesn't make sense. If $p_\omega$ is a density with respect to the surface measure on the upper unit hemisphere $\Omega_{p_i}$ over $p_i$, then $p_\omega$ is by definition a function $\Omega_{p_i}\to[0,\infty]$. Analogously, if $p_A$ is supposed to be a density with respect to the surface measure on the manifold $M$ build by all surfaces in the scene, then $p_A$ is by definition a function $M\to[0,\infty]$.
  
  So, there's at least an argument transformation missing. I would really appreciate, if anybody could clarify.
I don't understand what's the problem. We want to make a change of variables so that instead of integrating over the range of directions $d\omega_i$ in the upper hemisphere we integrate over patches of area $dA$ of a region $R$.

This whole thing is given in section 5.5 as Hubble said but let me post it again.

We know by definition that,

$\omega = A *  cos(\theta)/r^2$  

Now if you differentiate this, you get

$d\omega = dA * cos(\theta) / r^2$ 

where,

$Jacobian = cos(\theta) / r^2$

as Hubble said. I don't know what are you trying to do by calculating Jacobian through determinant but as far as I know, that Jacobian determinant is for multi-variables and involves partial derivatives. For change of variables involving a single variable the Jacobian can be gotten by simply derivating as I did.

Now this might seem the inverse but wait a moment. What you posted in the question is a relation between the PDFs.

Consider a single sample of the rendering equation integrated over the domain of $\omega_i$

$L_s = BRDF * L_f * (\omega_i . n) / p(\omega)  $

After changing the domain by above equations this would become

$L_s = BRDF * L_f * (\omega_i . n) (\omega_i . n_r)\, / \,  (p(A) * r^2)  $

$L_s = BRDF * L_f * (\omega_i . n) * JN\, / \,  p(A) $
where $Jacobian = JN$

This means, 
$ 1 / p(\omega) =  JN/p(A)$

$ p(A) = p(\omega) * JN$

$ p(A) = p(\omega) * cos(\theta)/r^2$

This is what you see above. Simple as ****.