Issue 1:
Hi,

I would like to find the best-fit axis of points that are on a cylindrical surface. Seems that scipy.linalg.svd is the function to look for.

So to test out, I decide to generate some points, function makeCylinder, from this thread https://stackoverflow.com/questions/22285994/how-to-generate-regular-points-on-cylindrical-surface, and estimate the axis.

This is the code:

def makeCylinder(radius, length, nlength, alpha, nalpha, center, orientation):

# Load
from numpy import array, allclose, linspace, tile, vstack
from numpy import pi, cos, sin, arccos, cross
from numpy.linalg import norm
from angoli import rotMatrixAxisAngle

# Create the length array
I = linspace(0, length, nlength)

# Create alpha array avoid duplication of endpoints
if int(alpha) == 360:
A = linspace(0, alpha, num=nalpha, endpoint=False)/180.0*pi
else:
A = linspace(0, alpha, num=nalpha)/180.0*pi

# Calculate X and Y
X = radius * cos(A)
Y = radius * sin(A)

# Tile/repeat indices so all unique pairs are present
pz = tile(I, nalpha)
px = X.repeat(nlength)
py = Y.repeat(nlength)

# Points
points = vstack(( pz, px, py )).T
## Shift to center
points += array(center) - points.mean(axis=0)

# Orient tube to new vector
ovec = orientation / norm(orientation)
cylvec = array([1,0,0])

if allclose(cylvec, ovec):
return points

# Get orthogonal axis and rotation
oaxis = cross(ovec, cylvec)
rot = arccos(ovec.dot(cylvec))

R = rotMatrixAxisAngle(oaxis, rot)

return points.dot(R)

from numpy.linalg import norm
from numpy.random import rand
from scipy.linalg import svd

for i in xrange(100):
orientation = rand(3)
orientation[0] = 0
orientation /= norm(orientation)

# Generate sample points
points = makeCylinder(radius = 3.0,
length = 20.0, nlength = 20,
alpha = 360, nalpha = 30,
center = [0,0,0],
orientation = orientation)
# Least Square
uu, dd, vv = svd(points - points.mean(axis=0))
asse = vv[0]

assert abs( abs(orientation.dot(asse)) - 1) <= 1e-4, orientation.dot(asse)

As you can see, I generate multiple cylinder whose axis is random (rand(3)).

The funny thing is that svd returns an axis that is absolutely perfect if the first component of orientation is zero (orientation[0] = 0).

If I comment this line the estimated axis is way off.

Even using leastsq on a cylinder equation returns the same behavior:

def bestLSQ1(points):
from numpy import array, sqrt
from scipy.optimize import leastsq

# Expand
points = array(points)
x = points[:,0]
y = points[:,1]
z = points[:,2]

# Calculate the distance of each points from the center (xc, yc, zc)
# http://geometry.puzzles.narkive.com/2HaVJ3XF/geometry-equation-of-an-arbitrary-orientated-cylinder
def calc_R(xc, yc, zc, u1, u2, u3):
return sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 - ( (x-xc)*u1 + (y-yc)*u2 + (z-zc)*u3 )**2 )

# Calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc, zc)
def dist(c):
Ri = calc_R(*c)
return Ri - Ri.mean()

# Axes - Minimize residu
xM, yM, zM = points.mean(axis=0)
# Calculate the center
center, ier = leastsq(dist, (xM, yM, zM, 0, 0, 1))
xc, yc, zc, u1, u2, u3 = center
asse = u1, u2, u3

return asse

This is the version that I am using:

import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
('0.19.0', '1.12.1', sys.version_info(major=2, minor=7, micro=13, releaselevel='final', serial=0))

Thanks

Issue 2:
My issue is about a possible memory leak in scipy.spatial.distance.pdist and scipy.spatial.distance.cdist

Reproducing code example:

import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
test = pdist(X/0.01, metric='sqeuclidean')

Error Description:

On my Mac OSX 10.14.6 the memory consumption of the snipped above quickly rises into the Gigabytes (after a few seconds). If I remove the division from the pdist argument (see below), the problem vanishes. (I don't think this is intended behavior, as it also leads to excessive memory consumption in the Gaussian Process Regression implementation of Sklearn, where they have a similar line in the implementation of the RBF kernel)

The code, which does not produce this problem (division removed):

import numpy as np
from scipy.spatial.distance import pdist
X = np.random.uniform(0, 100, size=10000).reshape(-1, 10)
for _ in range(10000):
test = pdist(X, metric='sqeuclidean')

Scipy/Numpy/Python version information:

1.7.0 1.21.0 sys.version_info(major=3, minor=8, micro=5, releaselevel='final', serial=0)