Calculate average spacing per direction

Question

I am interested in finding the average spacing between objects on a directional basis. I did my best to minimally implement this in Python and to conform to PEP8.

I try to explain what I'm doing in the code comments. The problem I want to solve is to find the average spacing between circles of known radius and locations (in my application, these locations are irregular) on a directional basis. For example, what is the average spacing between these circles in the x, diagonal, and y directions?

import numpy as np

# I have a set of x and y coordinates describing the
# locations of objects
px = np.array([0, 0, 0, 10, 10, 10, 20, 20, 20])
py = np.array([0, 10, 20, 0, 10, 20, 0, 10, 20])

# Each object is a circle with known radius size
point_diameter = 2

# I'll put my answers here
spacing = {}

# The challenge is to calculate the average spacing per point on a directional
#  basis, e.g., what is the average distance between points in the direction
#  of 45 degrees?


# To do this, I define a function to check if points are alligned in a
#  direction
def check_alligned(direction, p1x, p1y, p2x, p2y):

    # Since the points have some wisdth, I need to check the "edge cases"
    # as well as the centerlines.
    tolerance = 0.2
    pr = point_diameter / 2
    for offsetx1 in [-1, 0, 1]:
        for offsetx2 in [-1, 0, 1]:
            for offsety1 in [-1, 0, 1]:
                for offsety2 in [-1, 0, 1]:

                    # The scheme is to check if the direction leads to the object x or y location. 
                    # I chose or to handle cases where the x or y difference is zero.
                    distance = distance_between([p1x + pr*offsetx1,
                                                     p1y + pr*offsety1],
                                                    [p2x + pr*offsetx2,
                                                     p2y + pr*offsety2])
                    if abs(p1x + distance * np.cos(np.radians(direction)) - p2x) <= tolerance or abs(p1y + distance * np.sin(np.radians(direction)) - p2y) <= tolerance:
                         return True
    return False


def distance_between(v1, v2):
    return np.sqrt( (v1[0] - v2[0])**2 + (v1[1] - v2[1])**2)

# Initial Smallest Spacing - use to find if we get to a spacing in the loop
# and record the smallest one
ISP = 9e4

# Iterate through directions to calculate average distance in each direction.
# Directions go counterclockwise and x points in the direction of 0 degrees.
for direction in [0, 60, 45, 90]:

    # These are a running sum and counter of 1st nearest neighbor spacings
    # that are less than an initial smallest spacing value. They are used
    # to calculate the average spacing in each direction.
    total_spacing = 0.
    space_count = 0

    # Look through every point
    for ii in range(px.size):

        smallest_spacing = ISP
        # and every other point except ones we've already checked
        for jj in range(px.size - ii - 1):

            # If the jj point is not in the same line as the
            # ii point along this diretion, skip it
            idx2 = ii + jj + 1

            if check_alligned(direction, px[ii], py[ii], px[idx2], py[idx2]):
                space = np.sqrt((px[ii]-px[idx2])**2 +
                                          (py[ii] - py[idx2])**2)
                if space < smallest_spacing:
                    smallest_spacing = space

        # Record a distance if it is less than the initial
        # smallest spacing
        if smallest_spacing != ISP:
            total_spacing += smallest_spacing
            space_count += 1

    # record the average spacing in this direction
    if space_count:
        spacing[direction] = total_spacing / space_count
    else:
        spacing[direction] = 0.

print spacing

Answer 1

I don't fully understand your problem statement, could you elaborate? I asked about your definition of distance because unless the distance you are interested in is some projection of the euclidian distance, it should be constant regardless which direction you are looking.

In any case, I did your problem manually and found a different result than your program above reports. So here's my implementation with inline comments:

#!/usr/local/bin/python3

from sklearn.neighbors import NearestNeighbors
import numpy as np
import matplotlib.pyplot as plt

class Turbine():

    def __init__(self, x, y):
        self.initialx = x
        self.initialy = y
        self.x = x
        self.y = y

    def rotateBy(self, theta):
        newX = self.x*np.cos(theta)-self.y*np.sin(theta)
        newY = self.x*np.sin(theta)+self.y*np.cos(theta)
        self.x = newX
        self.y = newY

    def point(self):
        return [self.x, self.y]

class Farm():
    def __init__(self, turbineList):
        self.turbineList = turbineList

    def rotateFarmBy(self, theta):
        for turbine in self.turbineList:
            turbine.rotateBy(theta)

    def turbineCoordinates(self):
        return np.array([turbine.point() for turbine in self.turbineList])

    def turbineCount(self):
        return len(self.turbineList)


for theta in [0, np.pi/4., np.pi/3., np.pi/2.]:

    # establish the coordinates,
    # but first test with a small subset and compare with hand calculations
    coordinates = [[0, 0], [0, 10], [0, 20]]
    # coordinates = [[0, 0], [10, 10], [20, 20]]
    # coordinates = [[0, 0], [0, 10], [0, 20], [10, 0], [10, 10], [10, 20], [20, 0], [20, 10], [20, 20]]

    # instantiate the turbines and the container farm
    turbines = [Turbine(coord[0], coord[1]) for coord in coordinates]
    farm = Farm(turbines)
    farm.rotateFarmBy(theta)

    # get the N nearest neighbors, where N in this case is the number of turbines in the farm
    coordinates = farm.turbineCoordinates()
    neighbors = NearestNeighbors(n_neighbors=farm.turbineCount(), algorithm='ball_tree').fit(coordinates)
    distances, indices = neighbors.kneighbors(coordinates)

    # remove any zero distances so that you dont count each turbine twice
    nonzeros = distances[distances > 0]

    print("at {0:0.1f} degrees, the average distance is {1:0.3f}".format(theta*180/np.pi, sum(nonzeros)/farm.turbineCount()))

    # print(indices)
    # print(distances)
    # plt.scatter([c[0] for c in coordinates], [c[1] for c in coordinates])
    # plt.show()

Answer 2

Your "irregular" location setup is intriguing, if a bit vague.

alternatives

There are several approaches to identifying nearest neighbor, including exhaustively enumerating candidates as in your code, or sorting by X and/or Y to focus on fewer candidates, or sorting by bit-interleaved X+Y to give one kind of quadtree representation. One can also quantize to support hashing.

I'm going to forego such fancy techniques in favor of a simple, visual argument. It relies on the notion of "in expectation" when finding the average.

proposal

Based on your input data and desired accuracy, pick a pixel grid size, and render your circular objects on the grid.

Cycle through your desired angle settings. Uniformly pick a random location on the grid, or to the left on the Y-axis. Constrain to a conveniently small bounding box in the grid, or bounding interval on the Y-axis, so we are sure to run into objects.

From the uniformly distributed starting point, start drawing a line (a ray) in the desired direction. Note when you enter a rendered object, and when you exit an object. Once an "exit" event happens continue extending the ray until next "enter" event, and record the freespace distance.

At this point, your specification doesn't make it clear to me if you want to keep extending the ray to make additional freespace measurements since the irregular locations leave varying spaces, or if you prefer to end the iteration at this point.

Now iterate: go back to choosing another random starting point, and make additional distance measurements, until you have N of them.

Finally, report the average as sum(distances) / N.