<commit_before>#ifndef __elxEulerTransform_HXX_
#define __elxEulerTransform_HXX_

#include "elxEulerTransform.h"

namespace elastix
{
	using namespace itk;
	
	/**
	 * ********************* Constructor ****************************
	 */
	
	template <class TElastix>
		EulerTransformElastix<TElastix>
		::EulerTransformElastix()
	{
	} // end Constructor
	
	
	/**
	 * ******************* BeforeRegistration ***********************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::BeforeRegistration(void)
	{
		/** Task 1 - Set initial parameters.*/
		ParametersType dummyInitialParameters( this->GetNumberOfParameters() );
		dummyInitialParameters.Fill(0.0);

		/** And give it to this->m_Registration.*/
		this->m_Registration->GetAsITKBaseType()
			->SetInitialTransformParameters( dummyInitialParameters );
		
		/** Task 2 - Set the scales.*/
		/** Here is an heuristic rule for estimating good values for
		 * the rotation/translation scales.
		 *
		 * 1) Estimate the bounding box of your points (in physical units).
		 * 2) Take the 3D Diagonal of that bounding box
		 * 3) Multiply that by 10.0.
		 * 4) use 1.0 /[ value from (3) ] as the translation scaling value.
		 * 5) use 1.0 as the rotation scaling value.
		 *
		 * With this operation you bring the translation units
		 * to the range of rotations (e.g. around -1 to 1).
		 * After that, all your registration parameters are
		 * in the relaxed range of -1:1. At that point you
		 * can start setting your optimizer with step lengths
		 * in the ranges of 0.001 if you are conservative, or
		 * in the range of 0.1 if you want to live dangerously.
		 * (0.1 radians is about 5.7 degrees).
		 * 
		 * This heuristic rule is based on the naive assumption
		 * that your registration may require translations as
		 * large as 1/10 of the diagonal of the bounding box.
		 */
		ScalesType newscales( this->GetNumberOfParameters() );
		newscales.Fill(1.0);
		double scaler = 100000.0;
		this->m_Configuration->ReadParameter( scaler, "Scaler", 0 );
		/**
		 * - If the Dimension is 3, the first 3 parameters represent angles.
		 * - If the Dimension is 2, only the first parameter represent an angle.
		 */
		unsigned int AnglePartEnd = 3;
		if ( SpaceDimension == 2 ) AnglePartEnd = 1;
		for ( unsigned int i = 0; i < AnglePartEnd; i++ )
		{
			newscales[ i ] *= scaler;
		}

		/** And give it to the optimizer. */
		this->m_Registration->GetAsITKBaseType()->GetOptimizer()->SetScales( newscales );

