3D Line-of-Sight algorithm

Question

I'm currently implementing a line-of-sight algorithm that will tell me the points I can and cannot see along a line. So, if I am standing atop some hilly terrain, I can know where I can and cannot see.

The line that is generated from point A to point B includes a number of points evenly spaced between A and B. Starting from A, I see what the angle of elevation is between A and B. I note that as my alphaAngle.

Next, for each point between A and B, I get the angle of elevation between A and that point. If that point is the highest angle of elevation so far(excluding alphaAngle), then I mark it as such. Otherwise, it has a lower angle of elevation and therefor I should not be able to see it, and mark that point as having hasLOS = false.

Here are some object definitions I use:

struct TerrainProfilPnt
{
 double m_x_lon; //lon
 double m_y_lat; //lat
 double m_z_elev; //elev
 bool m_hasLOS; //Does this point have line of sight from another point?
};

class TerrainProfilePartitioner; // Holds a collection of points that make up the line

Here is my algorithm written out, however it doesn't return the correct results. Either it claims it has LOS when it shouldn't (like going from behind one hill to the opposite side of another hill, I shouldn't be able to see that). Or it claims that I cannot see the end point when I should (Top of a hill to the valley beneath it). So, I am suspecting that either my implementation of finding the line-of-sight is incorrect, or I am implementing it wrong in code.

using Point = TerrainProfilePnt;

    auto pointsInLine = terrainProfilePartitioner->GetPoints();
    auto& pointsVec = pointsInLine->GetVector();
    std::vector<Point> terrainProfileVec;
    terrainProfileVec.reserve(pointsVec.size());

    double start_xlon = 0.0f;
    double start_ylat = 0.0f;
    double start_zelev = 0.0f;

    double end_xlon = 0.0f;
    double end_ylat = 0.0f;
    double end_zelev = 0.0f;

    //The angle between the starting point and the ending point
    double initialElevAngle = 0.0f;


    //Assemble the first and last points
    start_xlon = pointsVec.front().x();
    start_ylat = pointsVec.front().y();
  GetPointElevation(start_xlon, start_ylat, start_zelev);    
    end_xlon = pointsVec.back().x();
    end_ylat = pointsVec.back().y();
    GetPointElevation(end_xlon, end_ylat, end_zelev);

    //Calculate the angle between the beginning and ending points
    initialElevAngle = atan2(start_zelev, end_zelev) * 180 / M_PI;

    Point initialPoint;
    initialPoint.m_x_lon = start_xlon;
    initialPoint.m_y_lat = start_ylat;
    initialPoint.m_z_elev = start_zelev;
    initialPoint.m_hasLOS = true;

    //pointsVec.push_back(initialPoint);
    terrainProfileVec.push_back(initialPoint);

    double oldAngle = 0.0f;;
    bool hasOldAngle = false;
    for (std::size_t i = 1; i < pointsVec.size(); ++i)
    {
        Point p;
        p.m_x_lon = pointsVec.at(i).x();
        p.m_y_lat = pointsVec.at(i).y();
        GetPointElevation(p.m_x_lon, p.m_y_lat, p.m_z_elev);

        double newAngle = atan2(start_zelev, p.m_z_elev) * 180 / M_PI;
        if (!hasOldAngle)
        {
            hasOldAngle = true;
            oldAngle = newAngle;
            p.m_hasLOS = true;
        }
        else
        {
            if (newAngle < oldAngle)
            {
                p.m_hasLOS = false;
            }
            else
            { 
                oldAngle = newAngle;
                p.m_hasLOS = true;
            }
        }

    }

    auto terrainPrfileSeq = new Common::ObjectRandomAccessSequence<TerrainProfilePnt>(terrainProfileVec);

    return terrainPrfile

score 2 · Accepted Answer · answered Oct 05 '16 at 20:26

2

atan2(start_zelev, p.m_z_elev) is meaningless. You need

    atan2(distance, p.m_z_elev - start_zelev);

where distance is horizontal distance between p and initialPoint.

answered Oct 05 '16 at 20:26

user58697

7,808
1
14
28

score 1 · Answer 2 · edited May 23 '17 at 12:00

I see few problems: You are not sorting the points by distance from A ascending and use only elevation. Also using goniometrics for each map location which is probably slow. I would use different approach instead:

set LOS=false for whole map
write a DDA line interpolation
```
P(t)=A+(B-A).t
```
where (B-A) is your viewing direction and t is parameter incrementing by step smaller then your grid size (so you do not miss anything). This would give you long,lat without any slow goniometrics.

Then find the object on your map. If it is topology sorted then this is either O(1) or O(log(n).log(m))
test each actual location P(t)

for visibility. If true set LOS=true and loop #2 otherwise stop the whole thing.

btw this is sometimes called ray-casting.

3D Line-of-Sight algorithm

2 Answers2