Кривые точки крепления в трехмерном пространстве - PullRequest
8 голосов
/ 06 декабря 2010

Попытка найти функции, которые помогут нам провести 3D-линию через серию точек.

Для каждой точки, которую мы знаем: дата и время, широта, долгота, высота, скорость и курс. Данные могут записываться каждые 10 секунд, и мы хотели бы иметь возможность угадывать точки между ними и увеличивать детализацию до 1 секунды. Таким образом создавая виртуальную траекторию полета в трехмерном пространстве.

Я нашел несколько алгоритмов подбора кривой, которые будут аппроксимировать линию через ряд точек, но они не гарантируют, что точки пересекаются. Они также не учитывают скорость и направление движения, чтобы определить наиболее вероятный путь, по которому объект достигнет следующей точки.

Ответы [ 5 ]

19 голосов
/ 06 декабря 2010

С точки зрения физики:

Вы должны предположить что-то относительно ускорения в ваших промежуточных точках, чтобы получить интерполяцию.

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

Векторное уравнение для постоянно меняющегося ускоренного движения:

 x''[t] = a t + b

где все величины, кроме t, являются векторами.

Для каждого сегмента, который вы уже знаете, v (t = t0) x (t = t0) tfinal и x (tfinal) v (tfinal)

Решая дифференциальное уравнение, вы получаете:

Eq 1:
x[t_] := (3 b t^2 Tf + a t^3 Tf - 3 b t Tf^2 - a t Tf^3 - 6 t X0 + 6 Tf X0 + 6 t Xf)/(6 Tf)  

И установив начальные и конечные ограничения для положения и скорости, вы получите:

Уравнения 2:

 a -> (6 (Tf^2 V0 - 2 T0 Tf Vf + Tf^2 Vf - 2 T0 X0 + 2 Tf X0 + 
        2 T0 Xf - 2 Tf Xf))/(Tf^2 (3 T0^2 - 4 T0 Tf + Tf^2))

 b -> (2 (-2 Tf^3 V0 + 3 T0^2 Tf Vf - Tf^3 Vf + 3 T0^2 X0 - 
        3 Tf^2 X0 - 3 T0^2 Xf + 3 Tf^2 Xf))/(Tf^2 (3 T0^2 - 4 T0 Tf + Tf^2))}}

Таким образом, вставляя значения для уравнений 2 в уравнение 1, вы получаете временную интерполяцию для ваших точек на основе начальной и конечной позиции и скоростей.

НТН!

Редактировать

Несколько примеров с резким изменением скорости в двух измерениях (в 3D точно то же самое). Если начальная и конечная скорости одинаковы, вы получите «более прямые» пути.

Предположим:

X0 = {0, 0}; Xf = {1, 1};
T0 = 0;      Tf = 1;  

Если

V0 = {0, 1}; Vf = {-1, 3};

alt text

V0 = {0, 1}; Vf = {-1, 5};

alt text

V0 = {0, 1}; Vf = {1, 3};

alt text

Вот анимация, где вы можете увидеть изменение скорости от V0 = {0, 1} до Vf = {1, 5}: alt text

Здесь вы можете увидеть ускоряющееся тело в 3D с позициями, взятыми с равными интервалами:

alt text

Редактировать

Полная проблема:

Для удобства я буду работать в декартовых координатах. Если вы хотите конвертировать из lat / log / alt в декартову, просто выполните:

x = rho sin(theta) cos(phi)
y = rho sin(theta) sin(phi)
z = rho cos(theta)

Где phi - долгота, theta - широта, а rho - ваша высота плюс радиус Земли.

Итак, предположим, что мы начинаем наш сегмент с:

 t=0 with coordinates (0,0,0) and velocity (1,0,0)

и заканчиваются на

 t=10 with coordinates (10,10,10) and velocity (0,0,1)  

Я четко внес изменение в начало координат, чтобы установить начало координат в моей начальной точке. Это просто для получения хороших круглых чисел ...

Таким образом, мы заменяем эти числа в формулах для a и b и получаем:

a = {-(3/50), -(3/25), -(3/50)}  b = {1/5, 3/5, 2/5}  

С этими значениями мы переходим к уравнению 1, и положение объекта определяется по формуле:

p[t] = {1/60 (60 t + 6 t^2 - (3 t^3)/5), 
        1/60 (18 t^2 - (6 t^3)/5), 
        1/60 (12 t^2 - (3 t^3)/5)}

