Commit b3fd8495 authored by hjohnson's avatar hjohnson
Browse files

ENH:

In response to numerous constructive critisms of the newly added suite of BRAINSTools, a code review was performed in order to make the applications work consistenently and with similar user interfaces.

Summary of ChangeLog:
-Instrumented xml files so that moving images are properly associated with their transforms in the Slicer3 MRML Interface.
-Made command line arguments consistent across tools for specifying images, transforms, similar parameters.
-Hid many advanced command line arguments from standard user interface when used from Slicer3
-Improved documentation of command line arguments to describe better what the intended purpose, and restrictions on use are.
-Added common set of image resample/warping options across all tools so that all tools provide common interface for choosing outputPixelType and interpolation mode
-Fixed reading and writing of transforms to use a common read/write paradigm, to be consistent across all tools, and to be compatible with Slicer3 (i.e. Write out Bspline transforms with the bulk transform given second).
-Increased code coverage by merging common functionality that existed in each tool separately into the BRAINSCommonLIb, thus removing code, and ensuring that exactly the same behavior was done across all tools.
-Worked around bug where ITK LSBFGSB optimizer does not properly return the number of iterations performed,  This improperly reported that Bspline registrations were never done when doing a Bspline registration alone.
-Improved source code documentation for commonly used functions.



