[Insight-users] Versor rotation
Luis Ibanez
luis.ibanez at kitware.com
Fri Sep 3 11:43:23 EDT 2010
Hi Arunachalam,
The Matrix of EigenVectors is exactly the rotation matrix
that you need to apply in order to align with the main
directions of the Hessian.
You can use that matrix to initialize a Versor, by taking
advantage of the versor method "Set( Matrix & )", as
described in:
http://www.itk.org/Doxygen312/html/classitk_1_1Versor.html#d22261bee493be34f2fc3941d24604dc
Regards,
Luis
-----------------------------------------------------------------------
On Thu, Sep 2, 2010 at 7:16 AM, Arunachalam Kana <
Kana.Arunachalam at fh-wels.at> wrote:
> Hi itkusers,
>
>
>
> I have cylinder which can be oriented in any angle in 3D space. I calculate
> the direction of cylinder using hessian and eigen value analysis.
>
> My goal is to draw the intensity profile of the plane perpendicular to the
> cylinder direction. The size of the plane is 30 x 30 with the current
>
> Point (point for which the hessian is calculate) as center.
>
>
>
> My pipeline is as follows: Read image -> cast image to float -> calculate
> second derivatives ->create hessian matrix ->calculate eigen values and
> vector.
>
>
>
> So now i have one vector showing the direction of cylinder (called normal)
> and the other two orthogonal to it (called ortho_1 and ortho_2). Till this
> point
>
> I am clear.
>
> As my goal is to extract the plane which is formed by ortho_1 and ortho_2
> and of size 30x30, i decided to do the following:
>
> 1. to align the plane create by ortho_1 and ortho_2 parallel to z
> axis - rotate the image using versor transform and align to z direction. I
> referred to versorrigidtransform example in itkguide.
>
> 2. Translate the point to be in first slice.
>
> 3. Select a region in the first slice and write the gray values as
> csv.
>
>
>
> I am stuck in rotation step. My rotation is not performing as expected. I
> think something is missing, don’t know where or what ? I have pasted my
> code below.
>
> If anyone thinks that this idea is a over kill for my goal, other ideas is
> also welcomed.
>
>
>
> int main(int argc, char* argv [] )
>
> {
>
> //image input filenames
>
> const char * inputFilename =
> "C:/LIBS/SW_Development/GradientAnalysis/testdata/itkVirtual-SingleFibresInXY-r3.mhd"
> ;
>
>
>
> //declare pixeltype and dimension
>
> typedef double InputPixelType;
>
> typedef float CastPixelType;
>
> typedef int CCPixelType;
>
>
>
> typedef itk::SymmetricSecondRankTensor< float, Dimension >
> T_PixelType;
>
> typedef itk::Vector< VectorPixelType, Dimension > EV_PixelType;
>
>
>
> //declare image type
>
> typedef itk::Image< InputPixelType, Dimension > InputImageType;
>
> typedef itk::Image< CastPixelType, Dimension > CastImageType;
>
> typedef itk::Image< CCPixelType, Dimension > CCImageType;
>
> typedef CCImageType::IndexType CCIndexType;
>
>
>
> //set image reader filter parameters
>
> typedef itk::ImageFileReader< InputImageType > ReaderType;
>
> ReaderType::Pointer reader_1 = ReaderType::New();
>
> reader_1->SetFileName( inputFilename );
>
> reader_1->Update();
>
>
>
> //set cast filter parameters to int as the derivative can be
> positive and negative
>
> typedef itk::CastImageFilter< InputImageType, CastImageType >
> CastType;
>
> CastType::Pointer castfilter = CastType::New();
>
> castfilter->SetInput( reader_1->GetOutput() );
>
> castfilter->Update();
>
>
>
> //create derivative
>
> typedef itk::DerivativeImageFilter<CastImageType,CastImageType>
> DerivativeFilterType;
>
>
>
> //create First derivative in X direction
>
> DerivativeFilterType::Pointer DXfilter = DerivativeFilterType::New();
>
> DXfilter->SetInput(castfilter->GetOutput());
>
> DXfilter->SetDirection(0);
>
> DXfilter->SetOrder(1);
>
> DXfilter->Update();
>
>
>
> //create second derivative in XX direction
>
> DerivativeFilterType::Pointer DXXfilter =
> DerivativeFilterType::New();
>
> DXXfilter->SetInput(DXfilter->GetOutput());
>
> DXXfilter->SetDirection(0);
>
> DXXfilter->SetOrder(1);
>
> DXXfilter->Update();
>
>
>
> //create second derivative in XY direction
>
> DerivativeFilterType::Pointer DXYfilter =
> DerivativeFilterType::New();
>
> DXYfilter->SetInput(DXfilter->GetOutput());
>
> DXYfilter->SetDirection(1);
>
> DXYfilter->SetOrder(1);
>
> DXYfilter->Update();
>
>
>
> //create second derivative in XZ direction
>
> DerivativeFilterType::Pointer DXZfilter =
> DerivativeFilterType::New();
>
> DXZfilter->SetInput(DXfilter->GetOutput());
>
> DXZfilter->SetDirection(2);
>
> DXZfilter->SetOrder(1);
>
> DXZfilter->Update();
>
>
>
> //create first derivative in Y direction
>
> DerivativeFilterType::Pointer DYfilter = DerivativeFilterType::New();
>
> DYfilter->SetInput(castfilter->GetOutput());
>
> DYfilter->SetDirection(1);
>
> DYfilter->SetOrder(1);
>
> DYfilter->Update();
>
>
>
> //create second derivative in YY direction
>
> DerivativeFilterType::Pointer DYYfilter =
> DerivativeFilterType::New();
>
> DYYfilter->SetInput(DYfilter->GetOutput());
>
> DYYfilter->SetDirection(1);
>
> DYYfilter->SetOrder(1);
>
> DYYfilter->Update();
>
>
>
> //create second derivative in YZ direction
>
> DerivativeFilterType::Pointer DYZfilter =
> DerivativeFilterType::New();
>
> DYZfilter->SetInput(DYfilter->GetOutput());
>
> DYZfilter->SetDirection(2);
>
> DYZfilter->SetOrder(1);
>
> DYZfilter->Update();
>
>
>
> //create first derivative in Z direction
>
> DerivativeFilterType::Pointer DZfilter = DerivativeFilterType::New();
>
> DZfilter->SetInput(castfilter->GetOutput());
>
> DZfilter->SetDirection(2);
>
> DZfilter->SetOrder(1);
>
> DZfilter->Update();
>
>
>
> //create second derivative in ZZ direction
>
> DerivativeFilterType::Pointer DZZfilter =
> DerivativeFilterType::New();
>
> DZZfilter->SetInput(DZfilter->GetOutput());
>
> DZZfilter->SetDirection(2);
>
> DZZfilter->SetOrder(1);
>
> DZZfilter->Update();
>
>
>
> //set eigen analysis
>
> typedef itk::SymmetricEigenAnalysis < T_PixelType, VectorPixelType,
> EV_PixelType > EigAnalysisType;
>
> EigAnalysisType eig;
>
> eig.SetDimension( Dimension );
>
> eig.SetOrderEigenMagnitudes( true );
>
>
>
> //derivative image iterator
>
> typedef itk::ImageRegionIterator <CastImageType>
> CastImageIteratorType;
>
> CastImageIteratorType DXXIter(DXXfilter->GetOutput(),
> DXXfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType DXYIter(DXYfilter->GetOutput(),
> DXYfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType DXZIter(DXZfilter->GetOutput(),
> DXZfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType DYYIter(DYYfilter->GetOutput(),
> DYYfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType DYZIter(DYZfilter->GetOutput(),
> DYZfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType DZZIter(DZZfilter->GetOutput(),
> DZZfilter->GetOutput()->GetLargestPossibleRegion());
>
> CastImageIteratorType InIter(castfilter->GetOutput(),
> castfilter->GetOutput()->GetLargestPossibleRegion());
>
>
>
> std::multimap<CCPixelType,CCIndexType>::iterator it;
>
> typedef itk::ResampleImageFilter<CastImageType,CastImageType>
> ResampleFilterType;
>
> ResampleFilterType::Pointer rfilter = ResampleFilterType::New();
>
>
>
> InIter.GoToBegin();
>
> While ( !InIter.IsAtEnd() )
>
> {
>
> IndexType currentindex = InIter.GetIndex();
>
>
>
> //initialise the variables
>
> T_PixelType hess_matrix;
>
> VectorPixelType eigen_values;
>
> EV_PixelType eigen_vectors;
>
>
>
> //set index for various iterators
>
> DXXIter.SetIndex(currentindex);
>
> DXYIter.SetIndex(currentindex);
>
> DXZIter.SetIndex(currentindex);
>
> DYYIter.SetIndex(currentindex);
>
> DYZIter.SetIndex(currentindex);
>
> DZZIter.SetIndex(currentindex);
>
> CCIter.SetIndex(currentindex);
>
> InIter.SetIndex(currentindex);
>
>
>
> //upload the hessian matrix from different second
> derivative
>
> hess_matrix[0] = DXXIter.Get();
>
> hess_matrix[1] = DXYIter.Get();
>
> hess_matrix[2] = DXZIter.Get();
>
> hess_matrix[3] = DYYIter.Get();
>
> hess_matrix[4] = DYZIter.Get();
>
> hess_matrix[5] = DZZIter.Get();
>
>
>
> //calculate the eigen values of the hessian matrix
>
> eig.ComputeEigenValuesAndVectors( hess_matrix,
> eigen_values, eigen_vectors );
>
>
>
> ///get the eigen vectors
>
> VectorPixelType normal, ortho_1, ortho_2;
>
> normal = eigen_vectors[0]; //direction of the fiber
>
> ortho_1 = eigen_vectors[1]; //one the orthogonal axis
>
> ortho_2 = eigen_vectors[2]; //other orthogonal axis
>
>
>
> double AdotB = normal[2] ;
>
> double ModA_ModB = sqrt(pow(normal[0],2) +
> pow(normal[1],2) + pow(normal[2],2));
>
> double theta = acos( AdotB / ModA_ModB ) ;
>
> double angle = sin ( theta / 2 );
>
>
>
> typedef itk::VersorTransform <double> TransformType;
>
> TransformType::Pointer transform = TransformType::New();
>
>
>
> typedef TransformType::VersorType VersorType;
>
> typedef TransformType::AxisType AxisType;
>
> typedef TransformType::CenterType CenterType;
>
>
>
> AxisType axis;
>
> axis[0]= normal[0];
>
> axis[1]= normal[1];
>
> axis[2]= normal[2];
>
>
>
> VersorType rotation;
>
> rotation.Set(axis, angle);
>
> transform->SetRotation(rotation);
>
>
>
> CenterType center;
>
> center[0] = normal[0];
>
> center[1] = normal[1];
>
> center[2] = normal[2];
>
> transform->SetCenter(center);
>
>
>
> typedefitk::LinearInterpolateImageFunction<CastImageType,
> double > InterpolatorType;
>
> InterpolatorType::Pointer interpolator =
> InterpolatorType::New();
>
>
>
> CastImageType::SizeType imgSize;
>
> imgSize =
> ccfilter->GetOutput()->GetBufferedRegion().GetSize();
>
>
>
> rfilter->SetTransform( transform );
>
> rfilter->SetInterpolator( interpolator );
>
> rfilter->SetSize( imgSize );
>
> rfilter->SetInput( castfilter->GetOutput() );
>
> rfilter->Update();
>
>
>
> }//if
>
>
>
> }
>
>
>
> //set image writer parameters
>
> const char * outputFilename =
> "C:/LIBS/SW_Development/GradientAnalysis/testdata/itkVirtual-SingleFibresInXY-r3-MAE-resampletest.mhd"
> ;
>
> typedef itk::ImageFileWriter< CastImageType > WriterType;
>
> WriterType::Pointer writer = WriterType::New();
>
> writer->SetFileName( outputFilename );
>
> writer->SetInput(rfilter->GetOutput());
>
>
>
> try
>
> {
>
> writer->Update();
>
> }
>
> catch( itk::ExceptionObject & err )
>
> {
>
> std::cerr << "ExceptionObject caught !" << std::endl;
>
> std::cerr << err << std::endl;
>
> return EXIT_FAILURE;
>
> }
>
>
>
> return EXIT_SUCCESS;
>
> }
>
>
>
> Any help will be appreciated.
>
>
>
> Thank you,
>
> Regards,
>
> Kana Arunachalam Kannappan
>
> Research Associate
>
> FH OÖ Forschungs & Entwicklungs GmbH
>
> Stelzhamer Strasse 23,
>
> 4600 Wels,
>
> Austria.
>
> Phone: +43 (0)7242 72811 -4420
>
> kana.arunachalam at fh-wels.at
>
> www.fh-ooe.at; www.3dct.at
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100903/13eb22a7/attachment-0001.htm>
More information about the Insight-users
mailing list