Как стать автором
Обновить

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

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

Статья рассчитана для начинающих (таких же как и я) программистов/математиков, которые столкнулись с такой проблемой - как же найти расстояние от точки к отрезку, имея их координаты? Да ещё и в трёхмерном пространстве! Подобной статьи, где есть чёткий ответ, мне не хватало, и информацию я собирал по крупицам, кое-что додумал сам. Сейчас же я хочу заполнить этот пробел. Ниже будет приведён сам алгоритм нахождения и примеры на 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));
    }
}

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

Теги:
Хабы:
Данная статья не подлежит комментированию, поскольку её автор ещё не является полноправным участником сообщества. Вы сможете связаться с автором только после того, как он получит приглашение от кого-либо из участников сообщества. До этого момента его username будет скрыт псевдонимом.