cppcode/miscIPOL.cpp

Go to the documentation of this file.

/*
 * Copyright (c) 2009-2010 Jose-Luis Lisani <joseluis.lisani@uib.es>
 * All rights reserved.
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>



void normalize_vector(float u[3])
{
  float d=sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
  u[0]/=d;
  u[1]/=d;
  u[2]/=d;
}

float scalar_product(float u[3], float v[3])
{
  float p=u[0]*v[0]+u[1]*v[1]+u[2]*v[2];
  return p;
}

void vectorial_product(float u[3], float v[3], float w[3])
{
  w[0]=u[1]*v[2]-u[2]*v[1];
  w[1]=u[2]*v[0]-u[0]*v[2];
  w[2]=u[0]*v[1]-u[1]*v[0];
}

float vector_norm(float u[3])
{
  float d=sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
  return d;
}

float **sample_arclength(float A[3], float B[3], float O[3], int npoints)
{
  float **P;
  float u[3], v[3], vv[3], w[3];
  float unorm, vvnorm, R, x, y;
  float angmax, stepang;

  if (npoints == 1) return NULL;

  //initialize list
  P=new float*[npoints];
  for (int n=0; n < npoints; n++) P[n]=new float[3];

  //u: vector in the direction OA
  u[0]=A[0]-O[0];
  u[1]=A[1]-O[1];
  u[2]=A[2]-O[2];
  unorm=vector_norm(u);
  //vv: vector in the direction OB
  vv[0]=B[0]-O[0];
  vv[1]=B[1]-O[1];
  vv[2]=B[2]-O[2];
  vvnorm=vector_norm(vv);
  //Since OA and OB belong to the same circunference centered at O, 
  //unorm=vvnorm=R (radius of circunference)
  R=unorm;

  //angle from OA to OB
  angmax=acos((double) (scalar_product(u, vv)/(unorm*vvnorm)));

  //w: normal vector to the plane containing A-O-B
  normalize_vector(u);
  normalize_vector(vv);
  vectorial_product(u, vv, w);
  normalize_vector(w);
  //v: vector on the plane, normal to u
  vectorial_product(w, u, v);
  normalize_vector(v);
  //u-v define an orthogonal coordinates system on the plane A-O-B

  //angular step
  stepang=angmax/(npoints-1);

  //get samples
  //first point
  P[0][0]=A[0];
  P[0][1]=A[1];
  P[0][2]=A[2];
  for (int n=1; n < npoints-1; n++) {
    x=R*cos((double) (n*stepang));
    y=R*sin((double) (n*stepang));
    P[n][0]=O[0]+x*u[0]+y*v[0];
    P[n][1]=O[1]+x*u[1]+y*v[1];
    P[n][2]=O[2]+x*u[2]+y*v[2];
  }
  //last point
  P[npoints-1][0]=B[0];
  P[npoints-1][1]=B[1];
  P[npoints-1][2]=B[2];

  //return list
  return P;
}

void delete_list3D(float **points, int npoints)
{
  for (int n=0; n < npoints; n++) delete[] points[n];
  delete[] points;
}

---