Pull to refresh

Как найти расстояние от точки к отрезку в трёхмерном пространстве, имея координаты трёх точек?

Небольшой дисклеймер

Статья рассчитана для начинающих (таких же как и я) программистов/математиков, которые столкнулись с такой проблемой - как же найти расстояние от точки к отрезку, имея их координаты? Да ещё и в трёхмерном пространстве! Подобной статьи, где есть чёткий ответ, мне не хватало, и информацию я собирал по крупицам, кое-что додумал сам. Сейчас же я хочу заполнить этот пробел. Ниже будет приведён сам алгоритм нахождения и примеры на Java. Так что, если нужен только код - сразу листайте вниз. Если хотите разобраться - добро пожаловать!

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

Постановка задачи

Предположим, что мы имеем три точки в трёхмерном пространстве. Точка B - начало отрезка с координатами (x1, y1, z1). Точка Е - конец отрезка с координатами (x2, y2, z2). Начало и конец, конечно, понятия более подходящие для вектора, уж простите меня. Соответсвенно, ВЕ - отрезок, к которому мы ищем расстояние. А точка Р(xp, yp, zp) - заданная точка, от которой и нужно найти расстояние к отрезку ВЕ. Ниже приведён рисунок одного из вероятных случаев расположения (далее будет подробней об этом):

Частные случаи расположения точки и отрезка в пространстве

1) Точка лежит на отрезке (в том числе, является одним из концов):

2) Точка лежит вне отрезка, но возможности провести перпендикуляр мы не имеем:

3) Точка лежит вне отрезка, и мы имеем возможность провести перпендикуляр:

4) Небольшой случай, который тоже нужно учесть - координаты начала и конца отрезков совпадают.

5) Все точки имеют одинаковые координаты, но данный случай перекроется в дальнейшей програмной реализации пунктом 1 и пунктом 4.

Начало програмной реализации

Хранить данные о координатах точки удобнее в классе. Поэтому создадим класс Point3D (поля оставлены публичными для удобства, но так делать не стоит):

class Point3D {
    public double x;
    public double y;
    public double z;

    Point3D(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    @Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;

        Point3D point3D = (Point3D) o;

        if (Double.compare(point3D.x, x) != 0) return false;
        if (Double.compare(point3D.y, y) != 0) return false;
        return Double.compare(point3D.z, z) == 0;
    }
}

Также создадим enum PointDispositionCase, хранящий возможные случаи:

enum PointDispositionCase {
    ON_SEGMENT,
    PERPENDICULAR_POSSIBLE,
    PERPENDICULAR_IMPOSSIBLE,
    BEGIN_AND_END_ARE_SAME,
}

А так же, сам класс для решения поставленной задачи - DistanceBetweenPointAndSegmentIn3D (далее просто - главный класс) и метод findDistance():

class DistanceBetweenPointAndSegmentIn3D {
    private Point3D P;
    private Point3D B;
    private Point3D E;

    DistanceBetweenPointAndSegmentIn3D(Point3D P, Point3D B, Point3D E) {
        this.P = P;
        this.B = B;
        this.E = E;
    }
    
    public double findDistance() {}
}

Подход к каждому случаю

Создадим в главном классе метод, который в зависимости от координат данных точек, будет возвращать соответсвующий член перечисления Cases:

static private PointDispositionCase findPointDispositionCase(
  Point3D P, Point3D B, Point3D E) {};

В первом случае определяем, принадлежит ли точка Р отрезку BE. Сделать это просто - если точка принадлежит отрезку, то сумма расстояния от начала отрезка к данной точке и расстояния от конца отрезка к данной точке равна длине самого отрезка. Если коротко, данное тождество должно быть верным:

BE = BP + PE

Напишем метод в главном классе, который находит длину отрезка:

static private double segmentLength(Point3D a, Point3D b) {
        return Math.sqrt(
                (b.x - a.x) * (b.x - a.x)
                        + (b.y - a.y) * (b.y - a.y)
                        + (b.z - a.z) * (b.z - a.z)
        );
    }

Не отходя от кассы, в методе findDespositionCase() найдём длину отрезков BE, BP и PE (это всё равно нам пригодится. и да, если рассматривать отрезки как вектора, эта формула справедлива и для них):

static private PointDispositionCase findPointDispositionCase(
  Point3D P, Point3D B, Point3D E) {
        double BE = segmentLength(B, E);
        double PE = segmentLength(P, E);
        double BP = segmentLength(B, P);
}

Если координаты одного из концов отрезка совпадают с заданной точкой Р, данная проверка учтёт и это. Если хотите, добавьте отдельную проверку, которая сверяет координаты заданных точек, но я ограничусь этой:

