00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkJointDomainImageToListAdaptor_h
00018 #define __itkJointDomainImageToListAdaptor_h
00019
00020 #include "itkMacro.h"
00021 #include "itkFixedArray.h"
00022 #include "itkPoint.h"
00023 #include "itkPixelTraits.h"
00024 #include "itkImageToListAdaptor.h"
00025 #include "itkImageRegionConstIteratorWithIndex.h"
00026 #include "itkEuclideanDistance.h"
00027 #include "itkListSample.h"
00028
00029 namespace itk {
00030 namespace Statistics {
00031
00039 template< class TImage >
00040 struct ImageJointDomainTraits
00041 {
00042 typedef ImageJointDomainTraits Self;
00043 typedef PixelTraits< typename TImage::PixelType > PixelTraitsType;
00044 typedef typename PixelTraitsType::ValueType RangeDomainMeasurementType;
00045 typedef typename TImage::IndexType::IndexValueType IndexValueType;
00046
00047 itkStaticConstMacro(ImageDimension, unsigned int, TImage::ImageDimension);
00048 itkStaticConstMacro(Dimension,
00049 unsigned int,
00050 TImage::ImageDimension +
00051 PixelTraitsType::Dimension );
00052
00053 typedef float CoordinateRepType;
00054 typedef Point< CoordinateRepType, itkGetStaticConstMacro(ImageDimension) >
00055 PointType;
00056 typedef JoinTraits< RangeDomainMeasurementType, CoordinateRepType >
00057 JoinTraitsType;
00058 typedef typename JoinTraitsType::ValueType MeasurementType;
00059 typedef FixedArray< MeasurementType, itkGetStaticConstMacro(Dimension) >
00060 MeasurementVectorType;
00061 };
00062
00092 template < class TImage >
00093 class ITK_EXPORT JointDomainImageToListAdaptor
00094 : public ImageToListAdaptor<
00095 TImage,
00096 typename ImageJointDomainTraits< TImage >::MeasurementVectorType >
00097 {
00098 public:
00099 typedef ImageJointDomainTraits< TImage > ImageJointDomainTraitsType;
00100 typedef typename ImageJointDomainTraitsType::MeasurementVectorType
00101 MeasurementVectorType;
00102 typedef typename ImageJointDomainTraitsType::MeasurementType
00103 MeasurementType;
00104 typedef typename ImageJointDomainTraitsType::RangeDomainMeasurementType
00105 RangeDomainMeasurementType;
00106 typedef typename ImageJointDomainTraitsType::PointType
00107 PointType;
00108 typedef typename ImageJointDomainTraitsType::CoordinateRepType
00109 CoordinateRepType;
00111 typedef JointDomainImageToListAdaptor Self;
00112 typedef ImageToListAdaptor< TImage, MeasurementVectorType >
00113 Superclass;
00114 typedef SmartPointer< Self > Pointer;
00115 typedef SmartPointer<const Self> ConstPointer;
00116
00118 itkTypeMacro(JointDomainImageToListAdaptor, ImageToListAdaptor);
00119
00121 itkNewMacro(Self);
00122
00124 itkStaticConstMacro(MeasurementVectorSize,
00125 unsigned int,
00126 ImageJointDomainTraitsType::Dimension);
00127
00128 typedef typename Superclass::MeasurementVectorSizeType MeasurementVectorSizeType;
00129
00130 virtual void SetMeasurementVectorSize( const MeasurementVectorSizeType s )
00131 {
00132
00133
00134 if( s != MeasurementVectorSize )
00135 {
00136 itkExceptionMacro( << "Cannot set measurement vector size of "
00137 << " JointDomainImageToListAdaptor to " << s );
00138 }
00139 }
00140
00141 MeasurementVectorSizeType GetMeasurementVectorSize() const
00142 {
00143 return MeasurementVectorSize;
00144 }
00145
00148 typedef typename Superclass::FrequencyType FrequencyType;
00149 typedef typename Superclass::InstanceIdentifier InstanceIdentifier;
00150
00151 typedef typename TImage::IndexType ImageIndexType;
00152 typedef typename TImage::IndexType::IndexValueType ImageIndexValueType;
00153 typedef typename TImage::SizeType ImageSizeType;
00154 typedef typename TImage::RegionType ImageRegionType;
00155 typedef ImageRegionConstIteratorWithIndex< TImage > ImageIteratorType;
00156
00157 typedef MeasurementVectorType ValueType;
00158
00159 itkStaticConstMacro(RangeDomainDimension,
00160 unsigned int,
00161 itk::PixelTraits<
00162 typename TImage::PixelType >::Dimension);
00163
00164 typedef FixedArray< RangeDomainMeasurementType,
00165 itkGetStaticConstMacro( RangeDomainDimension ) >
00166 RangeDomainMeasurementVectorType;
00167
00168 typedef std::vector< InstanceIdentifier > InstanceIdentifierVectorType;
00169
00170 typedef FixedArray< float, itkGetStaticConstMacro(MeasurementVectorSize) >
00171 NormalizationFactorsType;
00172
00174 void SetNormalizationFactors(NormalizationFactorsType& factors);
00175
00178 inline const MeasurementVectorType & GetMeasurementVector(const InstanceIdentifier &id) const;
00179
00182 inline void ComputeRegion(const MeasurementVectorType& mv,
00183 const double radius,
00184 ImageRegionType& region) const;
00185
00189 inline void Search(const MeasurementVectorType& mv,
00190 const double radius,
00191 InstanceIdentifierVectorType& result) const;
00192
00193 protected:
00194 JointDomainImageToListAdaptor();
00195 virtual ~JointDomainImageToListAdaptor() {}
00196 void PrintSelf(std::ostream& os, Indent indent) const;
00197
00198 private:
00199 JointDomainImageToListAdaptor(const Self&);
00200 void operator=(const Self&);
00201
00202 NormalizationFactorsType m_NormalizationFactors;
00203
00204 mutable MeasurementVectorType m_TempVector;
00205 mutable PointType m_TempPoint;
00206 mutable ImageIndexType m_TempIndex;
00207 mutable RangeDomainMeasurementVectorType m_TempRangeVector;
00208 };
00209
00210 template < class TImage >
00211 inline const typename JointDomainImageToListAdaptor< TImage >::MeasurementVectorType &
00212 JointDomainImageToListAdaptor< TImage >
00213 ::GetMeasurementVector(const InstanceIdentifier &id) const
00214 {
00215 m_TempIndex = this->GetImage()->ComputeIndex( id );
00216
00217 this->GetImage()->TransformIndexToPhysicalPoint( m_TempIndex, m_TempPoint );
00218
00219 for ( unsigned int i = 0; i < TImage::ImageDimension; ++i )
00220 {
00221 m_TempVector[i] = m_TempPoint[i] / m_NormalizationFactors[i];
00222 }
00223
00224 if( this->GetUseBuffer() )
00225 {
00226 m_TempRangeVector =
00227 *(reinterpret_cast<const RangeDomainMeasurementVectorType* >
00228 (&(*this->GetPixelContainer())[id]));
00229 }
00230 else
00231 {
00232 m_TempRangeVector =
00233 *(reinterpret_cast< const RangeDomainMeasurementVectorType* >
00234 (&(this->GetImage()->GetPixel( m_TempIndex ) ) ) );
00235 }
00236
00237 for ( unsigned int i = TImage::ImageDimension; i < MeasurementVectorType::Length; ++i )
00238 {
00239 m_TempVector[i] = m_TempRangeVector[i - TImage::ImageDimension]
00240 / m_NormalizationFactors[i];
00241 }
00242
00243 return m_TempVector;
00244 }
00245
00246 template < class TImage >
00247 inline void
00248 JointDomainImageToListAdaptor< TImage >
00249 ::ComputeRegion(const MeasurementVectorType& mv,
00250 const double radius,
00251 ImageRegionType& region) const
00252 {
00253 ImageIndexType beginIndex;
00254 ImageSizeType size;
00255
00256 for ( unsigned int i = 0; i < TImage::ImageDimension; ++i )
00257 {
00258 m_TempPoint[i] = m_NormalizationFactors[i] * (mv[i] - radius);
00259 size[i] = (unsigned long)(2.0 * m_NormalizationFactors[i] * radius
00260 / this->GetImage()->GetSpacing()[i]);
00261 }
00262
00263 this->GetImage()->TransformPhysicalPointToIndex(m_TempPoint , beginIndex );
00264
00265 for ( unsigned int i = 0; i < TImage::ImageDimension; ++i )
00266 {
00267 if ( beginIndex[i] < this->GetImageBeginIndex()[i] )
00268 {
00269 beginIndex[i] = this->GetImageBeginIndex()[i];
00270 size[i] -= (this->GetImageBeginIndex()[i] - beginIndex[i]);
00271 }
00272
00273 if ( static_cast<typename ImageIndexType::IndexValueType>(beginIndex[i] + size[i] - 1) > this->GetImageEndIndex()[i] )
00274 {
00275 size[i] = this->GetImageEndIndex()[i] - beginIndex[i] + 1;
00276 }
00277 }
00278
00279 region.SetIndex( beginIndex );
00280 region.SetSize( size );
00281 }
00282
00283 template < class TImage >
00284 inline void
00285 JointDomainImageToListAdaptor< TImage >
00286 ::Search(const MeasurementVectorType& mv,
00287 const double radius,
00288 InstanceIdentifierVectorType& result) const
00289 {
00290 ImageRegionType region;
00291 this->ComputeRegion( mv, radius, region );
00292
00293 result.clear();
00294 ImageIteratorType iter( this->GetImage(), region );
00295 iter.GoToBegin();
00296 double squaredRadius = radius * radius;
00297 InstanceIdentifier instanceID;
00298 while ( !iter.IsAtEnd() )
00299 {
00300 instanceID = this->GetImage()->ComputeOffset(iter.GetIndex());
00301 m_TempVector = this->GetMeasurementVector( instanceID );
00302 bool isWithinRange = true;
00303 double sum = 0.0;
00304 for ( unsigned int i = 0; i < MeasurementVectorSize; ++i )
00305 {
00306 const double temp = (m_TempVector[i] - mv[i]);
00307 sum += temp * temp;
00308 if ( sum > squaredRadius )
00309 {
00310 isWithinRange = false;
00311 break;
00312 }
00313 }
00314
00315 if ( isWithinRange )
00316 {
00317 result.push_back(instanceID);
00318 }
00319 ++iter;
00320 }
00321 }
00322
00323 }
00324 }
00325
00326 #ifndef ITK_MANUAL_INSTANTIATION
00327 #include "itkJointDomainImageToListAdaptor.txx"
00328 #endif
00329
00330 #endif
00331