Minimum Distance of a point to the boundary of an ellipse

ellipse.cpp

//
//  ellipse.cpp
//  EllipticalCylinder
//
//  Created by Yi Hu on 4/6/18.
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <nlopt.h>

struct opt_data_pack{
  double x, y, a, b;
};

// distance from a point (x, y) to another point (r, phi) on the ellipse
double distance2_point_to_ellipse_angle(unsigned n, const double *cx, double *grad, void *my_func_data)
{
  struct opt_data_pack *data = (struct opt_data_pack *)my_func_data;
  double xd = data->x - data->a*cos(*cx);
  double yd = data->y - data->b*sin(*cx);
  *grad = 2.0*(xd*data->a*sin(*cx) - yd*data->b*cos(*cx));
  return xd*xd+yd*yd;
}

// The minimum distance from a point to an ellipse
// point(x, y); ellipse: (x/a)^2+(y/b)^2=1
// if phi is NOT NULL, store the angular coordinate of the tangent point
double distance_point_to_ellipse(double x,double y,double a,double b, double* phi)
{
  double phi0, dsq, lb=0, ub=M_PI/2, tol=1e-12;
  double absx=fabs(x), absy=fabs(y);
  nlopt_opt optobj;
  nlopt_result optrst;
  struct opt_data_pack my_data = {.x=absx, .y=absy, .a=a, .b=b};
  if (x!=0.0)
  {
    phi0 = atan(absy/absx);
  }
  else
  {
    phi0 = M_PI/2;
  }
  optobj = nlopt_create(NLOPT_LD_MMA, 1);
  optrst = nlopt_set_min_objective(optobj, &distance2_point_to_ellipse_angle, &my_data);
  if(NLOPT_SUCCESS != optrst)
  {
    fprintf(stderr, "Error in set min objective. Error code: %d\n", optrst);
    return NAN;
  }
  nlopt_set_upper_bounds(optobj, &ub);
  nlopt_set_lower_bounds(optobj, &lb);
  nlopt_set_ftol_rel(optobj, tol);
  optrst = nlopt_optimize(optobj, &phi0, &dsq);
  if(NLOPT_SUCCESS != optrst && NLOPT_FTOL_REACHED != optrst)
  {
    if(optrst<0)
      fprintf(stderr, "Error in optimization. Error code: %d\n", optrst);
    else
      fprintf(stderr, "Warning in optimization. Warning code: %d\n", optrst);
  }
  nlopt_destroy(optobj);
  if(phi != NULL)
  {
    if(x<0)
      phi0 = M_PI - phi0;
    if(y<0)
      phi0 = 2.0*M_PI - phi0;
    *phi = phi0;
  }
  return sqrt(dsq);
}

double distance_point_to_ellipse(double x,double y,double a,double b)
{
  return distance_point_to_ellipse(x, y, a, b, NULL);
}