static private PointDispositionCase findPointDispositionCase(
  Point3D P, Point3D B, Point3D E) {
        double BE = segmentLength(B, E);
        double PE = segmentLength(P, E);
        double BP = segmentLength(B, P);

        if (BP + PE == BE) {
            return PointDispositionCase.ON_SEGMENT;
        }
}

Во втором случае поработаем немного с косинусами. Для этого придётся рассмотреть треугольник ВРЕ и некоторые вектора. Тут мы имеем три подслучая:

1) Угол между заданной точкой и одним из концов отрезка тупой. Соответсвенно, его косинус меньше 0. Соответсвенно, если это угол PBE, то искомая длина - отрезок PB, а если PEB - то PE.

2) Угол между заданной точкой и одним из концов отрезка прямой. Соответственно, cos = 0. Аналогично предыдущему пункту, если это угол PBE, то искомая длина - отрезок PB, а если PEB - то PE.

3) Углы между заданной точкой и обоими концами отрезка острые. Это переадресовывает нас к третьему случаю, о котором дальше.

Повторюсь, для нахождения косинуса угла нужно рассматривать стороны треугольника РВЕ как вектора. Создадим класс Vector:

class Vector {
    double x;
    double y;
    double z;

    Vector (Point3D a, Point3D b) {
        this.x = b.x - a.x;
        this.y = b.y - a.y;
        this.z = b.z - a.z;
    }

    public double vectorLength() {
        return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z);
    }
}

Теперь в главный класс добавим метод для нахождения косинуса между векторами:

static private double findCos(Vector a, Vector b) {
        return (a.x * b.x + a.y * b.y + a.z * b.z) /
                (a.vectorLength() + b.vectorLength());
}

И обновим проверку (случаи с прямым и тупым углом объединим):

static private PointDispositionCase findPointDispositionCase(Point3D P, Point3D B, Point3D E) {
        double BE = segmentLength(B, E);
        double PE = segmentLength(P, E);
        double BP = segmentLength(B, P);

        if (BP + PE == BE) {
            return PointDispositionCase.ON_SEGMENT;
        }

        double cosPEB = findCos(
                new Vector(E, P),
                new Vector(E, B)
        );

        double cosPBE = findCos(
                new Vector(B, P),
                new Vector(B, E)
        );

        if (cosPBE <= 0 || cosPEB <= 0) {
            return PointDispositionCase.PERPENDICULAR_IMPOSSIBLE;
        }
}

Рассмотрим теперь третий случай. Имеем треугольник BEP. Пускай BН - это перпендикуляр к BP (а значит, искомое расстояние). Вспоминаем, что есть формула вычисления площади треугольника через высоту: S = (a * ha) / 2

А так как мы имеем длины всех сторон треугольника, то эту же площадь можно найти по формуле Герона: S = sqrt(p * (p - a) * (p - b) * (p - c))

Соответсвенно, с предыдущих двух формул выводим формулу для нахождения высоты:

ha = 2 / a * Math.sqrt(p * (p - a) * (p - b) * (p - c))

Это значение и будет расстоянием к отрезку в данном случае. А теперь в главном классе создадим метод для нахождения высоты через площадь:

static private double findDistanceByTriangleSquare(
  Point3D P, Point3D B, Point3D E) {
        double a = segmentLength(B, E);
        double b = segmentLength(B, P);
        double c = segmentLength(E, P);
        double p = (a + b + c) / 2;

        return 2 / a * Math.sqrt(p * (p - a) * (p - b) * (p - c));
}

Дело за малым - четвёртый случай. Можем либо сравнить координаты точек, либо просто найти расстояние отрезка. Если оно равно 0, то начало и конец совпадают. Добавим этот и предыдущий случаи в метод проверки:

static private PointDispositionCase findPointDispositionCase(Point3D P, Point3D B, Point3D E) {
        double BE = segmentLength(B, E);
        double PE = segmentLength(P, E);
        double BP = segmentLength(B, P);

        if (BP + PE == BE) {
            return PointDispositionCase.ON_SEGMENT;
        }

        if (B.equals(E)) {
            return PointDispositionCase.BEGIN_AND_END_ARE_SAME;
        }

        double cosPEB = findCos(
                new Vector(E, P),
                new Vector(E, B)
        );

        double cosPBE = findCos(
                new Vector(B, P),
                new Vector(B, E)
        );

        if (cosPBE <= 0 || cosPEB <= 0) {
            return PointDispositionCase.PERPENDICULAR_IMPOSSIBLE;
        }

        return PointDispositionCase.PERPENDICULAR_POSSIBLE;
 }

