/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
#include <fstream>
// Include for "rusage"
-#include <ctime>
+#include <ctime>
#if defined(unix) || defined(__APPLE__)
# include <sys/time.h>
# include <sys/resource.h>
typedef unsigned char uchar;
typedef unsigned short ushort;
typedef unsigned int uint;
-
+
#define CLITK_TRY_CATCH_EXIT(func) \
try { \
func; \
std::cout << "Unknown excpetion" << std::endl; \
exit(-3); \
}
-
+
//--------------------------------------------------------------------
// when everything goes wrong
#define WHEREAMI "[ " << __FILE__ << " ] line " << __LINE__
-#define FATAL(a) std::cerr << "ERROR in " << WHEREAMI << ": " << a; exit(0);
-
+#define FATAL(a) { std::cerr << "ERROR in " << WHEREAMI << ": " << a; exit(0); }
+
//--------------------------------------------------------------------
// GGO with modified struct name
#define GGO(ggo_filename, args_info) \
cmdline_parser_##ggo_filename##2(argc, argv, &args_info, 1, 1, 0); \
if (args_info.config_given) \
cmdline_parser_##ggo_filename##_configfile (args_info.config_arg, &args_info, 0, 0, 1); \
- else cmdline_parser_##ggo_filename(argc, argv, &args_info);
+ else cmdline_parser_##ggo_filename(argc, argv, &args_info);
//--------------------------------------------------------------------
// skip line with #
//--------------------------------------------------------------------
// Return filename extension
std::string GetExtension(const std::string& filename);
-
+
//--------------------------------------------------------------------
// Convert float, double ... to string
template<class T> std::string toString(const T & t);
template<class T> std::string toStringVector(const T * t, const int n);
template<class T> std::string toStringVector(const T & t, const int n);
template<class T> std::string toStringVector(const std::vector<T> & t);
- template <class T> bool fromString(T& t,
- const std::string& s,
+ template <class T> bool fromString(T& t,
+ const std::string& s,
std::ios_base& (*f)(std::ios_base&)=std::dec);
//--------------------------------------------------------------------
// Return the name of a type as a string
template<class TPixel>
std::string GetTypeAsString();
-
+
//--------------------------------------------------------------------
// Convert radian / degree
double rad2deg(double anglerad);
std::string CreateListOfTypes(bool last=true) {
return GetTypeAsString<T1>();
}
-
+
template<class T1, class T2>
std::string CreateListOfTypes(bool last=true) {
if (last) return CreateListOfTypes<T1>()+" and "+CreateListOfTypes<T2>();
else return CreateListOfTypes<T1>()+", "+CreateListOfTypes<T2>();
}
-
+
template<class T1, class T2, class T3>
std::string CreateListOfTypes(bool last=true) {
if (last) return CreateListOfTypes<T1,T2>(false)+" and "+CreateListOfTypes<T3>();
else return CreateListOfTypes<T1,T2,T3,T4,T5,T6,T7>(false)+", "+CreateListOfTypes<T8>();
}
//--------------------------------------------------------------------
-
+
//--------------------------------------------------------------------
void FindAndReplace(std::string & line, const std::string & tofind, const std::string & replacement);
void FindAndReplace(std::string & line, const std::vector<std::string> & tofind, const std::vector<std::string> & toreplace);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
- double ComputeEuclideanDistanceFromPointToPlane(const itk::ContinuousIndex<double, 3> point,
- const itk::ContinuousIndex<double, 3> pointInPlane,
+ double ComputeEuclideanDistanceFromPointToPlane(const itk::ContinuousIndex<double, 3> point,
+ const itk::ContinuousIndex<double, 3> pointInPlane,
const itk::ContinuousIndex<double, 3> normalPlane);
//--------------------------------------------------------------------
//--------------------------------------------------------------------
// Convert a map to a vector
- template <typename M, typename V>
+ template <typename M, typename V>
void MapToVecFirst(const M & m, V & v);
- template <typename M, typename V>
+ template <typename M, typename V>
void MapToVecSecond(const M & m, V & v);
//--------------------------------------------------------------------
} // end namespace
#endif /* end #define CLITKCOMMON_H */
-
template<class T>
void _print_container(T const& a)
{ for(typename T::const_iterator i=a.begin();i!=a.end();++i) { std::cout << *i << " "; };}
-#define DDS(a) { std::cout << #a " = [ "; _print_container(a) ; std::cout << " ]" << std::endl;std::cout.flush():}
+#define DDS(a) { std::cout << #a " = [ "; _print_container(a) ; std::cout << " ]" << std::endl;std::cout.flush();}
#endif
mElapsed += (mEnd.ru_utime.tv_usec - mBegin.ru_utime.tv_usec)+
(mEnd.ru_utime.tv_sec - mBegin.ru_utime.tv_sec)*1000000;
}
+ else
#elif defined(_WIN32)
QueryPerformanceCounter((LARGE_INTEGER*)&mEnd);
if (accumulate) {
mElapsed += ((mEnd-mBegin)*1000000)/(long double)mFrequency;
}
+ else
#endif
- else {
+ {
mNumberOfCall--;
}
}
// #endif // If UNIX
#endif /* end #define CLITKTIMER_CXX */
-
typedef typename OutputImageType::PixelType DisplacementType;
DisplacementType displacement;
inputIt.GoToBegin();
-
+
typename OutputImageType::SizeType size = outputPtr->GetLargestPossibleRegion().GetSize();
//define some temp variables
overlap *= 1.0 - distance[dim];
}
upper >>= 1;
-
+
if (neighIndex[dim] >= size[dim])
neighIndex[dim] = size[dim] - 1;
}
//Empty constructor
template<class InputImageType, class OutputImageType > HelperClass2<InputImageType, OutputImageType>::HelperClass2()
{
- m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
+ PixelType zero;
+ for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
+ m_EdgePaddingValue=zero;
+ //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero;
}
#endif
{
// std::cout << "HelperClass2::ThreadedGenerateData - IN " << threadId << std::endl;
-
+
//Get pointer to the input
typename InputImageType::ConstPointer inputPtr = this->GetInput();
++inputIt;
}//end while
-
+
// std::cout << "HelperClass2::ThreadedGenerateData - OUT " << threadId << std::endl;
-
+
}//end member
template <class InputImageType, class OutputImageType>
InvertVFFilter<InputImageType, OutputImageType>::InvertVFFilter()
{
- m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
+
+ //m_EdgePaddingValue=itk::NumericTraits<PixelType>::Zero; //no other reasonable value?
+ PixelType zero;
+ for(unsigned int i=0;i <PixelType::Dimension; i++) zero[i] = 0.0;
+ m_EdgePaddingValue=zero; //no other reasonable value?
+
m_ThreadSafe=false;
m_Verbose=false;
}
//Set the output
this->SetNthOutput(0, helper2->GetOutput());
-
+
//std::cout << "InvertVFFilter::GenerateData - OUT" << std::endl;
}
filter->SetIOVerbose(args_info.verbose_flag);
filter->SetOutputFilename(output);
filter->SetVV(args_info.vv_flag);
+ filter->SetCorrectNegativeSpacingFlag(args_info.correct_flag);
filter->EnableWriteCompression(args_info.compression_flag);
if (args_info.type_given) filter->SetOutputPixelType(args_info.type_arg);
//-------------------------------------------------------------------=
#endif /* end #define CLITKIMAGECONVERT_CXX */
-
option "verbose" v "Verbose" flag off
option "compression" c "Compress output" flag off
option "vv" - "Read image as in vv and save transform in meta information" flag off
+option "correct" - "Correct dicom with negative Z spacing" flag off
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/
+
#ifndef CLITKIMAGECONVERTGENERICFILTER_CXX
#define CLITKIMAGECONVERTGENERICFILTER_CXX
#include "clitkImageConvertGenericFilter.h"
#include "vvImageReader.h"
#include "vvImageWriter.h"
+#include "itkFlipImageFilter.h"
+#include "itkGDCMImageIO.h"
+
+#include "gdcmReader.h"
+#include "gdcmAttribute.h"
+#include "gdcmPrinter.h"
+#include "gdcmDict.h"
+#include "gdcmStringFilter.h"
//--------------------------------------------------------------------
clitk::ImageConvertGenericFilter::ImageConvertGenericFilter():
mDisplayWarning = true;
mWarning = "";
mWarningOccur = false;
-
+ SetCorrectNegativeSpacingFlag(false);
+
InitializeImageType<2>();
InitializeImageType<3>();
InitializeImageType<4>();
return;
}
else if ((m_PixelTypeName == mOutputPixelTypeName) || (mOutputPixelTypeName == "NotSpecified")) {
+
+ // Get input image
typename InputImageType::Pointer input = this->template GetInput<InputImageType>(0);
- this->SetNextOutput<InputImageType>(input);
+
+ if (mCorrectNegativeSpacingFlag) {
+ // Read dicom
+ gdcm::Reader reader;
+ reader.SetFileName(m_InputFilenames[0].c_str());
+ // if (!reader.CanReadFile(m_InputFilenames[0])) {
+ // std::cout << "Error: " << m_InputFilenames[0] << " is not a dicom file. Abort." << std::endl;
+ // exit(0);
+ // }
+ reader.Read();
+
+ // the dataset is the the set of element we are interested in:
+ gdcm::DataSet & ds = reader.GetFile().GetDataSet();
+
+ // Read the attribute SpacingBetweenSlices, check if negative and replace
+ gdcm::Attribute<0x0018,0x0088> SpacingBetweenSlices;
+ SpacingBetweenSlices.SetFromDataSet(ds);
+ double s = SpacingBetweenSlices.GetValue();
+ if (s >=0) {
+ std::cout << "Error: no negative spacing found SpacingBetweenSlices = " << s << " Abort. " << std::endl;
+ exit(0);
+ }
+ s = -s;
+
+ // Set spacing
+ typename InputImageType::SpacingType spacing = input->GetSpacing();
+ spacing[2] = s;
+ input->SetSpacing(spacing);
+
+ // Flip
+ typedef itk::FlipImageFilter< InputImageType > FilterType;
+ typename FilterType::Pointer filter = FilterType::New();
+ typedef typename FilterType::FlipAxesArrayType FlipAxesArrayType;
+ FlipAxesArrayType flipArray;
+ flipArray[0] = false;
+ flipArray[1] = false;
+ flipArray[2] = true;
+ filter->SetFlipAxes(flipArray);
+ filter->SetInput(input);
+ filter->Update();
+
+ // Read the attribute Image Position (Patient)
+ gdcm::Tag DetectorInformationSequenceTag(0x0054,0x0022);
+ const gdcm::DataElement & DIS = ds.GetDataElement(DetectorInformationSequenceTag);
+ gdcm::SmartPointer<gdcm::SequenceOfItems> sqf = DIS.GetValueAsSQ();
+ gdcm::Item & item = sqf->GetItem(1);
+ gdcm::DataSet & ds_position = item.GetNestedDataSet();
+ gdcm::Attribute<0x0020,0x0032> ImagePositionPatient;
+ ImagePositionPatient.SetFromDataSet(ds_position);
+ double x = ImagePositionPatient.GetValue(0);
+ double y = ImagePositionPatient.GetValue(1);
+ double z = ImagePositionPatient.GetValue(2);
+
+ // Set offset
+ typename InputImageType::PointType origin = input->GetOrigin();
+ origin[0] = x;
+ origin[1] = y;
+ origin[2] = z;
+ input->SetOrigin(origin);
+
+ // Orientation
+ typename InputImageType::DirectionType direction = input->GetDirection();
+ direction[2][2] = -1;
+ input->SetDirection(direction);
+
+ // Empty meta info
+ itk::MetaDataDictionary dict;// = new itk::MetaDataDictionary;
+ input->SetMetaDataDictionary(dict);
+
+ this->SetNextOutput<InputImageType>(input);
+ }
+
} else {
- // "trick" to call independent versions of update according to the
+ // "trick" to call independent versions of update according to the
// pixel type (vector or scalar), using partial specializations
- if (!UpdateWithSelectiveOutputType<InputImageType, ImageConvertTraits<typename InputImageType::PixelType>::IS_VECTOR>::Run(*this, mOutputPixelTypeName))
+ if (!UpdateWithSelectiveOutputType<InputImageType, ImageConvertTraits<typename InputImageType::PixelType>::IS_VECTOR>::Run(*this, mOutputPixelTypeName))
exit(-1);
}
}
#endif /* end #define CLITKIMAGECONVERTGENERICFILTER_CXX */
-
/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
* @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
* @date 05 May 2008 10:40:24
- * @brief
+ * @brief
===================================================================*/
namespace clitk {
-
+
template <class TPixel>
class ImageConvertTraits
{
}
};
- template < class TPixel, unsigned int Comp >
- class ImageConvertTraits< itk::Vector<TPixel, Comp> >
- {
- public:
- enum { IS_VECTOR = true };
+ template < class TPixel, unsigned int Comp >
+ class ImageConvertTraits< itk::Vector<TPixel, Comp> >
+ {
+ public:
+ enum { IS_VECTOR = true };
};
- class ImageConvertGenericFilter:
+ class ImageConvertGenericFilter:
public clitk::ImageToImageGenericFilter<ImageConvertGenericFilter> {
-
- public:
+
+ public:
// constructor - destructor
ImageConvertGenericFilter();
// New
itkNewMacro(Self);
-
+
// Members functions
std::string GetInputPixelTypeName() { return m_PixelTypeName; }
std::string GetOutputPixelTypeName() { return mOutputPixelTypeName; }
bool IsWarningOccur() { return mWarningOccur; }
std::string & GetWarning() { return mWarning; }
void EnableDisplayWarning(bool b) { mDisplayWarning = b; }
+ void SetCorrectNegativeSpacingFlag(bool b) { mCorrectNegativeSpacingFlag = b; }
//--------------------------------------------------------------------
// Main function called each time the filter is updated
- template<class InputImageType>
+ template<class InputImageType>
void UpdateWithInputImageType();
template<class PixelType, class OutputPixelType>
void CheckTypes(std::string inType, std::string outType);
-
+
protected:
template<unsigned int Dim> void InitializeImageType();
bool mWarningOccur;
bool mDisplayWarning;
bool mVV;
+ bool mCorrectNegativeSpacingFlag;
private:
template <class InputImageType, bool isVector>
else
{
std::string list = CreateListOfTypes<float, double>();
- std::cerr << "Error, I don't know the vector output type '" << outputPixelType
+ std::cerr << "Error, I don't know the vector output type '" << outputPixelType
<< "'. " << std::endl << "Known types are " << list << "." << std::endl;
return false;
}
-
+
return true;
}
private:
-
- template <class OutputPixelType>
+
+ template <class OutputPixelType>
static void UpdateWithOutputType(ImageConvertGenericFilter& filter)
{
// Read
filter.SetNextOutput<OutputImageType>(cast_filter->GetOutput());
}
};
-
+
template <class InputImageType>
class UpdateWithSelectiveOutputType<InputImageType, true>
{
static bool Run(ImageConvertGenericFilter& filter, std::string outputPixelType)
{
/*
- // RP: future conversions?
+ // RP: future conversions?
if (IsSameType<char>(outputPixelType))
UpdateWithOutputVectorType<char>();
else if (IsSameType<uchar>(outputPixelType))
UpdateWithOutputVectorType<ushort>();
else if (IsSameType<int>(outputPixelType))
UpdateWithOutputVectorType<int>();
- else
+ else
*/
if (IsSameType<float>(outputPixelType))
UpdateWithOutputVectorType<float>(filter);
else
{
std::string list = CreateListOfTypes<float, double>();
- std::cerr << "Error, I don't know the vector output type '" << outputPixelType
+ std::cerr << "Error, I don't know the vector output type '" << outputPixelType
<< "'. " << std::endl << "Known types are " << list << "." << std::endl;
return false;
}
-
+
return true;
}
-
+
private:
-
- template <class OutputPixelType>
+
+ template <class OutputPixelType>
static void UpdateWithOutputVectorType(ImageConvertGenericFilter& filter)
{
// Read
// Warning
filter.CheckTypes<PixelType, OutputPixelType>(filter.GetInputPixelTypeName(), filter.GetOutputPixelTypeName());
-
+
// Cast
typedef itk::Image<itk::Vector<OutputPixelType, InputImageType::PixelType::Dimension>, InputImageType::ImageDimension> OutputImageType;
typedef itk::VectorCastImageFilter<InputImageType, OutputImageType> FilterType;
} // end namespace
#endif /* end #define CLITKIMAGECONVERTGENERICFILTER_H */
-
}
else {
- std::cerr << "Mask image has a different size/spacing than input. Abort" << std::endl;
+ std::cerr << "Mask image has a different size/spacing than input. Abort. (Use option to resize)" << std::endl;
exit(-1);
}
}
labelImage->SetRegions(input->GetLargestPossibleRegion());
labelImage->SetOrigin(input->GetOrigin());
labelImage->SetSpacing(input->GetSpacing());
+ labelImage->SetDirection(input->GetDirection());
labelImage->Allocate();
labelImage->FillBuffer(m_ArgsInfo.label_arg[0]);
}
std::cout<<std::endl;
if (m_Verbose) std::cout<<"-------------"<<std::endl;
- if (m_Verbose) std::cout<<"| Label: "<<label<<" |"<<std::endl;
+ if (m_Verbose) std::cout<<"| Label: "<< (int) label<<" |"<<std::endl;
if (m_Verbose) std::cout<<"-------------"<<std::endl;
// Histograms
// Average
VFPixelType vector;
- VFPixelType zeroVector=itk::NumericTraits<VFPixelType>::Zero;
+ VFPixelType zeroVector;//=itk::NumericTraits<VFPixelType>::Zero;
+ for(unsigned int i=0;i <VFPixelType::Dimension; i++) zeroVector[i] = 0.0;
while (!(iterators[0]).IsAtEnd()) {
vector=zeroVector;
p_bar.setValue(progress);
p_bar.show();
}
-