И это все. Вы получаете позицию от 1 до 10 секунд, заменяя t на ее значение в приведенном выше уравнении.
Анимация запускается:

alt text

Редактировать 2

Если вы не хотите связываться с вертикальным ускорением (возможно, из-за того, что ваш «спидометр» его не читает), вы можете просто назначить постоянную скорость оси z (есть очень незначительная ошибка, если учесть ее параллельно оси Rho), равным (Zfinal - Zinit) / (Tf-T0), а затем решаем задачу в плоскости, забывая высоту.

1 голос
/ 06 декабря 2010

Что вам нужно (вместо того, чтобы моделировать физику), это разместить сплайн в данных. Я использовал книгу числовых получателей (http://www.nrbook.com/a имеет бесплатные алгоритмы C и FORTRAN. Посмотрите раздел F77 3.3, чтобы получить необходимую математику). Если вы хотите быть простым, тогда просто проведите линии через точки, но это не приведет к гладкой траектории полета вообще. Время будет вашим значением x, и каждый зарегистрированный параметр будет иметь свои собственные общедоступные параметры сплайна.

Поскольку нам нравятся длинные сообщения по этому вопросу, вот полный код:

// драйвер программы

static void Main(string[] args)
{
    double[][] flight_data = new double[][] {
        new double[] { 0, 10, 20, 30 }, // time in seconds
        new double[] { 14500, 14750, 15000, 15125 }, //altitude in ft
        new double[] { 440, 425, 415, 410 }, // speed in knots
    };

    CubicSpline altitude = new CubicSpline(flight_data[0], flight_data[1]);
    CubicSpline speed = new CubicSpline(flight_data[0], flight_data[2]);

    double t = 22; //Find values at t
    double h = altitude.GetY(t);
    double v = speed.GetY(t);

    double ascent = altitude.GetYp(t); // ascent rate in ft/s

}

// Определение кубического сплайна

using System.Linq;

/// <summary>
///     Cubic spline interpolation for tabular data
/// </summary>
/// <remarks>
///     Adapted from numerical recipies for FORTRAN 77 
///     (ISBN 0-521-43064-X), page 110, section 3.3.
///     Function spline(x,y,yp1,ypn,y2) converted to 
///     C# by jalexiou, 27 November 2007.
///     Spline integration added also Nov 2007.
/// </remarks>
public class CubicSpline 
{
    double[] xi;
    double[] yi;
    double[] yp;
    double[] ypp;
    double[] yppp;
    double[] iy;

    #region Constructors

    public CubicSpline(double x_min, double x_max, double[] y)
        : this(Sequence(x_min, x_max, y.Length), y)
    { }

    public CubicSpline(double x_min, double x_max, double[] y, double yp1, double ypn)
        : this(Sequence(x_min, x_max, y.Length), y, yp1, ypn)
    { }

    public CubicSpline(double[] x, double[] y)
        : this(x, y, double.NaN, double.NaN)
    { }

    public CubicSpline(double[] x, double[] y, double yp1, double ypn)
    {
        if( x.Length == y.Length )
        {
            int N = x.Length;
            xi = new double[N];
            yi = new double[N];
            x.CopyTo(xi, 0);
            y.CopyTo(yi, 0);

            if( N > 0 )
            {
                double p, qn, sig, un;
                ypp = new double[N];
                double[] u = new double[N];

                if( double.IsNaN(yp1) )
                {
                    ypp[0] = 0;
                    u[0] = 0;
                }
                else
                {
                    ypp[0] = -0.5;
                    u[0] = (3 / (xi[1] - xi[0])) * 
                          ((yi[1] - yi[0]) / (x[1] - x[0]) - yp1);
                }
                for (int i = 1; i < N-1; i++)
                {
                    double hp = x[i] - x[i - 1];
                    double hn = x[i + 1] - x[i];
                    sig = hp / hn;
                    p = sig * ypp[i - 1] + 2.0;
                    ypp[i] = (sig - 1.0) / p;
                    u[i] = (6 * ((y[i + 1] - y[i]) / hn) - (y[i] - y[i - 1]) / hp)
                        / (hp + hn) - sig * u[i - 1] / p;
                }
                if( double.IsNaN(ypn) )
                {
                    qn = 0;
                    un = 0;
                }
                else
                {
                    qn = 0.5;
                    un = (3 / (x[N - 1] - x[N - 2])) * 
                          (ypn - (y[N - 1] - y[N - 2]) / (x[N - 1] - x[N - 2]));
                }
                ypp[N - 1] = (un - qn * u[N - 2]) / (qn * ypp[N - 2] + 1.0);
                for (int k = N-2; k > 0; k--)
                {
                    ypp[k] = ypp[k] * ypp[k + 1] + u[k];
                }

                // Calculate 1st derivatives
                yp = new double[N];
                double h;
                for( int i = 0; i < N - 1; i++ )
                {
                    h = xi[i + 1] - xi[i];
                    yp[i] = (yi[i + 1] - yi[i]) / h
                           - h / 6 * (ypp[i + 1] + 2 * ypp[i]);
                }
                h = xi[N - 1] - xi[N - 2];
                yp[N - 1] = (yi[N - 1] - yi[N - 2]) / h
                             + h / 6 * (2 * ypp[N - 1] + ypp[N - 2]);
                // Calculate 3rd derivatives as average of dYpp/dx
                yppp = new double[N];
                double[] jerk_ij = new double[N - 1];
                for( int i = 0; i < N - 1; i++ )
                {
                    h = xi[i + 1] - xi[i];
                    jerk_ij[i] = (ypp[i + 1] - ypp[i]) / h;
                }
                Yppp = new double[N];
                yppp[0] = jerk_ij[0];
                for( int i = 1; i < N - 1; i++ )
                {
                    yppp[i] = 0.5 * (jerk_ij[i - 1] + jerk_ij[i]);
                }
                yppp[N - 1] = jerk_ij[N - 2];
                // Calculate Integral over areas
                iy = new double[N];
                yi[0] = 0;  //Integration constant
                for( int i = 0; i < N - 1; i++ )
                {
                    h = xi[i + 1] - xi[i];
                    iy[i + 1] = h * (yi[i + 1] + yi[i]) / 2
                             - h * h * h / 24 * (ypp[i + 1] + ypp[i]);
                }
            }
            else
            {
                yp = new double[0];
                ypp = new double[0];
                yppp = new double[0];
                iy = new double[0];
            }
        }
        else
            throw new IndexOutOfRangeException();
    }

    #endregion

    #region Actions/Functions
    public int IndexOf(double x)
    {
        //Use bisection to find index
        int i1 = -1;
        int i2 = Xi.Length;
        int im;
        double x1 = Xi[0];
        double xn = Xi[Xi.Length - 1];
        bool ascending = (xn >= x1);
        while( i2 - i1 > 1 )
        {
            im = (i1 + i2) / 2;
            double xm = Xi[im];
            if( ascending & (x >= Xi[im]) ) { i1 = im; } else { i2 = im; }
        }
        if( (ascending && (x <= x1)) || (!ascending & (x >= x1)) )
        {
            return 0;
        }
        else if( (ascending && (x >= xn)) || (!ascending && (x <= xn)) )
        {
            return Xi.Length - 1;
        }
        else
        {
            return i1;
        }
    }

    public double GetIntY(double x)
    {
        int i = IndexOf(x);
        double x1 = xi[i];
        double x2 = xi[i + 1];
        double y1 = yi[i];
        double y2 = yi[i + 1];
        double y1pp = ypp[i];
        double y2pp = ypp[i + 1];
        double h = x2 - x1;
        double h2 = h * h;
        double a = (x-x1)/h;
        double a2 = a*a;
        return h / 6 * (3 * a * (2 - a) * y1 
                 + 3 * a2 * y2 - a2 * h2 * (a2 - 4 * a + 4) / 4 * y1pp
                 + a2 * h2 * (a2 - 2) / 4 * y2pp);
    }


    public double GetY(double x)
    {
        int i = IndexOf(x);
        double x1 = xi[i];
        double x2 = xi[i + 1];
        double y1 = yi[i];
        double y2 = yi[i + 1];
        double y1pp = ypp[i];
        double y2pp = ypp[i + 1];
        double h = x2 - x1;
        double h2 = h * h;
        double A = 1 - (x - x1) / (x2 - x1);
        double B = 1 - A;

        return A * y1 + B * y2 + h2 / 6 * (A * (A * A - 1) * y1pp
                     + B * (B * B - 1) * y2pp);
    }

    public double GetYp(double x)
    {
        int i = IndexOf(x);
        double x1 = xi[i];
        double x2 = xi[i + 1];
        double y1 = yi[i];
        double y2 = yi[i + 1];
        double y1pp = ypp[i];
        double y2pp = ypp[i + 1];
        double h = x2 - x1;
        double A = 1 - (x - x1) / (x2 - x1);
        double B = 1 - A;

        return (y2 - y1) / h + h / 6 * (y2pp * (3 * B * B - 1)
                   - y1pp * (3 * A * A - 1));
    }

    public double GetYpp(double x)
    {
        int i = IndexOf(x);
        double x1 = xi[i];
        double x2 = xi[i + 1];
        double y1pp = ypp[i];
        double y2pp = ypp[i + 1];
        double h = x2 - x1;
        double A = 1 - (x - x1) / (x2 - x1);
        double B = 1 - A;

        return A * y1pp + B * y2pp;
    }

    public double GetYppp(double x)
    {
        int i = IndexOf(x);
        double x1 = xi[i];
        double x2 = xi[i + 1];
        double y1pp = ypp[i];
        double y2pp = ypp[i + 1];
        double h = x2 - x1;

        return (y2pp - y1pp) / h;
    }

    public double Integrate(double from_x, double to_x)
    {
        if( to_x < from_x ) { return -Integrate(to_x, from_x); }
        int i = IndexOf(from_x);
        int j = IndexOf(to_x);
        double x1 = xi[i];
        double xn = xi[j];
        double z = GetIntY(to_x) - GetIntY(from_x);    // go to nearest nodes (k) (j)
        for( int k = i + 1; k <= j; k++ )
        {
            z += iy[k];                             // fill-in areas in-between
        }            
        return z;
    }


    #endregion

    #region Properties
    public bool IsEmpty { get { return xi.Length == 0; } }
    public double[] Xi { get { return xi; } set { xi = value; } }
    public double[] Yi { get { return yi; } set { yi = value; } }
    public double[] Yp { get { return yp; } set { yp = value; } }
    public double[] Ypp { get { return ypp; } set { ypp = value; } }
    public double[] Yppp { get { return yppp; } set { yppp = value; } }
    public double[] IntY { get { return yp; } set { iy = value; } }
    public int Count { get { return xi.Length; } }
    public double X_min { get { return xi.Min(); } }
    public double X_max { get { return xi.Max(); } }
    public double Y_min { get { return yi.Min(); } }
    public double Y_max { get { return yi.Max(); } }
    #endregion

    #region Helpers

    static double[] Sequence(double x_min, double x_max, int double_of_points)
    {
        double[] res = new double[double_of_points];
        for (int i = 0; i < double_of_points; i++)
        {
            res[i] = x_min + (double)i / (double)(double_of_points - 1) * (x_max - x_min);
        }
        return res;
    }

    #endregion
}
1 голос
/ 06 декабря 2010

То, что вы спрашиваете, является общей интерполяцией проблемы. Полагаю, ваша настоящая проблема не в том, что используется алгоритм подбора кривой, а в том, что вы применяете его ко всем дискретным значениям, записанным системой, а не к соответствующему набору значений.

Давайте разберем вашу проблему. В настоящее время вы рисуете точку в сферически отображаемом трехмерном пространстве, корректируя линейные и изогнутые контуры. Если мы дискретизируем операции, выполняемые объектом с шестью степенями свободы ( roll, pitch и yaw ), то единственные операции, которые вас особенно интересуют, - это линейные и изогнутые пути, учитывающие высоту и отклонение в любое направление. Учет ускорения и замедления также возможен при понимании базовой физики.

Работать со сферическим отображением легко. Просто разверните ваши точки относительно их положения на плоскости, корректируя широту, долготу и высоту. Это должно позволить вам сгладить данные, которые в противном случае могли бы существовать по изогнутой траектории, хотя это не может быть строго необходимым для решения вашей проблемы (см. Ниже).

Линейная интерполяция проста. Учитывая произвольное количество точек назад во времени, которые соответствуют линии в пределах ошибки n , как определено вашей системой, * построить линию и вычислить расстояние во времени между каждой точкой. Отсюда попытка подогнать временные точки к одному из двух случаев: постоянная скорость или постоянное ускорение.

Кривая интерполяция немного сложнее, но все же правдоподобно. Для случаев тангажа, рыскания или комбинированного тангажа и рыскания создайте плоскость, содержащую произвольное количество точек назад во времени, в пределах ошибки m для изогнутых показаний из вашей системы. * Из этих данных постройте плоская кривая и еще раз учитывается постоянная скорость или ускорение по кривой .

Вы можете добиться большего, чем это, пытаясь предсказать ожидаемые операции самолета в полете как часть дерева решений или нейронной сети относительно траектории полета. Я оставлю это как упражнение для читателя.

Удачи в разработке вашей системы.

-

* Ожидается, что оба показания ошибки будут получены из данных GPS, учитывая описание проблемы. Учет и корректировка ошибок в этих данных - отдельная интересная проблема.

0 голосов
/ 06 декабря 2010

Я предполагаю, что серьезное применение этого будет использовать http://en.wikipedia.org/wiki/Kalman_filter. Кстати, это, вероятно, не будет гарантировать, что сообщаемые точки также были пересечены, если только вы не поиграли с параметрами.Он будет ожидать некоторой степени ошибки в каждой заданной точке данных, поэтому, если он считает, что объект находится в момент времени T, он не обязательно будет там, где он был в момент времени T. Конечно, вы можете установить распределение ошибок, чтобы сказать, что вы былиАбсолютно уверен, что вы знали, где это было в момент времени T.

Если не использовать фильтр Калмана, я бы попытался превратить его в проблему оптимизации.Работайте с гранулярностью 1 с и запишите уравнения вроде x_t '= x_t + (Vx_t + Vx_t') / 2 + e,

