title: Fitting an Orthogonal Grid to Noisy Coordinates
body: <h3>Problem</h3>
  I have a list of coordinates that are meant to form a grid. Each coordinate has a random error component and some of the coordinates are missing. <i>Grid could be rotated (update).</i> I want to fit a orthogonal grid to the data points and return a list of the grid's vertices. For example:<br/><br/>
  <p><img src="https://i.stack.imgur.com/LatDf.jpg" alt="Grid Problem" /></p>
  <h3>Application</h3>
  The purpose is to find a grid in a scanned image. The data points come from the results of contour or edge detection in OpenCV. An example is image with a grid of photos.<br/><br/>
  <h3>Goal</h3>
  I wrote some Python code that works, but would like to find a linear algebra algorithm using SciPy, statsmodels or other modules that would be more robust and handle a small rotation of the grid (less than 10°).
  <br/><br/>
  <h3>Python Code Using Lists Only</h3>
  <pre><code># Noisy [x, y] coordinates (origin is upper-left corner)
  pts = [[103,101],
         [198,103],
         [300, 99],
         [ 97,205],
         [304,202],
         [102,295],
         [200,303],
         [104,405],
         [205,394],
         [298,401]]
  
  def row_col_avgs(num_list, ratio):
      # Finds the average of each row and column. Coordinates are
      # assigned to a row and column by specifying an error ratio.
      last_num, sum_nums, count_nums, avgs = 0, 0, 0, []
      num_list.sort()
      for num in num_list:
          # Calculate average for last row or column and begin new row or column
          if num &gt; (1+ratio)*last_num and count_nums != 0:
              avgs.append(int(round(sum_nums/count_nums,0)))
              sum_nums = num
              count_nums = 1
          # Or continue with current row or column
          else:
              sum_nums += num
              count_nums += 1
          last_num = num
      avgs.append(int(round(sum_nums/count_nums,0)))
      return avgs
  
  # Split coordinates into two lists of x's and y's
  xs, ys = map(list, zip(*pts))
  
  # Find averages of each row and column of the grid
  x_avgs = row_col_avgs(xs, 0.1)
  y_avgs = row_col_avgs(ys, 0.1)
  
  # Return vertices of completed averaged grid
  avg_grid = []
  for y_avg in y_avgs:
      avg_row = []
      for x_avg in x_avgs:
          avg_row.append([int(x_avg), int(y_avg)])
      avg_grid.append(avg_row)
  
  print(avg_grid)
  </code></pre>
  <h3>Output</h3>
  <pre><code>[[[102, 101], [201, 101], [301, 101]], 
   [[102, 204], [201, 204], [301, 204]], 
   [[102, 299], [201, 299], [301, 299]], 
   [[102, 400], [201, 400], [301, 400]]]
  </code></pre>