Astroua/TurbuStat

SCF distance metric is incorrectly weighted.

Closed this issue · 1 comments

In the distance_metric method of the SCF_Distance class, the difference between SCF surfaces is weighted by the inverse distance to center of the SCF surface. This is calculated here as dist_weight, where a and b are the x and y coordinates of each pixel in the SCF surface.

if self.weighted:
# Centre pixel set to 1
a[np.where(a == 0)] = 1.
b[np.where(b == 0)] = 1.
dist_weight = 1 / np.sqrt(a ** 2 + b ** 2)

I believe the central pixel is weighted incorrectly. The central pixel should have a weight of 1. But setting a and b of the central pixel to 1 means that dist_weight of the central pixel is dist_weight = 1 / sqrt(1^2 + 1^2) = 1 / sqrt(2). This calculation is not only incorrect for the central pixel, but for all pixels in the SCF surface with a or b equal to 0. Take for example, the pixel just 'above' the center, with a = 0 and b = 1. This pixel has a distance of 1 from the center. However, since the above code replaces a = 0 with a = 1, the weight for this pixel is also 1 / sqrt(2).

The following is a suggestion that correctly calculates the weight for each pixel as the inverse of its distance from center, except in the case of the central pixel where this would give infinite weight. The central pixel weight is replaced to 1.

if self.weighted:
            # Centre pixel set to 1
            a[np.where((a == 0) & (b == 0))] = 1.
            dist_weight = 1 / np.sqrt(a ** 2 + b ** 2)

Another more intuitive method which gives the same result, but gives a divide by zero warning:

if self.weighted:
            dist_weight = 1 / np.sqrt(a ** 2 + b ** 2)
            # Centre pixel set to 1
            dist_weight[np.where((a == 0) & (b == 0))] = 1.

The difference between the array of weights calculated by the method in turbustat and my method is shown below:
scf_weight_turbustat
scf_weight_correct
scf_weight_compare

Thanks for finding this @jrobbfed! Would you be able to issue a PR with the fix? Your second solution looks like the way to go.