My Algorithm

The location of O was fixed at (0,0) but P and Q can be anywhere (but O!=P!=Q).
To simplify my algorithm I define that the right angle is always at O or P and never at Q.
The size of the grid is 50x50, each side has length size=50.

I identified 4 cases:
1. The right angle is at the origin O. P and Q can be anywhere on the x-/y-axis:
P(0<p_x<=size, 0) and Q(0, 0<q_y<=size) → that's size^2 different triangles.

2. The right angle is on the x-axis at P(0<p_x<=size, 0). Any Q with the same x-value makes a valid triangle:
P(0<p_x<=size, 0) and Q(p_x, 0<q_y<=size) → that's size^2 different triangles.

3. The right angle is on the y-axis at P(0, 0<p_y<=size). Any Q with the same y-value makes a valid triangle:
P(0, 0<p_y<=size) and Q(0<q_x<=size, p_y) → that's size^2 different triangles.

4. The right angle is "inside" the grid at P(0<p_x<=size, 0<p_y<=size).

The fourth case is a more challenging case and will be discussed in detail:
The right angle is enclosed by two vectors \vec{a} = \vec{PO} and \vec{b} = \vec{PQ}.
If they are perpendicular \vec{a} \bot \vec{b}, then their dot product is zero: \vec{a} \cdotp \vec{b} = a_{x}b_x + a_{y}b_y = 0

The task is to find for a given point P(p_x, p_y) each Q(q_x, q_y) such that p_x(q_x - p_x) + p_y(q_y - p_y) = 0.
For example P(2,2) and Q(3,1) create a valid right angle:
p_x(q_x - p_x) + p_y(q_y - p_y) = 2(3-2) + 2(1-2) = 2 - 2 = 0

The dot product can be transformed to
p_x(q_x - p_x) = -p_y(q_y - p_y)
This formula is obivously true for q_x = p_x and q_y = p_y.
Even more, it's also true for q_x = p_x + k \cdot p_y and q_y = p_y - k \cdot p_x (last line was just the special case k=0).
p_x(p_x + k p_y - p_x) = -p_y(p_y - k p_x - p_y)
p_x(k p_y) = -p_y(-k p_x)
kp_{x}p_y = kp_{x}p_y

Two nested for-loops produce all integer combinations of 0 < p_x <= size and 0 < p_y <= size.
We only want integer solutions for q_x and q_y, too.
But if k is an integer, too, then we would miss some valid Q - e.g. for P(2,2) and Q(3,1) it follows that k=frac{1}{2}.
We get all integer solutions if k is a multiple of dfrac{1}{gcd(p_x,p_y)} (the Greatest Common Divisor, see en.wikipedia.org/wiki/Greatest_common_divisor)
For P(2,2) we get dfrac{1}{gcd(2,2)}=dfrac{1}{2}

Enough mathematics, let's dive into the code ... cases 1 to 3 don't need much of an explanation.
Then follows two nested for-loop that analyze each P(p_x,p_y) of case 4.
factor is the gcd(p_x,p_y) such that k_{n+1} = k_n + 1/factor
The point Q is found by repeated adding p_y/factor to p_x and subtracting p_x/factor from p_y (for k > 0)
until any coordinate leaves the grid. The same procedure is done with inverse sign for k < 0.

A simple speed-up was achieved by observing that the number of valid Q is the same for P(p_x, p_y) and P(p_y, p_x).
That means that the two nested loops only need to considered the bottom-right triangle of the grid: P(p_x, p_y<=p_x).
For P(p_x, p_y<p_x) I multiply the number of Qs by 2 - but not when p_x=p_y.

Alternative Approaches

Most other solutions replace my inner-most while-loops by a closed formula:
min(y * factor / x, (size - x) * factor / y) * 2
Noone really explains that formula and surprisingly it isn't faster on my machine for size <= 2500.

The original problem can be solved by brute-force, too (four nested loops for p_x, p_y, q_x, q_y).

Modifications by HackerRank

Grids can have an edge length of up to 2500.
The brute-force approach fails of course (time-out) while my code finishes in less than 0.22 seconds.