		/** Task 3 - Set center of rotation. */
		InputPointType rotationPoint;
		this->CalculateRotationPoint( rotationPoint );
		this->Superclass1::SetCenter( rotationPoint );
		
	} // end BeforeRegistration
	
	
	
	
	/**
	 * ***************** CalculateRotationPoint *********************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::CalculateRotationPoint( InputPointType & rotationPoint )
	{
		/** Fill rotationPoint.
		 * The CenterOfRotation is set by default to the middle
		 * of the fixed image. If CenterOfRotation is specified
		 * in the Parameterfile, then the default value is overwritten.
		 */
		
		/** Get fixed Image size.*/
		SizeType fixedSize = this->m_Registration->GetAsITKBaseType()->
			GetFixedImage()->GetLargestPossibleRegion().GetSize();
		
		/** Fill CenterOfRotationIndices,
		 * which is the rotationPoint, expressed in index-values.
		 */
		IndexType CenterOfRotationIndices;
		bool CORInImage = true;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			CenterOfRotationIndices[ i ] = static_cast<IndexValueType>( fixedSize[ i ] / 2 );
			this->m_Configuration->ReadParameter( CenterOfRotationIndices[ i ], "CenterOfRotation", i, true );
			/** Check if CenterOfRotation has index-values within image.*/
			if ( CenterOfRotationIndices[ i ] < 0 || CenterOfRotationIndices[ i ] > fixedSize[ i ] )
			{
				CORInImage = false;
			}
		}
		
		/** Give a warning if necessary.*/
		if ( !CORInImage )
		{
			xl::xout["warning"] << "WARNING: Center of Rotation is not within image boundaries!" << std::endl;
		}
		
		/** Convert from index-value to physical-point-value.*/
		this->m_Registration->GetAsITKBaseType()->GetFixedImage()->
			TransformIndexToPhysicalPoint( CenterOfRotationIndices, rotationPoint );

	} // end CalculateRotationPoint
	
	
	/**
	 * ************************* ReadFromFile ************************
	 */

	template <class TElastix>
	void EulerTransformElastix<TElastix>::
		ReadFromFile(void)
	{
		/** Read the center of rotation.*/
		IndexType rotationIndex;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			rotationIndex[ i ] = 0;
			this->m_Configuration->ReadParameter( rotationIndex[ i ], "CenterOfRotation", i, true );		
		}
		
		/** Get spacing, origin and size of the fixed image.
		 * We put this in a dummy image, so that we can correctly
		 * calculate the center of rotation in world coordinates.
		 */
		SpacingType		spacing;
		IndexType			index;
		PointType			origin;
		SizeType			size;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			/** No default size. Read size from the parameter file. */
			this->m_Configuration->ReadParameter(	size[ i ], "Size", i );

			/** Default index. Read index from the parameter file. */
			index[ i ] = 0;
			this->m_Configuration->ReadParameter(	index[ i ], "Index", i );

			/** Default spacing. Read spacing from the parameter file. */
			spacing[ i ] = 1.0;
			this->m_Configuration->ReadParameter(	spacing[ i ], "Spacing", i );

			/** Default origin. Read origin from the parameter file. */
			origin[ i ] = 0.0;
			this->m_Configuration->ReadParameter(	origin[ i ], "Origin", i );
		}

		/** Check for image size. */
		unsigned int sum = 0;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			if ( size[ i ] == 0 ) sum++;
		}
		if ( sum > 0 )
		{
			xl::xout["error"] << "ERROR: One or more image sizes are 0!" << std::endl;
			/** \todo quit program nicely. */
		}
		
		/** Make a temporary image with the right region info,
		* so that the TransformIndexToPhysicalPoint-functions will be right.
		*/
		typename DummyImageType::Pointer dummyImage = DummyImageType::New();
		RegionType region;
		region.SetIndex( index );
		region.SetSize( size );
		dummyImage->SetRegions( region );
		dummyImage->SetOrigin( origin );
		dummyImage->SetSpacing( spacing );

		/** region.SetIndex( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputStartIndex() );
		region.SetSize( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetSize() );
		dummyImage->SetRegions( region );
		dummyImage->SetOrigin( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputOrigin() );
		dummyImage->SetSpacing( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputSpacing() );*/
				
		/** Convert center of rotation from index-value to physical-point-value.*/
		InputPointType rotationPoint;
		dummyImage->TransformIndexToPhysicalPoint( rotationIndex, rotationPoint );

		/** Set it in this Transform.*/
		this->SetCenter( rotationPoint );

		/** Call the ReadFromFile from the TransformBase.
		 * BE AWARE: Only call Superclass2::ReadFromFile() after CenterOfRotation
		 * is set, because it is used in the SetParameters()-function of this transform.
		 */
		this->Superclass2::ReadFromFile();

	} // end ReadFromFile


	/**
	 * ************************* WriteToFile ************************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::WriteToFile( const ParametersType & param )
	{
		/** Call the WriteToFile from the TransformBase.*/
		this->Superclass2::WriteToFile( param );

		/** Write AffineTransform specific things.*/
		xout["transpar"] << std::endl << "// EulerTransform specific" << std::endl;

		/** Get the center of rotation and convert it from 
		 * physical-point-value to index-value.
		 */
		IndexType rotationIndex;
		InputPointType rotationPoint = this->GetCenter();
		this->m_Registration->GetAsITKBaseType()->GetFixedImage()->
			TransformPhysicalPointToIndex( rotationPoint, rotationIndex );

		/** Write the center of rotation.*/
		xout["transpar"] << "(CenterOfRotation ";
		for ( unsigned int i = 0; i < SpaceDimension - 1; i++ )
		{
			xout["transpar"] << rotationIndex[ i ] << " ";
		}
		xout["transpar"] << rotationIndex[ SpaceDimension - 1 ] << ")" << std::endl;
	
	} // end WriteToFile
	

} // end namespace elastix


#endif // end #ifndef __elxEulerTransform_HXX_

<commit_msg>SKlinux:<commit_after>#ifndef __elxEulerTransform_HXX_
#define __elxEulerTransform_HXX_

#include "elxEulerTransform.h"

namespace elastix
{
	using namespace itk;
	
	/**
	 * ********************* Constructor ****************************
	 */
	
