Commit 83939e01 authored by hjohnson's avatar hjohnson
Browse files

ENH: First submission to include the suite of BRAINS registraiton tools to...

ENH:  First submission to include the suite of BRAINS registraiton tools to Slicer3.  All the test tree has been removed (due to the large size of the regression test data that is usually included.

git-svn-id: http://svn.slicer.org/Slicer4/trunk@12920 3bd1e089-480b-0410-8dfb-8563597acbee
parent 0dd08054
#include <itkObjectFactory.h>
#include <itkObject.h>
#include <vnl/vnl_matrix_fixed.h>
#include <vnl/vnl_vector.h>
#include <iostream>
#include <fstream>
#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)
{
FourByteSwap ( &( air->s.bits ) );
FourByteSwap ( &( air->s.x_dim ) );
FourByteSwap ( &( air->s.y_dim ) );
FourByteSwap ( &( air->s.z_dim ) );
EightByteSwap ( &( air->s.x_size ) );
EightByteSwap ( &( air->s.y_size ) );
EightByteSwap ( &( air->s.z_size ) );
FourByteSwap ( &( air->r.bits ) );
FourByteSwap ( &( air->r.x_dim ) );
FourByteSwap ( &( air->r.y_dim ) );
FourByteSwap ( &( air->r.z_dim ) );
EightByteSwap ( &( air->r.x_size ) );
EightByteSwap ( &( air->r.y_size ) );
EightByteSwap ( &( air->r.z_size ) );
for ( int i = 0; i < 4; i++ )
{
for ( int j = 0; j < 4; j++ )
{
EightByteSwap ( &( air->e[i][j] ) );
}
}
}
CrossOverAffineSystem<double, 3> ::AffineTransformPointer ReadAir16File (
std::string airFile)
{
FILE *AirFilePointer = fopen(airFile.c_str(), "rb");
if ( AirFilePointer == NULL )
{
std::cout << "Air file '" << airFile << "' couldn't be opened 'rb'."
<< std::endl;
exit(-1);
}
struct AIR_Air16 airBlk;
fread(&airBlk, sizeof( struct AIR_Air16 ), 1, AirFilePointer);
fclose(AirFilePointer);
// 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.
{
swapAir16(&airBlk);
}
if ( ( airBlk.s.x_dim <= 0 ) || ( airBlk.s.x_dim > 2048 ) )
{
swapAir16(&airBlk);
fprintf(
stderr,
"Apparently, the s.bits field was garbage, following the brains2 heuristic on the range of x_dim.\n");
}
// Transfer values from air structures.
// Transfer to intermediate 4x4 matrix so that it can be checked easily in
// octave or matlab
CrossOverAffineSystem<double, 3> ::VnlTransformMatrixType44 Q_AIR(0.0);
for ( int j = 0; j < 4; j++ )
{
for ( int i = 0; i < 4; i++ )
{
Q_AIR.put(i, j, airBlk.e[j][i]); // NOTE: Need to transpose the "e"
// matrix
}
}
std::cout << "=======================Q_AIR input\n" << Q_AIR;
CrossOverAffineSystem<double,3> ::AffineTransformPointer result
= CrossOverAffineSystem<double,3> ::AffineTransformType::New();
AssignRigid::AssignConvertedTransform( result, Q_AIR );
return result;
}
#ifndef __AirAFFINE_RW_H__
#define __AirAFFINE_RW_H__
#include <itkMatrix.h>
#include <itkVector.h>
#include <itkAffineTransform.h>
#include <vnl/vnl_matrix.h>
#include <vnl/vnl_inverse.h>
#include <vnl/vnl_matrix_fixed.h>
#include <vnl/vnl_vector.h>
#include <vcl_cmath.h>
#include <stdio.h>
#include "iplUtils/iplUtils.h"
#include "iplUtils/b2Version.h"
#include "iplUtils/iplSystem.h"
#include "iplHeader/iplHeader.h"
#include <string>
#include "iplHeader/iplHeaderTypes.h"
#include "iplFilterBase/iplFilterBase.h"
#include "iplUtils/iplUtils.h"
#include "InsightToolkit/Utilities/itksys/SystemTools.hxx"
#include "itkImage.h"
#include "itkSpatialOrientation.h"
#include "itkSpatialOrientationAdapter.h"
#include "itkResampleImageFilter.h"
#include "CrossOverAffineSystem.h"
#include "ConvertToRigidAffine.h"
#ifndef __RPW_AIR_H__
#define __RPW_AIR_H__
#include <AIR.h>
#endif
#include "Air16_rw.h"
template <typename TImageType>
class AirAffineTransform
{
public:
typedef typename CrossOverAffineSystem<double,
3>::AffineTransformType TransformType;
typedef typename TransformType::Pointer
TransformPointer;
typedef typename TImageType::RegionType::SizeType
SizeType;
typedef typename TImageType::SpacingType
SpacingType;
AirAffineTransform()
{
// as a default to indicate failure, ReadAirAffineFile will return an IsNull
// Pointer.
m_AffineTransform = TransformType::New();
}
int Read (const std::string xfrmFile);
// HACK: Need corresponding Set functions, and then the write command needs
// to be pulled
// into the class with the folowing interface:
// void WriteAirAffineFile(const std::string xfrmFile);
int Write (const std::string & xfrmFile,
const std::string & PatientId,
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 )
)
{
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);
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 )
{
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();
}
private:
TransformPointer m_AffineTransform;
SizeType m_FixedImageSize;
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,
const std::string & PatientId,
const std::string & StudyId)
{
FILE *fp = fopen (AirFile.c_str(), "wb");
if ( fp == NULL )
{
std::cerr
<<
"Failed to open air file for writing -- is this a file or directory protection issue? "
<< __FILE__ << __LINE__ << std::endl;
return -1;
}
// Re-useable
struct AIR_Air16 airBlk;
strncpy (airBlk.s_file, "<fixed input file name here>", 127);
strncpy (airBlk.r_file, "<moving input file name here>", 127);
strncpy (airBlk.comment, "NoCommentProvided", 127);
strncpy (airBlk.reserved, "", 115);
airBlk.s.bits = 8; // Just assume unsigned char data.
airBlk.s.x_dim = this->m_FixedImageSize[0];
airBlk.s.y_dim = this->m_FixedImageSize[1];
airBlk.s.z_dim = this->m_FixedImageSize[2];
airBlk.s.x_size = this->m_FixedImageSpacing[0];
airBlk.s.y_size = this->m_FixedImageSpacing[1];
airBlk.s.z_size = this->m_FixedImageSpacing[2];
airBlk.r.bits = 8; // Just assume unsigned char data.
airBlk.r.x_dim = this->m_MovingImageSize[0];
airBlk.r.y_dim = this->m_MovingImageSize[1];
airBlk.r.z_dim = this->m_MovingImageSize[2];
airBlk.r.x_size = this->m_MovingImageSpacing[0];
airBlk.r.y_size = this->m_MovingImageSpacing[1];
airBlk.r.z_size = this->m_MovingImageSpacing[2];
airBlk.s_hash = 0;
airBlk.r_hash = 0;
airBlk.s_volume = 0;
airBlk.r_volume = 0;
// Transfer to intermediate 4x4 matrix so that it can be checked easily in
// octave or matlab
typename CrossOverAffineSystem<double, 3> ::VnlTransformMatrixType44 Q_AIR(
0.0);
AssignRigid::AssignConvertedTransform( Q_AIR, m_AffineTransform );
std::cout << "=======================Q_AIR output\n" << Q_AIR;
// Transfer values into air structures.
for ( int i = 0; i < 4; i++ )
{
for ( int j = 0; j < 4; j++ )
{
airBlk.e[j][i] = Q_AIR.get(i, j); // NOTE: Need to transpose the "e"
// 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.
//
if ( getMachineEndianess () != IPL_BIG_ENDIAN )
{
swapAir16(&airBlk);
}
fwrite(&airBlk, sizeof( struct AIR_Air16 ), 1, fp);
fclose(fp);
return 0;
}
template <typename TImageType>
int AirAffineTransform<TImageType>::
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;
return -1;
}
struct AIR_Air16 airBlk;
if ( fread(&airBlk, sizeof( struct AIR_Air16 ), 1, transformFilePointer) == 1 )
{
fclose(transformFilePointer);
}
else
{
std::cerr
<<
"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.
{
swapAir16(&airBlk);
}
if ( ( airBlk.s.x_dim <= 0 ) || ( airBlk.s.x_dim > 2048 ) )
{
swapAir16(&airBlk);
fprintf(
stderr,
"Apparently, the s.bits field was garbage, following the brains2 heuristic on the range of x_dim.\n");
}
this->m_FixedImageSize[0] = airBlk.s.x_dim;
this->m_FixedImageSize[1] = airBlk.s.y_dim;
this->m_FixedImageSize[2] = airBlk.s.z_dim;
this->m_FixedImageSpacing[0] = airBlk.s.x_size;
this->m_FixedImageSpacing[1] = airBlk.s.y_size;
this->m_FixedImageSpacing[2] = airBlk.s.z_size;
this->m_MovingImageSize[0] = airBlk.r.x_dim;
this->m_MovingImageSize[1] = airBlk.r.y_dim;
this->m_MovingImageSize[2] = airBlk.r.z_dim;
this->m_MovingImageSpacing[0] = airBlk.r.x_size;
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
// octave or matlab
CrossOverAffineSystem<double, 3> ::VnlTransformMatrixType44 Q_AIR(0.0);
for ( int j = 0; j < 4; j++ )
{
for ( int i = 0; i < 4; i++ )
{
Q_AIR.put(i, j, airBlk.e[j][i]); // NOTE: Need to transpose the "e"
// matrix
}
}
std::cout << "=======================Q_AIR input\n" << Q_AIR;
AssignRigid::AssignConvertedTransform( m_AffineTransform, Q_AIR );
return 0;
}
#endif
#ifndef _ApplicationBase_h
#define _ApplicationBase_h
#include "itkObjectFactory.h"
#include "itkObject.h"
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);