point to segment distance program having problems

[GonlyG]

problem: given N amount of x1,y1, x2, y2 and xp, yp, try to calculate perpendicular distance between (xp, yp) and the line formed by the x's and the y's.

For some reason my output is wrong by small amounts. Any help is appreciated.

  #include <math.h>
#include <iostream>
#include <iomanip>
using namespace std;
int main() {
int inputs;
cin>>inputs;
for (int i=0;i<inputs;i++){
double x1, y1, x2, y2, xp, yp;
cin>>x1>>y1>>x2>>y2>>xp>>yp;

if(x1>x2){
    double temp = x1; x1=x2; x2=temp;
    temp = y2; y1=y2; y2= temp;
}
double x =0.0; double y =0.0; //default
if (y2!=y1 && x2!=x1){
    //    slope formula
    double m = (y2-y1)/(x2-x1);
    double mp = (-1.0)/m;            //perpendicular lines 
    // y = mx + c, 
       double cs = y2 - m*x2;
       double cp = yp - mp*xp;
       
       x = (cp - cs)/(m - mp);
       y = (m*cp - mp*cs)/(m-mp);
}
else if (x2==x1){
    y = yp;
    x = x1;
}
else{
    x =xp;
    y= y1;
}

 if(x <= x2 && x >= x1){
        double d = sqrt((xp - x)*(xp - x) + (yp - y)*(yp - y));
        cout<<setprecision(12)<<d<<" ";}
    else{
        double d1 = sqrt((x1 - xp)*(x1 - xp) + (y1 - yp)*(y1 - yp));
        double d2 =  sqrt((x2 - xp)*(x2 - xp) + (y2 - yp)*(y2 - yp));
        if(d1<d2){
         cout<<setprecision(12)<<d1<<" ";}
         else{
         cout<<setprecision(12)<<d2<<" ";
         }
         
    }
         
}
    
}

sample input:
19
12 7 1 3 0 14
13 17 14 12 17 17
3 15 1 19 4 1
19 8 17 16 10 4
3 0 18 2 18 19
17 19 5 18 6 19
15 15 19 8 2 2
9 10 19 10 6 3
14 7 17 11 9 9
10 6 15 8 2 15
6 5 20 11 9 0
13 17 17 12 2 1
15 14 13 14 4 2
2 18 9 20 4 1
14 19 18 8 9 12
9 12 12 8 6 7
10 17 11 6 3 15
9 3 10 1 8 16
3 19 12 12 17 14

output: 11.0453610172 3.92232270276 18.0277563773 13.8924439894 17 1 17.7369670463 7.61577310586 5.38516480713 12.0415945788 5.77748304583 18.5846776647 15 17.1172427686 7.09116855905 5.4 7.15232393736 13.0384048104 5.38516480713

expected: 10.6794707215 3.92232270276 14.0356688476 9.8488578018 17 0.913500278391 17.7369670463 7.61577310586 5.38516480713 12.0415945788 5.77748304583 18.5846776647 15 17.1172427686 7.09116855905 5.4 7.15232393736 13.0384048104 5.38516480713

[lastchance]

Lines 41 to 47 aren't finding a perpendicular distance - they are finding the minimum distance to the truncated line segment ... which will be larger if the projection lies outside the two end points. You will see that, for the values you get wrong, you are always computing too high. If the projection lies between the two ends then you get it right.

If you want a perpendicular distance then use the magnitude of the vector product of displacement relative to one end point with a unit vector along the line. This will even work in 3d and is far quicker than finding the projection first.

[GonlyG]

hold on, my bad, I didn't include the important specific of the program:
"if the other end of the perpendicular belongs to our segment (i.e. lays between its ends), then the distance to this segment equals to distance to this line;
otherwise distance from point to segment equals to distance to nearest of its ends." (https://www.codeabbey.com/index/task_view/point-to-segment-distance)

does 41~47 seem correct for this? Someone used a similar algorithm in python and it worked.

[lastchance]

Your points swapping is wrong. Fix the first statement on line 14.