	template <class TElastix>
		EulerTransformElastix<TElastix>
		::EulerTransformElastix()
	{
	} // end Constructor
	
	
	/**
	 * ******************* BeforeRegistration ***********************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::BeforeRegistration(void)
	{
		/** Task 1 - Set initial parameters.*/
		ParametersType dummyInitialParameters( this->GetNumberOfParameters() );
		dummyInitialParameters.Fill(0.0);

		/** And give it to this->m_Registration.*/
		this->m_Registration->GetAsITKBaseType()
			->SetInitialTransformParameters( dummyInitialParameters );
		
		/** Task 2 - Set the scales.*/
		/** Here is an heuristic rule for estimating good values for
		 * the rotation/translation scales.
		 *
		 * 1) Estimate the bounding box of your points (in physical units).
		 * 2) Take the 3D Diagonal of that bounding box
		 * 3) Multiply that by 10.0.
		 * 4) use 1.0 /[ value from (3) ] as the translation scaling value.
		 * 5) use 1.0 as the rotation scaling value.
		 *
		 * With this operation you bring the translation units
		 * to the range of rotations (e.g. around -1 to 1).
		 * After that, all your registration parameters are
		 * in the relaxed range of -1:1. At that point you
		 * can start setting your optimizer with step lengths
		 * in the ranges of 0.001 if you are conservative, or
		 * in the range of 0.1 if you want to live dangerously.
		 * (0.1 radians is about 5.7 degrees).
		 * 
		 * This heuristic rule is based on the naive assumption
		 * that your registration may require translations as
		 * large as 1/10 of the diagonal of the bounding box.
		 */
		ScalesType newscales( this->GetNumberOfParameters() );
		newscales.Fill(1.0);
		double scaler = 100000.0;
		this->m_Configuration->ReadParameter( scaler, "Scaler", 0 );
		/**
		 * - If the Dimension is 3, the first 3 parameters represent angles.
		 * - If the Dimension is 2, only the first parameter represent an angle.
		 */
		unsigned int AnglePartEnd = 3;
		if ( SpaceDimension == 2 ) AnglePartEnd = 1;
		for ( unsigned int i = 0; i < AnglePartEnd; i++ )
		{
			newscales[ i ] *= scaler;
		}

		/** And give it to the optimizer. */
		this->m_Registration->GetAsITKBaseType()->GetOptimizer()->SetScales( newscales );

