coordTypes.h

Go to the documentation of this file.

/*=========================================================================
 
  Program:   Extract Skeleton
  Module:    $HeadURL$
  Language:  C++
  Date:      $Date$
  Version:   $Revision$
 
  Copyright (c) Brigham and Women's Hospital (BWH) All Rights Reserved.
 
  See License.txt or http://www.slicer.org/copyright/copyright.txt for details.
 
==========================================================================*/
 
#ifndef COORD_TYPES_H
#define COORD_TYPES_H
 
#include "math.h"
 
#include "misc.h"
 
class Coord3i
{
  int v[3];
 
public:
  Coord3i() { v[0] = v[1] = v[2] = -1; }
  inline int& operator[](int i) { return v[i]; }
  inline void conv(double* i)
  {
    i[0] = v[0];
    i[1] = v[1];
    i[2] = v[2];
  };
};
 
class Coord3f
{
  float v[3];
 
public:
  inline float& operator[](int i) { return v[i]; }
  inline void conv(float* i)
  {
    i[0] = v[0];
    i[1] = v[1];
    i[2] = v[2];
  };
  inline void conv(double* i)
  {
    i[0] = v[0];
    i[1] = v[1];
    i[2] = v[2];
  };
};
 
class Coord3d
{
  double v[3];
 
public:
  inline double& operator[](int i) { return v[i]; };
  inline void conv(int* i)
  {
    i[0] = (int)v[0];
    i[1] = (int)v[1];
    i[2] = (int)v[2];
  };
  inline void conv(float* i)
  {
    i[0] = static_cast<float>(v[0]);
    i[1] = static_cast<float>(v[1]);
    i[2] = static_cast<float>(v[2]);
  };
  inline void conv(double* i)
  {
    i[0] = v[0];
    i[1] = v[1];
    i[2] = v[2];
  };
};
 
// Euclidean distance between two points
// TODO: use image pixel spacing
inline double pointdistance(Coord3i& p1, Coord3i& p2, const double spacing[3])
{
  return sqrt(sqr((p1[0] - p2[0]) * spacing[0]) + sqr((p1[1] - p2[1]) * spacing[1]) + sqr((p1[2] - p2[2]) * spacing[2]));
}
 
inline void normcrossprod(double* v1, double* v2, double* norm)
// calculate normalized crossproduct
{
  double absval;
 
  norm[0] = v1[1] * v2[2] - v1[2] * v2[1];
  norm[1] = v1[2] * v2[0] - v1[0] * v2[2];
  norm[2] = v1[0] * v2[1] - v1[1] * v2[0];
 
  absval = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
  norm[0] /= absval;
  norm[1] /= absval;
  norm[2] /= absval;
}
 
// calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
// range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
inline double vectorangle(double* v1, double* v2)
{
  double prod = 0, length1 = 0, length2 = 0;
 
  for (int k = 0; k < 3; k++)
  {
    prod += v1[k] * v2[k];
    length1 += v1[k] * v1[k];
    length2 += v2[k] * v2[k];
  }
  return acos(prod / sqrt(length1 * length2));
}
 
// calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
// range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
inline double vectorangle(Coord3d v1, Coord3d v2)
{
  double prod = 0, length1 = 0, length2 = 0;
 
  for (int k = 0; k < 3; k++)
  {
    prod += v1[k] * v2[k];
    length1 += v1[k] * v1[k];
    length2 += v2[k] * v2[k];
  }
  return acos(prod / sqrt(length1 * length2));
}
 
inline double vec_length(Coord3d x)
{
  return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]));
}
 
inline double vec_length(double* x)
{
  return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]));
}
 
inline double vec_length(double* x, double* y)
{
  return sqrt(sqr(x[0] - y[0]) + sqr(x[1] - y[1]) + sqr(x[2] - y[2]));
}
 
// transform and check index, returns 0 on success and 1 if corrected location
inline int transWorldToImage(Coord3d loc_world, int* loc_img, double* origin, int* dims, double voxelsize)
{
  int adjust = 0;
 
  for (int i = 0; i < 3; i++)
  {
    loc_img[i] = (int)((loc_world[i] - origin[i]) / voxelsize);
    if (loc_img[i] < 0)
    {
      adjust = 1;
      loc_img[i] = 0;
    }
    if (loc_img[i] >= dims[i])
    {
      loc_img[i] = dims[i] - 1;
      adjust = 1;
    }
  }
 
  return adjust;
}
 
inline int transWorldToImage(double* loc_world, int* loc_img, double* origin, int* dims, double voxelsize)
// transform and check index, returns 0 on success and 1 if corrected location
{
  int adjust = 0;
 
  for (int i = 0; i < 3; i++)
  {
    loc_img[i] = (int)((loc_world[i] - origin[i]) / voxelsize);
    if (loc_img[i] < 0)
    {
      adjust = 1;
      loc_img[i] = 0;
    }
    if (loc_img[i] >= dims[i])
    {
      loc_img[i] = dims[i] - 1;
      adjust = 1;
    }
  }
 
  return adjust;
}
 
#endif