You can also find the projection point using the dot product (which also works in 3-d) rather than the obscure version with slopes and point swapping (which doesn't work in 3-d).

[keskiverto]

Dot product, as in: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation

Have you had structs/classes and operator overloading yet?
They would be nice here.

[lastchance]

#include <iostream>
#include <sstream>
#include <cmath>
using namespace std;

//======================================================================

struct Pt { double x, y, z; };

Pt operator + ( Pt p, Pt q )     { return { p.x + q.x, p.y + q.y, p.z + q.z }; }
Pt operator - ( Pt p, Pt q )     { return { p.x - q.x, p.y - q.y, p.z - q.z }; }
Pt operator * ( double r, Pt p ) { return {   r * p.x,   r * p.y,   r * p.z }; }
Pt operator / ( Pt p, double r ) { return {   p.x / r,   p.y / r,   p.z / r }; }
double dot( Pt p, Pt q ) { return p.x * q.x + p.y * q.y + p.z * q.z; }
Pt cross( Pt p, Pt q ) { return { p.y * q.z - p.z * q.y, p.z * q.x - p.x * q.z, p.x * q.y - p.y * q.x }; }
double abs( Pt p ) { return sqrt( dot( p, p ) ); }
ostream &operator << ( ostream &out, const Pt &p ) { return out << p.x << '\t' << p.y << '\t' << p.z; }
istream &operator >> ( istream &in ,       Pt &p ) { return in  >> p.x         >> p.y         >> p.z; }
bool getPt2d         ( istream &in ,       Pt &p ) { p.z = 0.0;   return (bool)( in >> p.x >> p.y );  }

//======================================================================

double perpendicularDistance( Pt p1, Pt p2, Pt p )
{
   return abs( cross( p - p1, p2 - p1 ) ) / abs( p2 - p1 );   // distance is  r sin( angle )
}

//======================================================================

double distanceToLineSegment( Pt p1, Pt p2, Pt p )
{
   double segLength = abs( p2 - p1 );
   if ( segLength < 1.0e-30 ) return abs( p - p1 );           // paranoia
   double d = dot( p - p1, p2 - p1 ) / segLength;             // signed distance to projection
   if      ( d < 0.0       ) return abs( p - p1 );            // projection is to "left" of p1
   else if ( d > segLength ) return abs( p - p2 );            // projection is to "right" of p2
   else                      return abs( p - ( p1 + d * ( p2 - p1 ) / segLength ) );    // distance to projection
}

//======================================================================

int main()
{
   stringstream in( "12 7 1 3 0 14     \n"
                    "13 17 14 12 17 17 \n"
                    "3 15 1 19 4 1     \n"
                    "19 8 17 16 10 4   \n"
                    "3 0 18 2 18 19    \n"
                    "17 19 5 18 6 19   \n" );

   cout << "Distance to segment  vs  perpendicular distance ...\n";
   for ( Pt p1, p2, p; getPt2d( in, p1 ) && getPt2d( in, p2 ) && getPt2d( in, p ); )
      cout << distanceToLineSegment( p1, p2, p ) << "\t\t" << perpendicularDistance( p1, p2, p ) << '\n';
}

Distance to segment  vs  perpendicular distance ...
10.6795		10.6795
3.92232		3.92232
14.0357		5.36656
9.84886		9.70143
17		16.8509
0.9135		0.9135

[keskiverto]

You were less paranoid in perpendicular case. ;)

By "nice" I obviously mean that the "algorithm" part of code, for example:

double perpendicularDistance( Pt p1, Pt p2, Pt p )
{
   return abs( cross( p - p1, p2 - p1 ) ) / abs( p2 - p1 );   // distance is  r sin( angle )
}

remains compact, clear and as similar to math equations as possible.
Less clutter, less typos.


The downside is that one had to write lines 8-19 for that. Then again, each function there is small and specific (i.e. easy to make correct), and this block of code could be (re)used in all programs that operate on "points".

[lastchance]

keskiverto wrote:

You were less paranoid in perpendicular case. ;)

Ah, that's true! Not much error-checking in my codes.

kesiverto wrote:

as similar to math equations as possible.

Mmm, I really wanted to define operators rather than call functions dot() and cross() for those products, but sadly c++ doesn't allow you to define your own names or symbols and none of the allowable set look right. Since I tend to write a caret (^) for a cross product I'm probably out of luck anyway.

this block of code could be (re)used in all programs that operate on "points".

I confess! It was.