diff --git a/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp b/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp index 61e03d1982..ad2239fabf 100644 --- a/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp +++ b/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp @@ -1,202 +1,205 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPlaneFit.h" #include "mitkPlaneGeometry.h" #include "mitkGeometryData.h" +#include #include #include mitk::PlaneFit::PlaneFit() : m_PointSet( NULL ) { - m_TimeSlicedGeometry = mitk::TimeSlicedGeometry::New(); + m_TimeGeometry = mitk::ProportionalTimeGeometry::New(); } mitk::PlaneFit::~PlaneFit() { } void mitk::PlaneFit::GenerateOutputInformation() { mitk::PointSet::ConstPointer input = this->GetInput(); mitk::GeometryData::Pointer output = this->GetOutput(); itkDebugMacro(<<"GenerateOutputInformation()"); if (input.IsNull()) return; if ( m_PointSet == NULL ) { return; } bool update = false; if ( output->GetGeometry() == NULL || output->GetTimeGeometry() == NULL ) update = true; if ( ( ! update ) && ( output->GetTimeGeometry()->GetNumberOfTimeSteps() != input->GetTimeGeometry()->GetNumberOfTimeSteps() ) ) update = true; if ( update ) { mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); - m_TimeSlicedGeometry->InitializeEvenlyTimed( - planeGeometry, m_PointSet->GetPointSetSeriesSize() ); + ProportionalTimeGeometry::Pointer timeGeometry = dynamic_cast(m_TimeGeometry.GetPointer()); + timeGeometry->Initialize(planeGeometry, m_PointSet->GetPointSetSeriesSize()); + //m_TimeGeometry->InitializeEvenlyTimed( + // planeGeometry, m_PointSet->GetPointSetSeriesSize() ); - unsigned int t; - for ( t = 0; - (t < m_PointSet->GetPointSetSeriesSize()) - && (t < m_Planes.size()); - ++t ) + TimeStepType timeStep; + for ( timeStep = 0; + (timeStep < m_PointSet->GetPointSetSeriesSize()) + && (timeStep < m_Planes.size()); + ++timeStep ) { - m_TimeSlicedGeometry->SetGeometry3D( m_Planes[t], (int) t ); + timeGeometry->SetTimeStepGeometry( m_Planes[timeStep], timeStep ); } - output->SetGeometry( m_TimeSlicedGeometry ); + output->SetTimeGeometry( m_TimeGeometry ); } } void mitk::PlaneFit::GenerateData() { unsigned int t; for ( t = 0; t < m_PointSet->GetPointSetSeriesSize(); ++t ) { // check number of data points - less then 3points isn't enough if ( m_PointSet->GetSize( t ) >= 3 ) { this->CalculateCentroid( t ); this->ProcessPointSet( t ); this->InitializePlane( t ); } } } void mitk::PlaneFit::SetInput( const mitk::PointSet* pointSet ) { // Process object is not const-correct so the const_cast is required here this->ProcessObject::SetNthInput(0, const_cast< mitk::PointSet * >( pointSet ) ); m_PointSet = pointSet; unsigned int pointSetSize = pointSet->GetPointSetSeriesSize(); m_Planes.resize( pointSetSize ); m_Centroids.resize( pointSetSize ); m_PlaneVectors.resize( pointSetSize ); unsigned int t; for ( t = 0; t < pointSetSize; ++t ) { m_Planes[t] = mitk::PlaneGeometry::New(); } } const mitk::PointSet* mitk::PlaneFit::GetInput() { if (this->GetNumberOfInputs() < 1) { return 0; } return static_cast (this->ProcessObject::GetInput(0) ); } void mitk::PlaneFit::CalculateCentroid( int t ) { if ( m_PointSet == NULL ) return; int ps_total = m_PointSet->GetSize( t ); m_Centroids[t][0] = m_Centroids[t][1] = m_Centroids[t][2] = 0.0; for (int i=0; iGetPoint(i,t); m_Centroids[t][0] += p3d[0]; m_Centroids[t][1] += p3d[1]; m_Centroids[t][2] += p3d[2]; } // calculation of centroid m_Centroids[t][0] /= ps_total; m_Centroids[t][1] /= ps_total; m_Centroids[t][2] /= ps_total; } void mitk::PlaneFit::ProcessPointSet( int t ) { if (m_PointSet == NULL ) return; // int matrix with POINTS x (X,Y,Z) vnl_matrix dataM( m_PointSet->GetSize( t ), 3); int ps_total = m_PointSet->GetSize( t ); for (int i=0; iGetPoint(i,t); dataM[i][0] = p3d[0] - m_Centroids[t][0]; dataM[i][1] = p3d[1] - m_Centroids[t][1]; dataM[i][2] = p3d[2] - m_Centroids[t][2]; } // process the SVD (singular value decomposition) from ITK // the vector will be orderd descending vnl_svd svd(dataM, 0.0); // calculate the SVD of A vnl_vector v = svd.nullvector(); // Avoid erratic normal sign switching when the plane changes minimally // by negating the vector for negative x values. if ( v[0] < 0 ) { v = -v; } m_PlaneVectors[t][0] = v[0]; m_PlaneVectors[t][1] = v[1]; m_PlaneVectors[t][2] = v[2]; } mitk::PlaneGeometry::Pointer mitk::PlaneFit::GetPlaneGeometry( int t ) { return m_Planes[t]; } const mitk::Vector3D &mitk::PlaneFit::GetPlaneNormal( int t ) const { return m_PlaneVectors[t]; } const mitk::Point3D &mitk::PlaneFit::GetCentroid( int t ) const { return m_Centroids[t]; } void mitk::PlaneFit::InitializePlane( int t ) { m_Planes[t]->InitializePlane( m_Centroids[t], m_PlaneVectors[t] ); } diff --git a/Modules/MitkExt/Algorithms/mitkPlaneFit.h b/Modules/MitkExt/Algorithms/mitkPlaneFit.h index 976dc4daec..74fd2ec6ad 100644 --- a/Modules/MitkExt/Algorithms/mitkPlaneFit.h +++ b/Modules/MitkExt/Algorithms/mitkPlaneFit.h @@ -1,147 +1,147 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #if !defined(MITK_PLANEFIT_H__INCLUDED_) #define MITK_PLANEFIT_H__INCLUDED_ #include "mitkPointSet.h" #include "MitkExtExports.h" -#include "mitkTimeSlicedGeometry.h" +#include "mitkTimeGeometry.h" #include "mitkPlaneGeometry.h" #include "mitkGeometryDataSource.h" namespace mitk { //! // kind regards to dr. math! // function [x0, a, d, normd] = lsplane(X) // --------------------------------------------------------------------- // LSPLANE.M Least-squares plane (orthogonal distance // regression). // // Version 1.0 // Last amended I M Smith 27 May 2002. // Created I M Smith 08 Mar 2002 // --------------------------------------------------------------------- // Input // X Array [x y z] where x = vector of x-coordinates, // y = vector of y-coordinates and z = vector of // z-coordinates. // Dimension: m x 3. // // Output // x0 Centroid of the data = point on the best-fit plane. // Dimension: 3 x 1. // // a Direction cosines of the normal to the best-fit // plane. // Dimension: 3 x 1. // // // // [x0, a <, d, normd >] = lsplane(X) // --------------------------------------------------------------------- class MitkExt_EXPORT PlaneFit : public GeometryDataSource { public: mitkClassMacro( PlaneFit, GeometryDataSource); itkNewMacro(Self); typedef mitk::PointSet::PointDataType PointDataType; typedef mitk::PointSet::PointDataIterator PointDataIterator; virtual void GenerateOutputInformation(); virtual void GenerateData(); /*!Getter for point set. * */ const mitk::PointSet *GetInput(); /*! filter initialisation. * */ virtual void SetInput( const mitk::PointSet *ps ); /*! returns the center of gravity of the point set. * */ virtual const mitk::Point3D &GetCentroid( int t = 0 ) const; /*! returns the plane geometry which represents the point set. * */ virtual mitk::PlaneGeometry::Pointer GetPlaneGeometry( int t = 0 ); /*! returns the normal of the plane which represents the point set. * */ virtual const mitk::Vector3D &GetPlaneNormal( int t = 0 ) const; protected: PlaneFit(); virtual ~PlaneFit(); /*! Calculates the centroid of the point set. * the center of gravity is calculated through the mean value of the whole point set */ void CalculateCentroid( int t = 0 ); /*! working with an SVD algorithm form matrix dataM. * ITK suplies the vnl_svd to solve an plan fit eigentvector problem * points are processed in the SVD matrix. The normal vector is the * singular vector of dataM corresponding to its smalest singular value. * The mehtod uses VNL library from ITK and at least the mehtod nullvector() * to extract the normalvector. */ void ProcessPointSet( int t = 0 ); /*! Initialize Plane and configuration. * */ void InitializePlane( int t = 0 ); private: /*!keeps a copy of the pointset.*/ const mitk::PointSet* m_PointSet; /* output object - a time sliced geometry.*/ - mitk::TimeSlicedGeometry::Pointer m_TimeSlicedGeometry; + mitk::TimeGeometry::Pointer m_TimeGeometry; std::vector< mitk::PlaneGeometry::Pointer > m_Planes; /*! the calculatet center point of all points in the point set.*/ std::vector< mitk::Point3D > m_Centroids; /* the normal vector to descrie a plane gemoetry.*/ std::vector< mitk::Vector3D > m_PlaneVectors; }; }//namespace mitk #endif //MITK_PLANFIT_INCLUDE_