python optimise.leastsq: подгонка круга к 3D-набору точек

Я пытаюсь использовать код подгонки окружности для набора 3D-данных. Я изменил его для 3D-точек, просто добавив z-координату, где это необходимо. Моя модификация отлично работает для одного набора точек и плохо для другого. Пожалуйста, посмотрите код, если в нем есть ошибки.

import trig_items
import numpy as np
from trig_items import *
from numpy import *
from matplotlib import pyplot as p
from scipy import optimize

# Coordinates of the 3D points
##x = r_[36, 36, 19, 18, 33, 26]
##y = r_[14, 10, 28, 31, 18, 26]
##z = r_[0, 1, 2, 3, 4, 5]

x = r_[ 2144.18908574,  2144.26880854,  2144.05552972,  2143.90303742,  2143.62520676,
  2143.43628579,  2143.14005775,  2142.79919654,  2142.51436023,  2142.11240866,
  2141.68564346,  2141.29333828,  2140.92596405,  2140.3475612,   2139.90848046,
  2139.24661021,  2138.67384709,  2138.03313547,  2137.40301734,  2137.40908256,
  2137.06611224,  2136.50943781,  2136.0553113,   2135.50313189,  2135.07049922,
  2134.62098139,  2134.10459535,  2133.50838433,  2130.6600465,   2130.03537342,
  2130.04047644,  2128.83522468,  2127.79827542,  2126.43513385,  2125.36700593,
  2124.00350543,  2122.68564431,  2121.20709478,  2119.79047011,  2118.38417647,
  2116.90063343,  2115.52685778,  2113.82246629,  2112.21159431,  2110.63180117,
  2109.00713198,  2108.94434529,  2106.82777156,  2100.62343757,  2098.5090226,
  2096.28787738,  2093.91550703,  2091.66075061,  2089.15316429,  2086.69753869,
  2084.3002414,   2081.87590579,  2079.19141866,  2076.5394574,   2073.89128676,
  2071.18786213]
y = r_[ 725.74913818,  724.43874065,  723.15226506,  720.45950581,  717.77827954,
  715.07048092,  712.39633862,  709.73267688,  707.06039438,  704.43405908,
  701.80074596,  699.15371526,  696.5309022,   693.96109921,  691.35585501,
  688.83496327,  686.32148661,  683.80286662,  681.30705568,  681.30530975,
  679.66483676,  678.01922321,  676.32721779,  674.6667554,   672.9658024,
  671.23686095,  669.52021535,  667.84999077,  659.19757984,  657.46179949,
  657.45700508,  654.46901086,  651.38177517,  648.41739432,  645.32356976,
  642.39034578,  639.42628453,  636.51107198,  633.57732055,  630.63825133,
  627.75308356,  624.80162215,  622.01980232,  619.18814892,  616.37688894,
  613.57400131,  613.61535723,  610.4724493,   600.98277781,  597.84782844,
  594.75983001,  591.77946964,  588.74874068,  585.84525834,  582.92311166,
  579.99564481,  577.06666417,  574.30782762,  571.54115037,  568.79760614,
  566.08551098]
z = r_[ 339.77146775,  339.60021095,  339.47645894,  339.47130963,  339.37216218,
  339.4126132,   339.67942046,  339.40917728,  339.39500353,  339.15041461,
  339.38959195,  339.3358209,   339.47764895,  339.17854867,  339.14624071,
  339.16403926,  339.02308811,  339.27011082,  338.97684183,  338.95087698,
  338.97321177,  339.02175448,  339.02543922,  338.88725411,  339.06942374,
  339.0557553,   339.04414618,  338.89234303,  338.95572249,  339.00880416,
  339.00413073,  338.91080374,  338.98214758,  339.01135789,  338.96393537,
  338.73446188,  338.62784913,  338.72443217,  338.74880562,  338.69090173,
  338.50765186,  338.49056867,  338.57353355,  338.6196255,   338.43754399,
  338.27218569,  338.10587265,  338.43880881,  338.28962141,  338.14338705,
  338.25784154,  338.49792568,  338.15572139,  338.52967693,  338.4594245,
  338.1511823,   338.03711207,  338.19144663,  338.22022045,  338.29032321,
  337.8623197 ]