Vx_t_reported = Vx_t + f,

Vx_t '= Vx_t + gгде e, f и g представляют шум.Затем создайте штрафную функцию, такую ​​как e ^ 2 + f ^ 2 + g ^ 2 + ... или некоторую взвешенную версию, такую ​​как 1.5e ^ 2 + 3f ^ 2 + 2.6g ^ 2 + ... в соответствии с вашей идеейкакие ошибки на самом деле и насколько гладко вы хотите получить ответ, и найдите значения, которые делают функцию штрафа как можно меньше - с наименьшими квадратами, если уравнения получаются хорошими.

0 голосов
/ 06 декабря 2010

Вы можете найти аппроксимацию линии, которая пересекает точки в 3d и 2d пространстве, используя алгоритм Преобразование Хафа . Я только знаком с его использованием в 2d, но он все равно будет работать для трехмерных пространств, учитывая, что вы знаете, какую линию вы ищете. Существует базовое описание реализации. Вы можете использовать Google для готовых изделий, и вот ссылка на реализацию 2d C CImg .

Алгоритм процесса (примерно) ... Сначала вы найдете уравнение линии, которое, по вашему мнению, будет наилучшим образом приближено к форме линии (в 2-й параболической, логарифмической, экспоненциальной и т. Д.). Вы берете эту формулу и решаете один из параметров.

