2016-08-03 3 views
4

Я пытается найти плоскости в точке облака 3d, используя формулу регрессии Z= aX + bY +CСамолет установки в 3D облако точек

Я реализовал наименьших квадратов и RANSAC решения, но уравнение 3 параметра ограничивает плоскость установки на 2.5D - формула не может применяться на плоскостях, параллельных оси Z.

Мой вопрос: Как я могу обобщить плоскость фитинга на полный 3d? Я хочу добавить четвертый параметр, чтобы получить полное уравнение a X + b Y + c * Z + d Как я могу избежать тривиального (0,0,0,0) решения?

Спасибо!

Кодекс Я использую:

from sklearn import linear_model 

def local_regression_plane_ransac(neighborhood): 
    """ 
    Computes parameters for a local regression plane using RANSAC 
    """ 

    XY = neighborhood[:,:2] 
    Z = neighborhood[:,2] 
    ransac = linear_model.RANSACRegressor(
              linear_model.LinearRegression(), 
              residual_threshold=0.1 
             ) 
    ransac.fit(XY, Z) 

    inlier_mask = ransac.inlier_mask_ 
    coeff = model_ransac.estimator_.coef_ 
    intercept = model_ransac.estimator_.intercept_ 
+0

Вы можете добавить ограничение, что (а, Ь, с) должны быть нормированными. – Julien

+0

Спасибо, не могли бы вы подробнее объяснить, как применить ограничение? – Tom

+0

Я никогда не использовал sklearn, но я уверен, что у него должен быть вариант где-то. В противном случае я знаю, что scipy.optimize имеет некоторые хорошие процедуры, которые это сделают. – Julien

ответ

7

Update

Эта функциональность теперь интегрирована в https://github.com/daavoo/pyntcloud и делает самолет фитинг процесс намного проще:

Учитывая точку помутнения:

enter image description here

Вам просто нужно добавить скалярное поле, как это:

is_floor = cloud.add_scalar_field("plane_fit") 

Wich будет добавить новый столбец со значением 1 для точек плоскости установлена.

Вы можете визуализировать скалярное поле:

enter image description here


Старый ответ

Я думаю, что вы могли бы легко использовать PCA, чтобы соответствовать плоскости к 3D точек вместо регрессии.

Вот простая реализация PCA:

def PCA(data, correlation = False, sort = True): 
""" Applies Principal Component Analysis to the data 

Parameters 
----------   
data: array 
    The array containing the data. The array must have NxM dimensions, where each 
    of the N rows represents a different individual record and each of the M columns 
    represents a different variable recorded for that individual record. 
     array([ 
     [V11, ... , V1m], 
     ..., 
     [Vn1, ... , Vnm]]) 

correlation(Optional) : bool 
     Set the type of matrix to be computed (see Notes): 
      If True compute the correlation matrix. 
      If False(Default) compute the covariance matrix. 

sort(Optional) : bool 
     Set the order that the eigenvalues/vectors will have 
      If True(Default) they will be sorted (from higher value to less). 
      If False they won't. 
Returns 
------- 
eigenvalues: (1,M) array 
    The eigenvalues of the corresponding matrix. 

eigenvector: (M,M) array 
    The eigenvectors of the corresponding matrix. 

Notes 
----- 
The correlation matrix is a better choice when there are different magnitudes 
representing the M variables. Use covariance matrix in other cases. 

""" 

mean = np.mean(data, axis=0) 

data_adjust = data - mean 

#: the data is transposed due to np.cov/corrcoef syntax 
if correlation: 

    matrix = np.corrcoef(data_adjust.T) 

else: 
    matrix = np.cov(data_adjust.T) 

eigenvalues, eigenvectors = np.linalg.eig(matrix) 

if sort: 
    #: sort eigenvalues and eigenvectors 
    sort = eigenvalues.argsort()[::-1] 
    eigenvalues = eigenvalues[sort] 
    eigenvectors = eigenvectors[:,sort] 

return eigenvalues, eigenvectors 

А вот как можно подправить точки на плоскости:

def best_fitting_plane(points, equation=False): 
""" Computes the best fitting plane of the given points 

Parameters 
----------   
points: array 
    The x,y,z coordinates corresponding to the points from which we want 
    to define the best fitting plane. Expected format: 
     array([ 
     [x1,y1,z1], 
     ..., 
     [xn,yn,zn]]) 

equation(Optional) : bool 
     Set the oputput plane format: 
      If True return the a,b,c,d coefficients of the plane. 
      If False(Default) return 1 Point and 1 Normal vector.  
Returns 
------- 
a, b, c, d : float 
    The coefficients solving the plane equation. 

or 

point, normal: array 
    The plane defined by 1 Point and 1 Normal vector. With format: 
    array([Px,Py,Pz]), array([Nx,Ny,Nz]) 

""" 

w, v = PCA(points) 

#: the normal of the plane is the last eigenvector 
normal = v[:,2] 

#: get a point from the plane 
point = np.mean(points, axis=0) 


if equation: 
    a, b, c = normal 
    d = -(np.dot(normal, point)) 
    return a, b, c, d 

else: 
    return point, normal  

Однако поскольку этот метод чувствителен к выбросам вы могли бы использовать RANSAC чтобы они соответствовали требованиям к выбросам.

Существует реализация Python ransac here.

И вам нужно только определить класс модели самолета, чтобы использовать его для установки плоскостей в 3D-точки.

В любом случае, если вы можете очистить 3D-точки от выбросов (возможно, вы можете использовать фильтр S.O.R KD-Tree), вы должны получить неплохие результаты с помощью PCA.

Вот реализация S.O.R:

def statistical_outilier_removal(kdtree, k=8, z_max=2): 
""" Compute a Statistical Outlier Removal filter on the given KDTree. 

Parameters 
----------       
kdtree: scipy's KDTree instance 
    The KDTree's structure which will be used to 
    compute the filter. 

k(Optional): int 
    The number of nearest neighbors wich will be used to estimate the 
    mean distance from each point to his nearest neighbors. 
    Default : 8 

z_max(Optional): int 
    The maximum Z score wich determines if the point is an outlier or 
    not. 

Returns 
------- 
sor_filter : boolean array 
    The boolean mask indicating wherever a point should be keeped or not. 
    The size of the boolean mask will be the same as the number of points 
    in the KDTree. 

Notes 
-----  
The 2 optional parameters (k and z_max) should be used in order to adjust 
the filter to the desired result. 

A HIGHER 'k' value will result(normally) in a HIGHER number of points trimmed. 

A LOWER 'z_max' value will result(normally) in a HIGHER number of points trimmed. 

""" 

distances, i = kdtree.query(kdtree.data, k=k, n_jobs=-1) 

z_distances = stats.zscore(np.mean(distances, axis=1)) 

sor_filter = abs(z_distances) < z_max 

return sor_filter 

Вы можете кормить функцию с KDtree ваших 3D-точек, вычисленных, возможно с помощью this implementation