Hough transform: How to get lines from voting-matrix?

Question

so Im trying to implement the hough transform using python and c++ (using Pybind11 for interfacing between the two languages). When Im plotting the hough space it seems alright but I just can't get a line from the maximum of the voting matrix. Here is the C++ code (looks a bit different bc I use PyBind11):

py::array_t<int> houghTransform(py::array_t<int> image, int angleStep, int angleAmount) {
    auto imageBuf = image.mutable_unchecked<3>();
    int height = imageBuf.shape(0);
    int width = imageBuf.shape(1);

    py::array_t<int> edgeMatrix = edgeDetect(imageBuf, height, width);
    auto edgeMatrixBuf = edgeMatrix.mutable_unchecked<2>();

    int distanceAxis = 2 * sqrt(pow((float) height, 2.0) + pow((float) width, 2.0));
    int angleAxis = angleAmount;
    
    int angleDim = (int) angleAxis / angleStep;
    int distanceDim = (int) distanceAxis / 2;

    py::array_t<int> votingMatrix = py::array_t<int>({distanceAxis, angleDim});
    auto votingMatrixBuf = votingMatrix.mutable_unchecked<2>();

    // fill voting matrices with zeros
    for(int i=0; i<distanceDim; i++) {
        for(int j=0; j<angleDim; j++) {
            votingMatrixBuf(i, j) = 0;
        }
    }

    // vote
    for(int x=0; x<edgeMatrixBuf.shape(0); x++) {
        for(int y=0; y<edgeMatrixBuf.shape(1); y++) {
            if(edgeMatrixBuf(x, y) == 1) {

                int counter = 0;
                float theta;
                float ro;

                for(int thetaIdx=0; thetaIdx<=angleAxis; thetaIdx++) {
                    if(thetaIdx % angleStep == 0) {
                        counter++;
                        theta = (float) (thetaIdx) * (M_PI / 180);
                        ro = distanceDim + std::round((x * cos(theta)) + (y * sin(theta)));
                        votingMatrixBuf(ro, counter) += 1;
                    }
                }

            }
        }
    }

    return votingMatrix;
}

As you can see the arguments of the functions are the image matrix, which I transform to a matrix where the edges are 1 and the rest 0, so I got my pixels of interest. int angleAmount is what angle range I want to try outand int angleStep is how many of angles of that theta I really want to use (for example, skip every second theta). But for this example I will use angleAmount = 360 and angleStep = 1. So I will use all angles form 1 to 360.

Here is the python code:

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import time

from houghTransform import houghTransform


def apply_hough_transform(image_path: str=""):
    image = np.array(Image.open(image_path))
    lines = houghTransform(image, 1, 360)

    p = np.unravel_index(lines.argmax(), lines.shape)
    
    max_distance = 2 * np.sqrt(pow(image.shape[0], 2) + pow(image.shape[1], 2))

    ro = p[0] - (max_distance / 2)
    theta = p[1] * (np.pi / 180)

    a = np.cos(theta)
    b = np.sin(theta)
    x = a * ro
    y = b * ro

    pt1 = (int(x + 1000*(-b)), int(y + 1000*(a)))
    pt2 = (int(x - 1000*(-b)), int(y - 1000*(a)))
    
    fig, axs = plt.subplots(2)

    axs[0].matshow(lines)
    axs[0].scatter(p[1], p[0], facecolors="none", edgecolors="r")

    axs[1].plot([pt1[0], pt2[0]], [pt1[1], pt2[1]])
    axs[1].imshow(image)
    plt.show()

apply_hough_transform(image_path="images/black_line.png")

The function houghTransform is the same as in the c++ code which I exported to Python using PyBind11. Here are the images:

I also tried to create the line using this function:

def line(x):
    return -(1 / np.arctan(theta)) * (x - ro * np.cos(theta)) + ro * np.sin(theta)

But it also didnt work. Can you spot my error? Im sitting on this for quite some time so help is really appreciated!