# coordinates of the barycenter
xm = mean(x)
ym = mean(y)
zm = mean(z)

### Basic usage of optimize.leastsq

def calc_R(xc, yc, zc):
    """ calculate the distance of each 3D points from the center (xc, yc, zc) """
    return sqrt((x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2)

def func(c):
    """ calculate the algebraic distance between the 3D points and the mean circle centered at c=(xc, yc, zc) """
    Ri = calc_R(*c)
    return Ri - Ri.mean()

center_estimate = xm, ym, zm
center, ier = optimize.leastsq(func, center_estimate)
##print center

xc, yc, zc = center
Ri       = calc_R(xc, yc, zc)
R        = Ri.mean()
residu   = sum((Ri - R)**2)
print 'R =', R

Итак, для первого набора x, y, z (прокомментированного в коде) это работает хорошо: выход R = 39.0097846735. Если я запускаю код со вторым набором точек (без комментариев), результирующий радиус равен R = 108576.859834, что является почти прямой линией. Я нарисовал последний.

Синие точки — заданный набор данных, красные — дуга результирующего радиуса R = 108576.859834. Очевидно, что данный набор данных имеет гораздо меньший радиус, чем результат.

Вот еще набор точек.

Понятно, что метод наименьших квадратов работает некорректно.

Пожалуйста, помогите мне решить эту проблему.

ОБНОВЛЕНИЕ

Вот мое решение:

### fit 3D arc into a set of 3D points             ###
### output is the centre and the radius of the arc ###
def fitArc3d(arr, eps = 0.0001):
    # Coordinates of the 3D points
    x = numpy.array([arr[k][0] for k in range(len(arr))])
    y = numpy.array([arr[k][4] for k in range(len(arr))])
    z = numpy.array([arr[k][5] for k in range(len(arr))])
    # coordinates of the barycenter
    xm = mean(x)
    ym = mean(y)
    zm = mean(z)
    ### gradient descent minimisation method ###
    pnts = [[x[k], y[k], z[k]] for k in range(len(x))]
    meanP = Point(xm, ym, zm) # mean point
    Ri = [Point(*meanP).distance(Point(*pnts[k])) for k in range(len(pnts))] # radii to the points
    Rm = math.fsum(Ri) / len(Ri) # mean radius
    dR = Rm + 10 # difference between mean radii
    alpha = 0.1
    c = meanP
    cArr = []
    while dR  > eps:
        cArr.append(c)
        Jx = math.fsum([2 * (x[k] - c[0]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        Jy = math.fsum([2 * (y[k] - c[1]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        Jz = math.fsum([2 * (z[k] - c[2]) * (Ri[k] - Rm) / Ri[k] for k in range(len(Ri))])
        gradJ = [Jx, Jy, Jz] # find gradient
        c = [c[k] + alpha * gradJ[k] for k in range(len(c)) if len(c) == len(gradJ)] # find new centre point
        Ri = [Point(*c).distance(Point(*pnts[k])) for k in range(len(pnts))] # calculate new radii
        RmOld = Rm
        Rm = math.fsum(Ri) / len(Ri) # calculate new mean radius
        dR = abs(Rm - RmOld) # new difference between mean radii

    return Point(*c), Rm

Это не очень оптимальный код (у меня нет времени его дорабатывать), но он работает.


python optimization least-squares curve-fitting
person Community    schedule 18.03.2013    source источник
comment
Вы можете найти полезные решения в 2D-эквивалентном вопросе круг"> stackoverflow.com/questions/14834693/   -  person Hooked    schedule 03.04.2013
comment
Если ваши данные обычно выглядят как параболические синие точки здесь, подгонка плоскости должна быть очень простой. Следовательно, вы можете легко ограничить круговую подгонку 2D-версией. (см. мое обновление ниже. Хотя вы можете использовать разные данные для плоскости и для окружности)   -  person mikuszefski    schedule 17.04.2013

Ответы (1)


arrow_upward
6
arrow_downward

Я предполагаю, что проблема заключается в данных и соответствующем алгоритме. Метод наименьших квадратов отлично работает, если он дает локальный параболический минимум, так что метод простого градиента проходит примерно по направлению минимума. К сожалению, это не обязательно относится к вашим данным. Вы можете проверить это, сохранив некоторые грубые оценки для xc и yc фиксированными и построив сумму квадратов остатков как функцию zc и R. У меня получается минимум в форме бумеранга. В зависимости от ваших начальных параметров вы можете закончить одну из ветвей, уходящих от реального минимума. Оказавшись в долине, это может быть очень плоским, так что вы превысите максимальное количество итераций или получите что-то, что приемлемо в пределах допуска алгоритма. Как всегда, чем лучше ваши стартовые параметры, тем лучше. К сожалению, у вас есть только небольшая дуга круга, так что трудно поправиться. Я не специалист по Python, но думаю, что leastsq позволяет поиграться с методами Якобиана и Градиента. Попробуйте поиграть и с толерантностью. Вкратце: код мне в принципе кажется нормальным, но ваши данные патологические, и вам нужно адаптировать код к такого рода данным. Существует неитеративное решение в 2D от Karimäki, возможно, вы сможете адаптировать этим методом в 3D. Вы также можете посмотреть на это. Конечно, вы найдете больше литературы.

Я только что проверил данные с помощью симплексного алгоритма. Минимум, как я уже сказал, плохо себя ведет. См. здесь некоторые разрезы функции невязки. Только в xy-плоскости вы получаете разумное поведение. Свойства zr- и xr-плоскости делают процесс нахождения очень трудным.

введите здесь описание изображения

Итак, вначале симплекс-алгоритм находит несколько почти устойчивых решений. Вы можете увидеть их в виде плоских шагов на графике ниже (синий x, фиолетовый y, желтый z, зеленый R). В конце алгоритм должен пройти почти плоскую, но очень растянутую долину, что приводит к окончательному преобразованию z и R. Тем не менее, я ожидаю много областей, которые выглядят как решение, если допуск недостаточен. При стандартном допуске 10 ^ -5 алгоритм остановился примерно после 350 итераций. Мне пришлось установить его на 10 ^ -10, чтобы получить это решение, то есть [1899,32, 741,874, 298,696, 248,956], что выглядит вполне нормально.

введите здесь описание изображения

Обновить

Как упоминалось ранее, решение зависит от точности работы и требуемой точности. Таким образом, ваш ручной метод градиента, вероятно, работает лучше, поскольку эти значения отличаются от встроенного метода наименьших квадратов. Тем не менее, это моя версия, состоящая из двух шагов. Сначала я подгоняю самолет к данным. На следующем шаге я помещаю круг в эту плоскость. Оба шага используют метод наименьших квадратов. На этот раз это работает, так как каждый шаг позволяет избежать критических минимумов. (Естественно, с плоскостной подгонкой возникают проблемы, если сегмент дуги становится маленьким и данные лежат практически на прямой линии. Но это будет происходить для всех алгоритмов)

from math import *
from matplotlib import pyplot as plt
from scipy import optimize
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import pprint as pp

dataTupel=zip(xs,ys,zs) #your data from above

# Fitting a plane first
# let the affine plane be defined by two vectors, 
# the zero point P0 and the plane normal n0
# a point p is member of the plane if (p-p0).n0 = 0 

def distanceToPlane(p0,n0,p):
    return np.dot(np.array(n0),np.array(p)-np.array(p0))    

def residualsPlane(parameters,dataPoint):
    px,py,pz,theta,phi = parameters
    nx,ny,nz =sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta)
    distances = [distanceToPlane([px,py,pz],[nx,ny,nz],[x,y,z]) for x,y,z in dataPoint]
    return distances

estimate = [1900, 700, 335,0,0] # px,py,pz and zeta, phi
#you may automize this by using the center of mass data
# note that the normal vector is given in polar coordinates
bestFitValues, ier = optimize.leastsq(residualsPlane, estimate, args=(dataTupel))
xF,yF,zF,tF,pF = bestFitValues

point  = [xF,yF,zF]
normal = [sin(tF)*cos(pF),sin(tF)*sin(pF),cos(tF)]

# Fitting a circle inside the plane
#creating two inplane vectors
sArr=np.cross(np.array([1,0,0]),np.array(normal))#assuming that normal not parallel x!
sArr=sArr/np.linalg.norm(sArr)
rArr=np.cross(sArr,np.array(normal))
rArr=rArr/np.linalg.norm(rArr)#should be normalized already, but anyhow


def residualsCircle(parameters,dataPoint):
    r,s,Ri = parameters
    planePointArr = s*sArr + r*rArr + np.array(point)
    distance = [ np.linalg.norm( planePointArr-np.array([x,y,z])) for x,y,z in dataPoint]
    res = [(Ri-dist) for dist in distance]
    return res

estimateCircle = [0, 0, 335] # px,py,pz and zeta, phi
bestCircleFitValues, ier = optimize.leastsq(residualsCircle, estimateCircle,args=(dataTupel))

rF,sF,RiF = bestCircleFitValues
print bestCircleFitValues

# Synthetic Data
centerPointArr=sF*sArr + rF*rArr + np.array(point)
synthetic=[list(centerPointArr+ RiF*cos(phi)*rArr+RiF*sin(phi)*sArr) for phi in np.linspace(0, 2*pi,50)]
[cxTupel,cyTupel,czTupel]=[ x for x in zip(*synthetic)]

### Plotting
d = -np.dot(np.array(point),np.array(normal))# dot product
# create x,y mesh
xx, yy = np.meshgrid(np.linspace(2000,2200,10), np.linspace(540,740,10))
# calculate corresponding z
# Note: does not work if normal vector is without z-component
z = (-normal[0]*xx - normal[1]*yy - d)/normal[2]

# plot the surface, data, and synthetic circle
fig = plt.figure()
ax = fig.add_subplot(211, projection='3d')
ax.scatter(xs, ys, zs, c='b', marker='o')
ax.plot_wireframe(xx,yy,z)
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')
bx = fig.add_subplot(212, projection='3d')
bx.scatter(xs, ys, zs, c='b', marker='o')
bx.scatter(cxTupel,cyTupel,czTupel, c='r', marker='o')
bx.set_xlabel('X Label')
bx.set_ylabel('Y Label')
bx.set_zlabel('Z Label')
plt.show()

которые дают радиус 245. Это близко к тому, что дал другой подход (249). Таким образом, в пределах погрешности я получаю то же самое.

Полученный результат выглядит разумным. Надеюсь это поможет.

person mikuszefski    schedule 03.04.2013
comment
Что касается 2D-комментария от Hooked: это нормально, если вы можете свести свою проблему к 2D. Вы можете сделать это, сначала подобрав плоскость к своим данным и заставив центр находиться в этой плоскости. Подгонка плоскости будет использовать расстояние до плоскости в качестве невязок. Например: res=(p-p0).n0, где p — точка данных, p0 — точка плоскости и n0 — нормаль. p0 и n0 должны быть установлены. Здесь . скалярное векторное произведение. - person mikuszefski; 03.04.2013
comment
Спасибо за ваш ответ. Я думал так же. Алгоритм сделал сам, и сейчас он работает, хотя не могу сказать, чем он отличался от стандартного leastsq метода. Я также сейчас занят подгонкой самолета, поэтому я могу попробовать использовать его и для этой задачи и конвертировать в 2D. - person ; 08.04.2013
comment
Спасибо за обновление. Не все данные так хороши. Но это интересное решение для другой моей задачи, где все точки находятся примерно в плоскости по определению. - person ; 17.04.2013
comment
Я хотел бы использовать ваш код для своей задачи. Могу ли я это сделать, и как я должен ссылаться на вас в файле? Могу ли я отправить вам личное сообщение на этом сайте? - person ; 18.04.2013
comment
Привет, Александр, извини, что так долго не проверял. Конечно, вы можете использовать мой код, я сделал его общедоступным. В принципе, нет необходимости обращаться ко мне. Если хотите, можете связаться со мной в LinkedIn. Я обновил свою информацию LinkedIn до этого профиля пользователя. - person mikuszefski; 02.05.2013