1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkVector_hxx
19 #define itkVector_hxx
20
21 #include "itkMath.h"
22 #include "vnl/vnl_vector.h"
23 #include "itkObject.h"
24 #include "itkNumericTraitsVectorPixel.h"
25
26 namespace itk
27 {
28 template< typename T, unsigned int TVectorDimension >
29 Vector< T, TVectorDimension >
Vector(const ValueType & r)30 ::Vector(const ValueType & r)
31 {
32 for ( typename BaseArray::Iterator i = BaseArray::Begin(); i != BaseArray::End(); ++i )
33 {
34 *i = r;
35 }
36 }
37
38 template< typename T, unsigned int TVectorDimension >
39 Vector< T, TVectorDimension > &
40 Vector< T, TVectorDimension >
operator =(const ValueType r[TVectorDimension])41 ::operator=(const ValueType r[TVectorDimension])
42 {
43 BaseArray::operator=(r);
44 return *this;
45 }
46
47 template< typename T, unsigned int TVectorDimension >
48 const typename Vector< T, TVectorDimension >::Self &
49 Vector< T, TVectorDimension >
operator +=(const Self & vec)50 ::operator+=(const Self & vec)
51 {
52 for ( unsigned int i = 0; i < TVectorDimension; i++ )
53 {
54 ( *this )[i] += vec[i];
55 }
56 return *this;
57 }
58
59 template< typename T, unsigned int TVectorDimension >
60 const typename Vector< T, TVectorDimension >::Self &
61 Vector< T, TVectorDimension >
operator -=(const Self & vec)62 ::operator-=(const Self & vec)
63 {
64 for ( unsigned int i = 0; i < TVectorDimension; i++ )
65 {
66 ( *this )[i] -= vec[i];
67 }
68 return *this;
69 }
70
71 template< typename T, unsigned int TVectorDimension >
72 Vector< T, TVectorDimension >
73 Vector< T, TVectorDimension >
operator -() const74 ::operator-() const
75 {
76 Self result;
77
78 for ( unsigned int i = 0; i < TVectorDimension; i++ )
79 {
80 result[i] = -( *this )[i];
81 }
82 return result;
83 }
84
85 template< typename T, unsigned int TVectorDimension >
86 Vector< T, TVectorDimension >
87 Vector< T, TVectorDimension >
operator +(const Self & vec) const88 ::operator+(const Self & vec) const
89 {
90 Self result;
91
92 for ( unsigned int i = 0; i < TVectorDimension; i++ )
93 {
94 result[i] = ( *this )[i] + vec[i];
95 }
96 return result;
97 }
98
99 template< typename T, unsigned int TVectorDimension >
100 Vector< T, TVectorDimension >
101 Vector< T, TVectorDimension >
operator -(const Self & vec) const102 ::operator-(const Self & vec) const
103 {
104 Self result;
105
106 for ( unsigned int i = 0; i < TVectorDimension; i++ )
107 {
108 result[i] = ( *this )[i] - vec[i];
109 }
110 return result;
111 }
112
113 template< typename T, unsigned int TVectorDimension >
114 typename Vector< T, TVectorDimension >::RealValueType
115 Vector< T, TVectorDimension >
GetSquaredNorm() const116 ::GetSquaredNorm() const
117 {
118 typename NumericTraits< RealValueType >::AccumulateType sum = NumericTraits< T >::ZeroValue();
119 for ( unsigned int i = 0; i < TVectorDimension; i++ )
120 {
121 const RealValueType value = ( *this )[i];
122 sum += value * value;
123 }
124 return sum;
125 }
126
127 template< typename T, unsigned int TVectorDimension >
128 typename Vector< T, TVectorDimension >::RealValueType
129 Vector< T, TVectorDimension >
GetNorm() const130 ::GetNorm() const
131 {
132 return RealValueType( std::sqrt( double( this->GetSquaredNorm() ) ) );
133 }
134
135 template< typename T, unsigned int TVectorDimension >
136 typename Vector< T, TVectorDimension >::RealValueType
137 Vector< T, TVectorDimension >
Normalize()138 ::Normalize()
139 {
140 const RealValueType norm = this->GetNorm();
141 if (norm < NumericTraits< RealValueType >::epsilon())
142 {
143 return norm; // Prevent division by 0
144 }
145
146 const RealValueType inversedNorm = 1.0 / norm;
147 for ( unsigned int i = 0; i < TVectorDimension; i++ )
148 {
149 ( *this )[i] =
150 static_cast< T >( static_cast< RealValueType >( ( *this )[i] * inversedNorm ));
151 }
152 return norm;
153 }
154
155 template< typename T, unsigned int TVectorDimension >
156 vnl_vector_ref< T >
157 Vector< T, TVectorDimension >
GetVnlVector()158 ::GetVnlVector()
159 {
160 return vnl_vector_ref< T >( TVectorDimension, this->GetDataPointer() );
161 }
162
163 template< typename T, unsigned int TVectorDimension >
164 vnl_vector< T >
165 Vector< T, TVectorDimension >
GetVnlVector() const166 ::GetVnlVector() const
167 {
168 // Return a vector_ref<>. This will be automatically converted to a
169 // vnl_vector<>. We have to use a const_cast<> which would normally
170 // be prohibited in a const method, but it is safe to do here
171 // because the cast to vnl_vector<> will ultimately copy the data.
172 return vnl_vector_ref< T >( TVectorDimension,
173 const_cast< T * >( this->GetDataPointer() ) );
174 }
175
176 template< typename T, unsigned int TVectorDimension >
177 void
178 Vector< T, TVectorDimension >
SetVnlVector(const vnl_vector<T> & v)179 ::SetVnlVector(const vnl_vector< T > & v)
180 {
181 for ( unsigned int i = 0; i < v.size(); i++ )
182 {
183 ( *this )[i] = v(i);
184 }
185 }
186
187 template< typename T, unsigned int TVectorDimension >
188 std::ostream &
operator <<(std::ostream & os,const Vector<T,TVectorDimension> & vct)189 operator<<(std::ostream & os, const Vector< T, TVectorDimension > & vct)
190 {
191 os << "[";
192 if ( TVectorDimension == 1 )
193 {
194 os << vct[0];
195 }
196 else
197 {
198 for ( unsigned int i = 0; i + 1 < TVectorDimension; i++ )
199 {
200 os << vct[i] << ", ";
201 }
202 os << vct[TVectorDimension - 1];
203 }
204 os << "]";
205 return os;
206 }
207
208 template< typename T, unsigned int TVectorDimension >
209 std::istream &
operator >>(std::istream & is,Vector<T,TVectorDimension> & vct)210 operator>>(std::istream & is, Vector< T, TVectorDimension > & vct)
211 {
212 for ( unsigned int i = 0; i < TVectorDimension; i++ )
213 {
214 is >> vct[i];
215 }
216 return is;
217 }
218
219 template< typename T, unsigned int TVectorDimension >
220 typename Vector< T, TVectorDimension >::ValueType
221 Vector< T, TVectorDimension >
operator *(const Self & other) const222 ::operator*(const Self & other) const
223 {
224 typename NumericTraits< T >::AccumulateType value = NumericTraits< T >::ZeroValue();
225 for ( unsigned int i = 0; i < TVectorDimension; i++ )
226 {
227 value += ( *this )[i] * other[i];
228 }
229 return value;
230 }
231
232 } // end namespace itk
233 #endif
234