This is a neat observation. I was unable to find historical references for it in particular, but there is a rich history of solving this type of problems. The problem is clearly reminiscent of finding the Pythagorean triples, right triangles with integer sides, for which Elements X.29 gives a solution without a hint of how it was obtained (of course once the rule is known it is easy to check that it works). A method first appears in Diophantus's Arithmetica (c. 250 AD) when solving the famous problem II.8, "partition a given square into two squares", the one whose margins were too narrow in Fermat's copy of the book. Schappacher's article (pp.12-13) gives a nice commentary.
In modern retelling if $p,q,r$ are the sides then for $x:=p/r$, $y:=q/r$ we have $x^2+y^2=1$ by the Pythagorean theorem. Write $y^2=1-x^2$ and make substitution $y=t(1-x)$ then equation for $x$ is linear and gives $x=(t^2-1)/(t^2+1)$, $y=2t/(t^2+1)$. If $t$ is rational so are $x$ and $y$, and setting $t=n/m$ gives the Pythagorean triples $n^2-m^2$, $2nm$, $n^2-m^2$. The trick can be interpreted geometrically as intersecting the circle with lines of slope $-t$ passing through a point $(1,0)$, and it works because each line has exactly one other intersection point with the circle. Diophantus describes nothing of the sort of course. Schappacher's remarks (p.27) suggest that Poincare might have come up with this interpretation in 1901, it is now called rational parametrization in algebraic geometry.
To extend this to self-trisecting triangles we first need a "Pythagorean theorem". Just like having a right angle makes triangle's angles into $\alpha,\pi/2-\alpha,\pi/2$ the trisection condition makes them into $\alpha,3\alpha,\pi-4\alpha$. In the first case writing equations of the theorem of sines and eliminating $\alpha$ from them leads to the Pythagorean equation $p^2+q^2=r^2$, in the second case after some fun with more esoteric trig identities one gets $q^2r=(p+r)(p-r)^2$. This is the "Pythagorean theorem" for self-trisecting triangles. With the same notation as above the corresponding curve is a cubic $$y^2=(x+1)(x-1)^2.$$ Normally cubics do not have rational parametrizations, but this one is special because of the squared factor, it means that it has a node (self-intersection point) at $(1,0)$. Drawing lines through the node, i.e. making the substitution $y=t(x-1)$, immediately gives $x=t^2-1$ and $y=t(t^2-2)$. Setting $t=n/m$ gives us $$p=m(n^2-m^2),\ q=n(n^2-2m^2),\ r=m^3,$$ as in the OP up to notation. So not only does the formula work, it actually lists all self-trisecting triples. The nodal cubic looks something like the folium of Descartes (rotated by $135^\circ$ clockwise, scaled, and shifted by $1$ to the right so that the $x$-intercepts move to $(-1,0)$ and $(1,0)$), but it does not have an asymptote.