git-svn-id: http://svn.slicer.org/Slicer4/trunk@13566 3bd1e089-480b-0410-8dfb-8563597acbee
parent 60b6b354
......@@ -9,77 +9,6 @@
#include "Air16_rw.h"
#include "ConvertToRigidAffine.h"
#if 0
// HACK: TODO: All the byte swapping code should use methods for itksys
// routines.
static void
FourByteSwap (void *const pntr)
{
/*NOTE: This can probably be speed up by using byte swapping */
const unsigned char b0 = *( (unsigned char *)(pntr) + 0 );
const unsigned char b1 = *( (unsigned char *)(pntr) + 1 );
const unsigned char b2 = *( (unsigned char *)(pntr) + 2 );
const unsigned char b3 = *( (unsigned char *)(pntr) + 3 );
*( (unsigned char *)(pntr) + 0 ) = b3;
*( (unsigned char *)(pntr) + 1 ) = b2;
*( (unsigned char *)(pntr) + 2 ) = b1;
*( (unsigned char *)(pntr) + 3 ) = b0;
}
static void
EightByteSwap (void *const pntr)
{
const unsigned char b0 = *( (unsigned char *)(pntr) + 0 );
const unsigned char b1 = *( (unsigned char *)(pntr) + 1 );
const unsigned char b2 = *( (unsigned char *)(pntr) + 2 );
const unsigned char b3 = *( (unsigned char *)(pntr) + 3 );
const unsigned char b4 = *( (unsigned char *)(pntr) + 4 );
const unsigned char b5 = *( (unsigned char *)(pntr) + 5 );
const unsigned char b6 = *( (unsigned char *)(pntr) + 6 );
const unsigned char b7 = *( (unsigned char *)(pntr) + 7 );
*( (unsigned char *)(pntr) + 0 ) = b7;
*( (unsigned char *)(pntr) + 1 ) = b6;
*( (unsigned char *)(pntr) + 2 ) = b5;
*( (unsigned char *)(pntr) + 3 ) = b4;
*( (unsigned char *)(pntr) + 4 ) = b3;
*( (unsigned char *)(pntr) + 5 ) = b2;
*( (unsigned char *)(pntr) + 6 ) = b1;
*( (unsigned char *)(pntr) + 7 ) = b0;
}
IPL_ENDIAN_TYPE
getMachineEndianess (void)
{
int check;
/*00000000 */
( (unsigned char *)(&check) )[0] = 0;
/*00000001 */
( (unsigned char *)(&check) )[1] = 1;
/*00000010 */
( (unsigned char *)(&check) )[2] = 2;
/*00000011 */
( (unsigned char *)(&check) )[3] = 3;
switch ( check )
{
/*00000000 00000001 00000010 0000011 -> 66051 Mips/SGI */
case IPL_BIG_ENDIAN:
return ( IPL_BIG_ENDIAN );
/*00000011 00000010 00000001 0000000 -> 5046297 Intel x86 Alpha AMD */
case IPL_LITTLE_ENDIAN:
return ( IPL_LITTLE_ENDIAN );
/*00000001 00000000 00000011 0000010 -> 1677986 Vax */
case IPL_MIXED_ENDIAN:
return ( IPL_MIXED_ENDIAN );
default:
return ( IPL_UNKNOWN_ENDIAN );
}
}
#endif
void swapAir16(struct AIR_Air16 *air)
{
......@@ -114,7 +43,7 @@ CrossOverAffineSystem<double, 3> ::AffineTransformPointer ReadAir16File (
if ( AirFilePointer == NULL )
{
std::cout << "Air file '" << airFile << "' couldn't be opened 'rb'."
<< std::endl;
<< std::endl;
exit(-1);
}
struct AIR_Air16 airBlk;
......@@ -124,7 +53,7 @@ CrossOverAffineSystem<double, 3> ::AffineTransformPointer ReadAir16File (
// Need to implement byteswapping code here.
const double F64_UINT_MAX = UINT_MAX;
if ( airBlk.s.bits > vcl_sqrt(F64_UINT_MAX) ) // copied this test from
// scanair source code.
// scanair source code.
{
swapAir16(&airBlk);
}
......@@ -146,7 +75,7 @@ CrossOverAffineSystem<double, 3> ::AffineTransformPointer ReadAir16File (
for ( int i = 0; i < 4; i++ )
{
Q_AIR.put(i, j, airBlk.e[j][i]); // NOTE: Need to transpose the "e"
// matrix
// matrix
}
}
......
......@@ -36,16 +36,16 @@
template <typename TImageType>
class AirAffineTransform
{
{
public:
typedef typename CrossOverAffineSystem<double,
3>::AffineTransformType TransformType;
3>::AffineTransformType TransformType;
typedef typename TransformType::Pointer
TransformPointer;
TransformPointer;
typedef typename TImageType::RegionType::SizeType
SizeType;
SizeType;
typedef typename TImageType::SpacingType
SpacingType;
SpacingType;
AirAffineTransform()
{
......@@ -65,138 +65,138 @@ public:
const std::string & StudyId);
int Write (const std::string & xfrmFile)
{
{
const std::string dummy("UNKNOWN");
return Write(xfrmFile, dummy, dummy);
}
}
TransformPointer GetAffineTransformPointer()
{
{
return m_AffineTransform;
}
}
void SetAffineTransformPointer(TransformPointer & xfrm)
{
{
m_AffineTransform = xfrm;
}
}
// Check to make sure that input image
// and transform characteristics match, but then this needs to be templated.
bool CheckMovingImageCharacteristics(
typename TImageType::Pointer & ImageToDeform) const
{
{
if ( ( ImageToDeform->GetLargestPossibleRegion().GetSize() !=
m_MovingImageSize )
|| ( ImageToDeform->GetSpacing() != m_MovingImageSpacing )
)
m_MovingImageSize )
|| ( ImageToDeform->GetSpacing() != m_MovingImageSpacing )
)
{
return false;
}
return true;
}
}
SizeType GetFixedImageSize()
{
{
return m_FixedImageSize;
}
}
SpacingType GetFixedImageSpacing()
{
{
return m_FixedImageSpacing;
}
}
// Brains2 XFRM always assume origin is 0,0,0
typename TImageType::PointType GetFixedImageOrigin()
{
{
typename TImageType::PointType tempOrigin;
tempOrigin[0] = tempOrigin[1] = tempOrigin[2] = 0;
return tempOrigin;
}
}
// Brains2 XFRM always assume direction is coronal;
// REFACTOR: In general, Air is probably agnostic about orientation.
typename TImageType::DirectionType GetFixedImageDirection()
{
{
typename itk::SpatialOrientationAdapter::DirectionType CORdir
= itk::SpatialOrientationAdapter().ToDirectionCosines(
itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP);
itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP);
return CORdir;
}
}
void SetFixedImageSize(SizeType & size)
{
{
m_FixedImageSize = size;
}
}
void SetFixedImageSpacing(SpacingType & spacing)
{
{
m_FixedImageSpacing = spacing;
}
}
void SetFixedAttributesFromImage(typename TImageType::Pointer & fixedImage)
{
{
SizeType size = fixedImage->GetLargestPossibleRegion().GetSize();
this->SetFixedImageSize(size);
SpacingType spacing = fixedImage->GetSpacing();
this->SetFixedImageSpacing(spacing);
}
}
SizeType GetMovingImageSize()
{
{
return m_MovingImageSize;
}
}
SpacingType GetMovingImageSpacing()
{
{
return m_MovingImageSpacing;
}
}
void SetMovingImageSize(SizeType & size)
{
{
m_MovingImageSize = size;
}
}
void SetMovingImageSpacing(SpacingType & spacing)
{
{
m_MovingImageSpacing = spacing;
}
}
void SetMovingAttributesFromImage(typename TImageType::Pointer & movingImage)
{
{
SizeType size = movingImage->GetLargestPossibleRegion().GetSize();
this->SetMovingImageSize(size);
SpacingType spacing = movingImage->GetSpacing();
this->SetMovingImageSpacing(spacing);
}
}
typename TImageType::Pointer
ApplyTransform(typename TImageType::Pointer & ImageToDeform)
{
const bool ValidInput = this->CheckMovingImageCharacteristics(ImageToDeform);
if ( !ValidInput )
ApplyTransform(typename TImageType::Pointer & ImageToDeform)
{
return TImageType::New(); // Return a null pointer
const bool ValidInput = this->CheckMovingImageCharacteristics(ImageToDeform);
if ( !ValidInput )
{
return TImageType::New(); // Return a null pointer
}
typedef typename itk::ResampleImageFilter<TImageType,
TImageType> ResampleFilterType;
typename ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(m_AffineTransform);
resample->SetInput(ImageToDeform);
resample->SetSize( this->GetFixedImageSize() );
resample->SetOutputOrigin( this->GetFixedImageOrigin() );
resample->SetOutputSpacing( this->GetFixedImageSpacing() );
resample->SetOutputDirection( this->GetFixedImageDirection() );
typename TImageType::PixelType p(
static_cast<typename TImageType::PixelType>( 0 ) );
resample->SetDefaultPixelValue(p);
resample->Update();
return resample->GetOutput();
}
typedef typename itk::ResampleImageFilter<TImageType,
TImageType> ResampleFilterType;
typename ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(m_AffineTransform);
resample->SetInput(ImageToDeform);
resample->SetSize( this->GetFixedImageSize() );
resample->SetOutputOrigin( this->GetFixedImageOrigin() );
resample->SetOutputSpacing( this->GetFixedImageSpacing() );
resample->SetOutputDirection( this->GetFixedImageDirection() );
typename TImageType::PixelType p(
static_cast<typename TImageType::PixelType>( 0 ) );
resample->SetDefaultPixelValue(p);
resample->Update();
return resample->GetOutput();
}
private:
TransformPointer m_AffineTransform;
......@@ -204,14 +204,14 @@ private:
SpacingType m_FixedImageSpacing;
SizeType m_MovingImageSize;
SpacingType m_MovingImageSpacing;
};
};
// HACK: The following function should be put into the previous class, and all
// code should be modified to reflect the simplified interface.
template <typename TImageType>
int
AirAffineTransform<TImageType>::
Write (const std::string & AirFile,
Write (const std::string & AirFile,
const std::string & /* PatientId */,
const std::string & /* StudyId */)
{
......@@ -220,9 +220,9 @@ AirAffineTransform<TImageType>::
if ( fp == NULL )
{
std::cerr
<<
"Failed to open air file for writing -- is this a file or directory protection issue? "
<< __FILE__ << __LINE__ << std::endl;
<<
"Failed to open air file for writing -- is this a file or directory protection issue? "
<< __FILE__ << __LINE__ << std::endl;
return -1;
}
......@@ -264,29 +264,10 @@ AirAffineTransform<TImageType>::
for ( int j = 0; j < 4; j++ )
{
airBlk.e[j][i] = Q_AIR.get(i, j); // NOTE: Need to transpose the "e"
// air matrix
// air matrix
}
}
#if 0
double pixel_size_s = airBlk.s.x_size;
if ( airBlk.s.y_size < pixel_size_s )
{
pixel_size_s = airBlk.s.y_size;
}
if ( airBlk.s.z_size < pixel_size_s )
{
pixel_size_s = airBlk.s.z_size;
}
for ( int i = 0; i < 4; i++ )
{
airBlk.e[0][i] /= (F64)(airBlk.s.x_size / pixel_size_s);
airBlk.e[1][i] /= (F64)(airBlk.s.y_size / pixel_size_s);
airBlk.e[2][i] /= (F64)(airBlk.s.z_size / pixel_size_s);
}
#endif
// All ipl .air16 files need to be BigEndian.
//
......@@ -302,15 +283,15 @@ AirAffineTransform<TImageType>::
template <typename TImageType>
int AirAffineTransform<TImageType>::
Read (std::string xfrmFile)
Read (std::string xfrmFile)
{
FILE *transformFilePointer = fopen (xfrmFile.c_str(), "rb");
if ( transformFilePointer == NULL )
{
std::cerr
<< "Failed to open Air file for reading -- does your file exist as named? "
<< __FILE__ << __LINE__ << std::endl;
<< "Failed to open Air file for reading -- does your file exist as named? "
<< __FILE__ << __LINE__ << std::endl;
return -1;
}
struct AIR_Air16 airBlk;
......@@ -321,16 +302,16 @@ int AirAffineTransform<TImageType>::
else
{
std::cerr
<<
"Failed to read a whole air transform from Air file -- was your file really an Air file? "
<< __FILE__ << __LINE__ << std::endl;
<<
"Failed to read a whole air transform from Air file -- was your file really an Air file? "
<< __FILE__ << __LINE__ << std::endl;
return -1;
}
// Need to implement byteswapping code here.
const double F64_UINT_MAX = UINT_MAX;
if ( airBlk.s.bits > vcl_sqrt(F64_UINT_MAX) ) // copied this test from
// scanair source code.
// scanair source code.
{
swapAir16(&airBlk);
}
......@@ -355,25 +336,6 @@ int AirAffineTransform<TImageType>::
this->m_MovingImageSpacing[1] = airBlk.r.y_size;
this->m_MovingImageSpacing[2] = airBlk.r.z_size;
#if 0
double pixel_size_s = airBlk.s.x_size;
if ( airBlk.s.y_size < pixel_size_s )
{
pixel_size_s = airBlk.s.y_size;
}
if ( airBlk.s.z_size < pixel_size_s )
{
pixel_size_s = airBlk.s.z_size;
}
for ( int i = 0; i < 4; i++ )
{
airBlk.e[0][i] *= ( airBlk.s.x_size / pixel_size_s );
airBlk.e[1][i] *= ( airBlk.s.y_size / pixel_size_s );
airBlk.e[2][i] *= ( airBlk.s.z_size / pixel_size_s );
}
#endif
// Transfer values from air structures.
// Transfer to intermediate 4x4 matrix so that it can be checked easily in
......@@ -385,7 +347,7 @@ int AirAffineTransform<TImageType>::
for ( int i = 0; i < 4; i++ )
{
Q_AIR.put(i, j, airBlk.e[j][i]); // NOTE: Need to transpose the "e"
// matrix
// matrix
}
}
......
......@@ -6,83 +6,83 @@
namespace itk
{
/** \class ApplicationBase
*
* This class ties together an input parser, a preprocessor,
* a registrator components to
* form a deformable registration/atlas segmentation application.
*
*/
template <typename TParser,
typename TPreprocessor,
typename TRegistrator>
class ApplicationBase : public Object
{
public:
/** Standard class typedefs. */
typedef ApplicationBase Self;
typedef Object Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Run-time type information (and related methods). */
itkTypeMacro(MIMApplication, Object);
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Input parser type. */
typedef TParser ParserType;
typedef typename ParserType::Pointer ParserPointer;
/** Preprocessor type. */
typedef TPreprocessor PreprocessorType;
typedef typename PreprocessorType::Pointer PreprocessorPointer;
/** Registrator type. */
typedef TRegistrator RegistratorType;
typedef typename RegistratorType::Pointer RegistratorPointer;
/**Set Debug mode*/
itkSetMacro(OutDebug, bool);
itkGetConstMacro(OutDebug, bool);
RegistratorType * GetRegistratorType(void)
{
return m_Registrator;
}
/** Execute the application. */
virtual void Execute();
protected:
ApplicationBase();
virtual ~ApplicationBase()
{}
/** Initialize the input parser. */
virtual void InitializeParser()
{}
/*** Initialize the preprocessor */
virtual void InitializePreprocessor()
{}
/*** Initialize the registrator */
virtual void InitializeRegistrator()
{}
ParserPointer m_Parser;
PreprocessorPointer m_Preprocessor;
RegistratorPointer m_Registrator;
bool m_OutDebug;
};
/** \class ApplicationBase
*
* This class ties together an input parser, a preprocessor,
* a registrator components to
* form a deformable registration/atlas segmentation application.
*
*/
template <typename TParser,
typename TPreprocessor,
typename TRegistrator>
class ApplicationBase : public Object
{
public:
/** Standard class typedefs. */
typedef ApplicationBase Self;
typedef Object Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
/** Run-time type information (and related methods). */
itkTypeMacro(MIMApplication, Object);
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Input parser type. */
typedef TParser ParserType;
typedef typename ParserType::Pointer ParserPointer;
/** Preprocessor type. */
typedef TPreprocessor PreprocessorType;
typedef typename PreprocessorType::Pointer PreprocessorPointer;
/** Registrator type. */
typedef TRegistrator RegistratorType;
typedef typename RegistratorType::Pointer RegistratorPointer;
/**Set Debug mode*/
itkSetMacro(OutDebug, bool);
itkGetConstMacro(OutDebug, bool);
RegistratorType * GetRegistratorType(void)
{
return m_Registrator;
}
/** Execute the application. */
virtual void Execute();
protected:
ApplicationBase();
virtual ~ApplicationBase()
{}
/** Initialize the input parser. */
virtual void InitializeParser()
{}
/*** Initialize the preprocessor */
virtual void InitializePreprocessor()
{}