In a comment to the question $n^2+(n+1)^2+n^2\cdot(n+1)^2$ is a perfect square it is proved that:
$(1)\quad b=a+1$ give integer solutions to
$a^2+b^2+(ab)^2=c^2$ for all $a\in\mathbb N$.
In the answers to the question $\forall m\in\mathbb N\exists n>m+1\exists N\in\mathbb N:m^2+n^2+(mn)^2=N^2$ it is proved that:
$(2)\quad b=2a^2$ give integer soulutions for all $a\in\mathbb N$.
Conjecture:
For each $a\in\mathbb N$ there are integer solutions to $a^2+b^2+(ab)^2=c^2$ that are neither of the types $(1)$ or $(2)$ above.
Prove the conjecture or give a counter-example. It is tested for all $a<1000$.
The Pell type equation $$ x^2 - (1 + A^2) y^2 = A^2 $$ has an infinite set of solutions $(x,y).$ The first three predictable types are $$ \left( \begin{array}{c} A \\ 0 \end{array} \right), $$ $$ \left( \begin{array}{c} A^2 - A + 1 \\ A - 1 \end{array} \right), $$ $$ \left( \begin{array}{c} A^2 + A + 1 \\ A + 1 \end{array} \right). $$ After that, an infinite sequence of $(x,y)$ as column vectors may be found by multiplying by the generator of the (oriented) automorphism group of the binary quadratic form $ x^2 - (1 + A^2) y^2. $ That matrix is $$ M = \left( \begin{array}{cc} 2A^2 + 1 & 2 A^3 + 2 A \\ 2A & 2 A^2 + 1 \end{array} \right). $$
The answer you write as $b = 2 a^2$ comes up as $$ M = \left( \begin{array}{cc} 2A^2 + 1 & 2 A^3 + 2 A \\ 2A & 2 A^2 + 1 \end{array} \right) \left( \begin{array}{c} A \\ 0 \end{array} \right) = \left( \begin{array}{c} 2A^3 + A \\ 2 A^2 \end{array} \right). $$
The bad news is that, for a fixed $A,$ there may be others. These others show up as later entries in the Pell sequences for smaller values of $A.$ I do not have a genuinely sensible two-dimensional description for all answers.