[Insight-users] Versor rotation
Arunachalam Kana
Kana.Arunachalam at fh-wels.at
Thu Sep 2 07:16:24 EDT 2010
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);
typedef itk::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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100902/36546b95/attachment-0001.htm>
More information about the Insight-users
mailing list