0

In Python 3.6 using NumPy I have two integer values x,y, two arrays rho_values, theta_values and an equation rho = x*sin(theta) + y*cos(theta).

I want to find all pairs of (rho, theta) with rho in rho_values and theta in theta_values which satisfy the equation.

According to this I came up with the following solution which I expected to give me a boolean matrix with true on every (rho,theta) index if rho, theta solve the equation:

diag = np.linalg.norm((256,400))  # Length of diagonal for some image shape (256,400)
theta_bins = np.linspace(-np.pi / 2, np.pi / 2, 300)
rho_bins = np.linspace(-diag, diag, 300)
def solve(x, y):
    return rho_bins == y*np.sin(theta_bins[:,np.newaxis])+x*np.cos(theta_bins[:,np.newaxis])

for x in range(0, 255):
    for y in range(0, 399):
        print(solve(x,y))

It does not work. What would be a correct and fast way to do this? Iterating over both arrays is not an option for me since I have many x,y values for each of which iterating over both arrays would cost a lot of time.

lasbr
  • 79
  • 1
  • 8
  • Have you occasionally written `row` when you meant `rho` (or vice versa)? Also, you once write `tetha` instead of `theta`. – jmd_dk May 03 '22 at 13:50
  • @jmd_dk Yes, sorry. I got them mixed up. I edited the question. – lasbr May 03 '22 at 13:52
  • Have you checked https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html? – Dilara Gokay May 03 '22 at 14:05
  • @Dilara Looks interesting, but I would prefer a plain Python and/or Numpy only solution – lasbr May 03 '22 at 14:07
  • @flawr Added the values. The example comes from a part of an implemenation of a Hough Transformation for images. This context should although not be necessary for my question – lasbr May 03 '22 at 14:29
  • Notice you have for x twice, should be for x for y. – jmmcd May 03 '22 at 14:35
  • 2
    Does this answer your question? [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) . Use `np.isclose`. – Jérôme Richard May 03 '22 at 14:44
  • @JérômeRichard Yes, thank you. This actually does. It is also similar to what jmmcd's answer says – lasbr May 03 '22 at 14:48

1 Answers1

0

As a general rule, using == with floating-point numbers is a mistake. Consider using an idiom such as abs(x-y)<0.00000001 instead.

jmmcd
  • 731
  • 6
  • 15
  • This actually was the problem. I was not familiar with the floating point math in Python – lasbr May 03 '22 at 14:48
  • 1
    Great! Notice that pretty much all modern languages use the same floating-point math (IEEE 754) so the same general rule applies in other languages also. – jmmcd May 04 '22 at 15:02