		/** Task 3 - Set center of rotation. */
		InputPointType rotationPoint;
		this->CalculateRotationPoint( rotationPoint );
		this->Superclass1::SetCenter( rotationPoint );
		
	} // end BeforeRegistration
	
	
	
	
	/**
	 * ***************** CalculateRotationPoint *********************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::CalculateRotationPoint( InputPointType & rotationPoint )
	{
		/** Fill rotationPoint.
		 * The CenterOfRotation is set by default to the middle
		 * of the fixed image. If CenterOfRotation is specified
		 * in the Parameterfile, then the default value is overwritten.
		 */
		
		/** Get fixed Image size.*/
		SizeType fixedSize = this->m_Registration->GetAsITKBaseType()->
			GetFixedImage()->GetLargestPossibleRegion().GetSize();
		
		/** Fill CenterOfRotationIndices,
		 * which is the rotationPoint, expressed in index-values.
		 */
		IndexType CenterOfRotationIndices;
		bool CORInImage = true;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			CenterOfRotationIndices[ i ] = static_cast<IndexValueType>( fixedSize[ i ] / 2 );
			this->m_Configuration->ReadParameter( CenterOfRotationIndices[ i ], "CenterOfRotation", i, true );
			/** Check if CenterOfRotation has index-values within image.*/
			if ( CenterOfRotationIndices[ i ] < 0 || CenterOfRotationIndices[ i ] > fixedSize[ i ] )
			{
				CORInImage = false;
			}
		}
		
		/** Give a warning if necessary.*/
		if ( !CORInImage )
		{
			xl::xout["warning"] << "WARNING: Center of Rotation is not within image boundaries!" << std::endl;
		}
		
		/** Convert from index-value to physical-point-value.*/
		this->m_Registration->GetAsITKBaseType()->GetFixedImage()->
			TransformIndexToPhysicalPoint( CenterOfRotationIndices, rotationPoint );

	} // end CalculateRotationPoint
	
	
	/**
	 * ************************* ReadFromFile ************************
	 */

	template <class TElastix>
	void EulerTransformElastix<TElastix>::
		ReadFromFile(void)
	{
		/** Read the center of rotation.*/
		IndexType rotationIndex;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			rotationIndex[ i ] = 0;
			this->m_Configuration->ReadParameter( rotationIndex[ i ], "CenterOfRotation", i, true );		
		}
		
		/** Get spacing, origin and size of the fixed image.
		 * We put this in a dummy image, so that we can correctly
		 * calculate the center of rotation in world coordinates.
		 */
		SpacingType		spacing;
		IndexType			index;
		PointType			origin;
		SizeType			size;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			/** No default size. Read size from the parameter file. */
			this->m_Configuration->ReadParameter(	size[ i ], "Size", i );

			/** Default index. Read index from the parameter file. */
			index[ i ] = 0;
			this->m_Configuration->ReadParameter(	index[ i ], "Index", i );

			/** Default spacing. Read spacing from the parameter file. */
			spacing[ i ] = 1.0;
			this->m_Configuration->ReadParameter(	spacing[ i ], "Spacing", i );

			/** Default origin. Read origin from the parameter file. */
			origin[ i ] = 0.0;
			this->m_Configuration->ReadParameter(	origin[ i ], "Origin", i );
		}

		/** Check for image size. */
		unsigned int sum = 0;
		for ( unsigned int i = 0; i < SpaceDimension; i++ )
		{
			if ( size[ i ] == 0 ) sum++;
		}
		if ( sum > 0 )
		{
			xl::xout["error"] << "ERROR: One or more image sizes are 0!" << std::endl;
			/** \todo quit program nicely. */
		}
		
		/** Make a temporary image with the right region info,
		* so that the TransformIndexToPhysicalPoint-functions will be right.
		*/
		typedef FixedImageType DummyImageType; 
		typename DummyImageType::Pointer dummyImage = DummyImageType::New();
		RegionType region;
		region.SetIndex( index );
		region.SetSize( size );
		dummyImage->SetRegions( region );
		dummyImage->SetOrigin( origin );
		dummyImage->SetSpacing( spacing );

		/** region.SetIndex( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputStartIndex() );
		region.SetSize( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetSize() );
		dummyImage->SetRegions( region );
		dummyImage->SetOrigin( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputOrigin() );
		dummyImage->SetSpacing( (dynamic_cast< typename ElastixType::ResamplerBaseType *>( this->m_Elastix->GetResampler() ))
			->GetAsITKBaseType()->GetOutputSpacing() );*/
				
		/** Convert center of rotation from index-value to physical-point-value.*/
		InputPointType rotationPoint;
		dummyImage->TransformIndexToPhysicalPoint( rotationIndex, rotationPoint );

		/** Set it in this Transform.*/
		this->SetCenter( rotationPoint );

		/** Call the ReadFromFile from the TransformBase.
		 * BE AWARE: Only call Superclass2::ReadFromFile() after CenterOfRotation
		 * is set, because it is used in the SetParameters()-function of this transform.
		 */
		this->Superclass2::ReadFromFile();

	} // end ReadFromFile


	/**
	 * ************************* WriteToFile ************************
	 */
	
	template <class TElastix>
		void EulerTransformElastix<TElastix>
		::WriteToFile( const ParametersType & param )
	{
		/** Call the WriteToFile from the TransformBase.*/
		this->Superclass2::WriteToFile( param );

		/** Write AffineTransform specific things.*/
		xout["transpar"] << std::endl << "// EulerTransform specific" << std::endl;

		/** Get the center of rotation and convert it from 
		 * physical-point-value to index-value.
		 */
		IndexType rotationIndex;
		InputPointType rotationPoint = this->GetCenter();
		this->m_Registration->GetAsITKBaseType()->GetFixedImage()->
			TransformPhysicalPointToIndex( rotationPoint, rotationIndex );

		/** Write the center of rotation.*/
		xout["transpar"] << "(CenterOfRotation ";
		for ( unsigned int i = 0; i < SpaceDimension - 1; i++ )
		{
			xout["transpar"] << rotationIndex[ i ] << " ";
		}
		xout["transpar"] << rotationIndex[ SpaceDimension - 1 ] << ")" << std::endl;
	
	} // end WriteToFile
	

} // end namespace elastix


#endif // end #ifndef __elxEulerTransform_HXX_


{'fraction_non_alphanumeric': 0.14592571264036855, 'fraction_numerical': 0.008753239274402533, 'mean_word_length': 3.6695348211884915, 'pattern_counts': {'":': 0, '<': 99, '<?xml version=': 0, '>': 167, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1217229', 'n_tokens_mistral': 5933, 'n_tokens_neox': 5070, 'n_words': 1533}