y = ax + b

становится

b = y - ax

Далее, для каждой точки, которую вы пытаетесь сопоставить, вы подключаете точки к значениям y и x. С 3 баллами у вас будет 3 отдельных функции b по отношению к a.

(2, 3)  : b =  3 -  2a
(4, 1)  : b =  1 -  4a
(10, -5): b = -5 - 10a

Далее, теория состоит в том, что вы находите все возможные линии, которые проходят через каждую из точек, что является бесконечным числом для каждой отдельной точки, однако при объединении в пространство аккумулятора только несколько возможных параметров лучше всего подходят. На практике это делается путем выбора диапазона значений для параметров (я выбрал -2 <= a <= 1, 1 <= b <= 6) и начинаю подключать значения для параметра (ов) варианта и решать для другого , Вы подсчитываете количество пересечений каждой функции в аккумуляторе. Точки с самыми высокими значениями дают вам ваши параметры. </p>

Аккумулятор после обработки b = 3 - 2a

    a: -2   -1    0    1
b: 1                   
   2                  
   3              1   
   4                  
   5         1        
   6                  

аккумулятор после обработки b = 1 - 4a

    a: -2   -1    0    1
b: 1              1     
   2                  
   3              1   
   4                  
   4                  
   5         2        
   6                  

Аккумулятор после обработки b = -5 - 10a

    a: -2   -1    0    1
b: 1              1     
   2                  
   3              1   
   4                  
   5         3        
   6                  

Набор параметров с наибольшим накопленным значением - (b a) = (5 -1), а функция, наиболее подходящая для заданных точек, - y = 5 - x.

Удачи.

...