For a 1D conservative Hamiltonian like $H(q,p)$, we can transform to Angle-Action (AA) variables using the following formulae, which are based on area-conservation of this canonical transformation (Ref)
$$ I(E)=\frac{1}{\pi}\int\limits^{q_2}_{q_1}dq p(q,E), $$ $$ \theta(I)=\frac{\partial}{\partial I}\int\limits^{q}_{0}dq p(q,I), $$ $$ \theta(t)=\frac{\partial H(I)}{\partial I}t+\theta_{0}=\omega t+\theta_{0}, $$
where $E$ is any specific energy level (=$H$), $(q_{1},q_{2})$ are the lowest and highest points that are allowed by a potential function $V(q)$ for the given energy $E$ considered (covering the closed phase-space being integrated, for an oscillation or rotation).
My question now is how do we extend this to 2D or 3D cases (i.e. 4D and 6D phase spaces), where we have $H(q_{1},q_{2},p_{1},p_{2})$ or $H(q_{1},q_{2},q_{3},p_{1},p_{2},p_{3})$, respectively? Does every dimension pair, say $(q_{k},p_{k})$, map to a new angle-action pair, say $(\theta_{k},I_{k})$, and do we use the same formulae above for each pair in isolation of the others since the variables $q_{i}$ or $p_{i}$ are independent? Or do we try to somehow find them by making perhaps the integrals above multidimensional in 2D or 3D, to integrate the phase space in 2D or 3D now? Or is it some other concept?
Is there a clear and concise explanation for this extension and its calculation procedure, perhaps with some solved examples?