Updates will be applied on October 27th between 12pm - 12:45pm EDT (UTC-0400). Gitlab may be slow during the maintenance window.

Commit c064c934 authored by hjohnson's avatar hjohnson
Browse files

ENH: Merged recent changes contributed by Slicer4 developers into Slicer3 and BRAINS3.

git-svn-id: http://svn.slicer.org/Slicer4/trunk@15625 3bd1e089-480b-0410-8dfb-8563597acbee
parent 51c08257
#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 __BRAINSFitBSpline_h
#define __BRAINSFitBSpline_h
#include <itkMattesMutualInformationImageToImageMetric.h>
#include <itkImageToImageMetric.h>
#include <itkBSplineDeformableTransform.h>
#include "itkBSplineDeformableTransformInitializer.h"
#include <itkLBFGSBOptimizer.h>
......@@ -10,7 +11,7 @@
#include <itkLinearInterpolateImageFunction.h>
#include <itkResampleImageFilter.h>
#include "itkMultiModal3DMutualRegistrationHelper.h"
#include "genericRegistrationHelper.h"
/**
* This class is the BSpline component of the BRAINSFit program developed at
......@@ -27,8 +28,8 @@ DoBSpline(typename BSplineTransformType::Pointer InitializerBsplineTransform,
typename ImageMaskSpatialObjectType::Pointer m_FixedMask,
typename ImageMaskSpatialObjectType::Pointer m_MovingMask,
const int m_NumberOfSamples,
const bool m_UseCachingOfBSplineWeights,
const bool m_UseExplicitPDFDerivatives,
typename itk::ImageToImageMetric<
RegisterImageType, RegisterImageType >::Pointer CostMetricObject,
const double m_MaxBSplineDisplacement,
const float m_CostFunctionConvergenceFactor,
const float m_ProjectedGradientTolerance,
......@@ -59,9 +60,8 @@ DoBSpline(typename BSplineTransformType::Pointer InitializerBsplineTransform,
typedef typename itk::LBFGSBOptimizer OptimizerType;
typedef typename itk::MattesMutualInformationImageToImageMetric<