diff --git a/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp b/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp index de773d509f..0cd073efc7 100644 --- a/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp +++ b/Modules/MitkExt/Algorithms/mitkPlaneFit.cpp @@ -1,212 +1,229 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #include "mitkPlaneFit.h" #include "mitkPlaneGeometry.h" #include "mitkGeometryData.h" #include #include mitk::PlaneFit::PlaneFit() : m_PointSet( NULL ) { m_TimeSlicedGeometry = mitk::TimeSlicedGeometry::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->GetTimeSlicedGeometry() == NULL ) update = true; if ( ( ! update ) && ( output->GetTimeSlicedGeometry()->GetTimeSteps() != input->GetTimeSlicedGeometry()->GetTimeSteps() ) ) update = true; if ( update ) { mitk::PlaneGeometry::Pointer planeGeometry = mitk::PlaneGeometry::New(); m_TimeSlicedGeometry->InitializeEvenlyTimed( planeGeometry, m_PointSet->GetPointSetSeriesSize() ); unsigned int t; for ( t = 0; (t < m_PointSet->GetPointSetSeriesSize()) && (t < m_Planes.size()); ++t ) { m_TimeSlicedGeometry->SetGeometry3D( m_Planes[t], (int) t ); } output->SetGeometry( m_TimeSlicedGeometry ); } } 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; - // sum of all points - mitk::PointSet::PointsContainer::Iterator pit, end; - pit = m_PointSet->GetPointSet( t )->GetPoints()->Begin(); - end = m_PointSet->GetPointSet( t )->GetPoints()->End(); - for ( ; pit!=end; ++pit ) + //// sum of all points + //mitk::PointSet::PointsContainer::Iterator pit, end; + //pit = m_PointSet->GetPointSet( t )->GetPoints()->Begin(); + //end = m_PointSet->GetPointSet( t )->GetPoints()->End(); + //for ( ; pit!=end; ++pit ) + //{ + // mitk::Point3D p3d = pit.Value(); + // m_Centroids[t][0] += p3d[0]; + // m_Centroids[t][1] += p3d[1]; + // m_Centroids[t][2] += p3d[2]; + //} + + 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); - // calculate point distance to centroid and inserting it in the matrix - mitk::PointSet::PointsContainer::Iterator pit, end; - pit = m_PointSet->GetPointSet( t )->GetPoints()->Begin(); - end = m_PointSet->GetPointSet( t )->GetPoints()->End(); - - for (int p=0; pit!=end; pit++, p++) + //// calculate point distance to centroid and inserting it in the matrix + //mitk::PointSet::PointsContainer::Iterator pit, end; + //pit = m_PointSet->GetPointSet( t )->GetPoints()->Begin(); + //end = m_PointSet->GetPointSet( t )->GetPoints()->End(); + // + //for (int p=0; pit!=end; pit++, p++) + //{ + // mitk::Point3D p3d = pit.Value(); + // dataM[p][0] = p3d[0] - m_Centroids[t][0]; + // dataM[p][1] = p3d[1] - m_Centroids[t][1]; + // dataM[p][2] = p3d[2] - m_Centroids[t][2]; + //} + + + 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] ); }