Find intersection in noisy Hough transform

Question

I have various sets of points, which are supposed to be part of different lines (see image below). What I want to do is to group the points which are part of the same line. To do so, I made a Hough transform that turns cartesian coordinates into polar one. What I need to do is to locate the position of the intersection in the Hough plane; do you have any suggestion on how to do it?

This is what I am doing now:

1. Apply a threshold to the image of the Hough transform, so to consider only the pixels with at least half of the maximum intensity

   max_Hough = max(max(Hough));
   for ii = 1:dim_X
       for jj = 1:dim_Y
           if Hough(ii,jj) > max_Hough*0.5
               Hough_bin(ii,jj) = Hough(ii,jj);
           end
       end
   end
   

2. Locate the blobs in the corresponding image

   se = strel('disk',1);
   Hough_final = imclose(Hough_bin,se);
   [L,num] = bwlabel(Hough_final);   
   % Let's exclude the single points
   counts = sum(bsxfun(@eq,L(:),1:num));
   valid_blobs = find(counts>threshold);
   number_valid = length(find(counts>threshold));
   
   if (number_valid~=0)
      for blob_num = (1:number_valid)
          for x = (1:dim_X)
              for y = (1:dim_X)
                  if (L(x,y) == valid_blobs(blob_num))     
                      hist_bin_clean(x,y) = Hough_final(x,y);
                  end
              end
          end
      end
   end
   

3. Find the centroids of the blobs

   Ilabel = bwlabel(hist_bin_clean);
   stat = regionprops(Ilabel,'centroid');
   

You can see the result in the second image below, with the centroid pointed by the white arrow. I would expect it to be a few pixels higher and to the left... How would you improve the result?