Теперь всё, что нам осталось - это написать метод в главном классе, который в соответсвии с положением точки вернёт расстояние от точки к отрезку:

static public double findDistance(Point3D P, Point3D B, Point3D E) {
        PointDispositionCase pointDispositionCase =
                findPointDispositionCase(P, B , E);

        switch (pointDispositionCase) {
            case ON_SEGMENT:
                return 0;
            case BEGIN_AND_END_ARE_SAME:
                return segmentLength(P, B);
            case PERPENDICULAR_IMPOSSIBLE:
                return Math.min(
                        segmentLength(P, E),
                        segmentLength(P, B)
                );
            case PERPENDICULAR_POSSIBLE:
                return findDistanceByTriangleSquare(P, B, E);
        }

        return -1;
    }

Финальный результат (Java)

Вот, что у нас получилось. На идеал не претендую, поправите под свои нужды) Мне в своё время хотя бы такой статьи не хватало.

import java.util.Scanner;

class Vector {
    double x;
    double y;
    double z;

    Vector (Point3D a, Point3D b) {
        this.x = b.x - a.x;
        this.y = b.y - a.y;
        this.z = b.z - a.z;
    }

    public double vectorLength() {
        return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z);
    }
}

class Point3D {
    public double x;
    public double y;
    public double z;

    Point3D(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    @Override
    public boolean equals(Object o) {
        if (this == o) return true;
        if (o == null || getClass() != o.getClass()) return false;

        Point3D point3D = (Point3D) o;

        if (Double.compare(point3D.x, x) != 0) return false;
        if (Double.compare(point3D.y, y) != 0) return false;
        return Double.compare(point3D.z, z) == 0;
    }
}

enum PointDispositionCase {
    ON_SEGMENT,
    PERPENDICULAR_POSSIBLE,
    PERPENDICULAR_IMPOSSIBLE,
    BEGIN_AND_END_ARE_SAME,
}

class DistanceBetweenPointAndSegmentIn3D {
    static public double findDistance(Point3D P, Point3D B, Point3D E) {
        PointDispositionCase pointDispositionCase =
                findPointDispositionCase(P, B , E);

        switch (pointDispositionCase) {
            case ON_SEGMENT:
                return 0;
            case BEGIN_AND_END_ARE_SAME:
                return segmentLength(P, B);
            case PERPENDICULAR_IMPOSSIBLE:
                return Math.min(
                        segmentLength(P, E),
                        segmentLength(P, B)
                );
            case PERPENDICULAR_POSSIBLE:
                return findDistanceByTriangleSquare(P, B, E);
        }

        return -1;
    }

    static private PointDispositionCase findPointDispositionCase(Point3D P, Point3D B, Point3D E) {
        double BE = segmentLength(B, E);
        double PE = segmentLength(P, E);
        double BP = segmentLength(B, P);

        if (BP + PE == BE) {
            return PointDispositionCase.ON_SEGMENT;
        }

        if (B.equals(E)) {
            return PointDispositionCase.BEGIN_AND_END_ARE_SAME;
        }

        double cosPEB = findCos(
                new Vector(E, P),
                new Vector(E, B)
        );

        double cosPBE = findCos(
                new Vector(B, P),
                new Vector(B, E)
        );

        if (cosPBE <= 0 || cosPEB <= 0) {
            return PointDispositionCase.PERPENDICULAR_IMPOSSIBLE;
        }

        return PointDispositionCase.PERPENDICULAR_POSSIBLE;
    }

    static private double segmentLength(Point3D a, Point3D b) {
        return Math.sqrt(
                (b.x - a.x) * (b.x - a.x)
                        + (b.y - a.y) * (b.y - a.y)
                        + (b.z - a.z) * (b.z - a.z)
        );
    }

    static private double findCos(Vector a, Vector b) {
        return (a.x * b.x + a.y * b.y + a.z * b.z) /
                (a.vectorLength() + b.vectorLength());
    }

    static private double findDistanceByTriangleSquare(Point3D P, Point3D B, Point3D E) {
        double a = segmentLength(B, E);
        double b = segmentLength(B, P);
        double c = segmentLength(E, P);
        double p = (a + b + c) / 2;

        return 2 / a * Math.sqrt(p * (p - a) * (p - b) * (p - c));
    }
}

Спасибо всем, кто уделил время данной статье! Надеюсь, кому-то это было полезно и сэкономило хотя бы чуточку времени. Удачи!

Tags:
Hubs:
You can’t comment this publication because its author is not yet a full member of the community. You will be able to contact the author only after he or she has been invited by someone in the community. Until then, author’s username will be hidden by an alias.