Slicer  5.1
Slicer is a multi-platform, free and open source software package for visualization and medical image computing
coordTypes.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Extract Skeleton
4  Module: $HeadURL$
5  Language: C++
6  Date: $Date$
7  Version: $Revision$
8 
9  Copyright (c) Brigham and Women's Hospital (BWH) All Rights Reserved.
10 
11  See License.txt or http://www.slicer.org/copyright/copyright.txt for details.
12 
13 ==========================================================================*/
14 
15 #ifndef COORD_TYPES_H
16 #define COORD_TYPES_H
17 
18 #include "math.h"
19 
20 #include "misc.h"
21 
22 class Coord3i
23 {
24  int v[3];
25 public:
27  {
28  v[0] = v[1] = v[2] = -1;
29  }
30  inline int & operator[](int i)
31  {
32  return v[i];
33  }
34  inline void conv(double * i)
35  {
36  i[0] = v[0]; i[1] = v[1]; i[2] = v[2];
37  };
38 };
39 
40 class Coord3f
41 {
42  float v[3];
43 public:
44  inline float & operator[](int i)
45  {
46  return v[i];
47  }
48  inline void conv(float * i)
49  {
50  i[0] = v[0]; i[1] = v[1]; i[2] = v[2];
51  };
52  inline void conv(double * i)
53  {
54  i[0] = v[0]; i[1] = v[1]; i[2] = v[2];
55  };
56 };
57 
58 class Coord3d
59 {
60  double v[3];
61 public:
62  inline double & operator[](int i)
63  {
64  return v[i];
65  };
66  inline void conv(int * i)
67  {
68  i[0] = (int) v[0]; i[1] = (int) v[1]; i[2] = (int) v[2];
69  };
70  inline void conv(float * i)
71  {
72  i[0] = static_cast<float>(v[0]); i[1] = static_cast<float>(v[1]); i[2] = static_cast<float>(v[2]);
73  };
74  inline void conv(double * i)
75  {
76  i[0] = v[0]; i[1] = v[1]; i[2] = v[2];
77  };
78 };
79 
80 // Euclidean distance between two points
81 // TODO: use image pixel spacing
82 inline double pointdistance(Coord3i &p1, Coord3i &p2, const double spacing[3])
83 {
84  return sqrt(
85  sqr((p1[0] - p2[0]) * spacing[0])
86  + sqr((p1[1] - p2[1]) * spacing[1])
87  + sqr((p1[2] - p2[2]) * spacing[2])
88  );
89 }
90 
91 inline void normcrossprod(double *v1, double *v2, double *norm)
92 // calculate normalized crossproduct
93 {
94  double absval;
95 
96  norm[0] = v1[1] * v2[2] - v1[2] * v2[1];
97  norm[1] = v1[2] * v2[0] - v1[0] * v2[2];
98  norm[2] = v1[0] * v2[1] - v1[1] * v2[0];
99 
100  absval = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
101  norm[0] /= absval;
102  norm[1] /= absval;
103  norm[2] /= absval;
104 }
105 
106 // calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
107 // range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
108 inline double vectorangle(double *v1, double *v2)
109 {
110  double prod = 0, length1 = 0, length2 = 0;
111 
112  for( int k = 0; k < 3; k++ )
113  {
114  prod += v1[k] * v2[k];
115  length1 += v1[k] * v1[k];
116  length2 += v2[k] * v2[k];
117  }
118  return acos(prod / sqrt(length1 * length2) );
119 }
120 
121 // calculate angle between two vectors (0..M_PI), you might want to adjust to 0..M_PI/2
122 // range after call via 'if (angle > M_PI/2) angle = M_PI-angle;'
123 inline double vectorangle(Coord3d v1, Coord3d v2)
124 {
125  double prod = 0, length1 = 0, length2 = 0;
126 
127  for( int k = 0; k < 3; k++ )
128  {
129  prod += v1[k] * v2[k];
130  length1 += v1[k] * v1[k];
131  length2 += v2[k] * v2[k];
132  }
133  return acos(prod / sqrt(length1 * length2) );
134 }
135 
136 inline double vec_length(Coord3d x)
137 {
138  return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]) );
139 }
140 
141 inline double vec_length(double *x)
142 {
143  return sqrt(sqr(x[0]) + sqr(x[1]) + sqr(x[2]) );
144 }
145 
146 inline double vec_length(double *x, double *y)
147 {
148  return sqrt(sqr(x[0] - y[0]) + sqr(x[1] - y[1]) + sqr(x[2] - y[2]) );
149 }
150 
151 // transform and check index, returns 0 on success and 1 if corrected location
152 inline int transWorldToImage(Coord3d loc_world, int *loc_img,
153  double *origin, int *dims, double voxelsize)
154 {
155  int adjust = 0;
156 
157  for( int i = 0; i < 3; i++ )
158  {
159  loc_img[i] = (int) ( (loc_world[i] - origin[i]) / voxelsize);
160  if( loc_img[i] < 0 )
161  {
162  adjust = 1; loc_img[i] = 0;
163  }
164  if( loc_img[i] >= dims[i] )
165  {
166  loc_img[i] = dims[i] - 1; adjust = 1;
167  }
168  }
169 
170  return adjust;
171 }
172 
173 inline int transWorldToImage(double *loc_world, int *loc_img,
174  double *origin, int *dims, double voxelsize)
175 // transform and check index, returns 0 on success and 1 if corrected location
176 {
177  int adjust = 0;
178 
179  for( int i = 0; i < 3; i++ )
180  {
181  loc_img[i] = (int) ( (loc_world[i] - origin[i]) / voxelsize);
182  if( loc_img[i] < 0 )
183  {
184  adjust = 1; loc_img[i] = 0;
185  }
186  if( loc_img[i] >= dims[i] )
187  {
188  loc_img[i] = dims[i] - 1; adjust = 1;
189  }
190  }
191 
192  return adjust;
193 }
194 
195 #endif
void conv(double *i)
Definition: coordTypes.h:34
double pointdistance(Coord3i &p1, Coord3i &p2, const double spacing[3])
Definition: coordTypes.h:82
void conv(float *i)
Definition: coordTypes.h:48
void conv(int *i)
Definition: coordTypes.h:66
double vectorangle(double *v1, double *v2)
Definition: coordTypes.h:108
double vec_length(Coord3d x)
Definition: coordTypes.h:136
void conv(double *i)
Definition: coordTypes.h:74
double & operator[](int i)
Definition: coordTypes.h:62
void conv(double *i)
Definition: coordTypes.h:52
void normcrossprod(double *v1, double *v2, double *norm)
Definition: coordTypes.h:91
int & operator[](int i)
Definition: coordTypes.h:30
int transWorldToImage(Coord3d loc_world, int *loc_img, double *origin, int *dims, double voxelsize)
Definition: coordTypes.h:152
Coord3i()
Definition: coordTypes.h:26
void conv(float *i)
Definition: coordTypes.h:70
float & operator[](int i)
Definition: coordTypes.h:44
T sqr(T x)
Definition: misc.h:64