template <typename M, typename V>
void MapToVecSecond(const M & m, V & v);
+ //--------------------------------------------------------------------
+ // Find/replace string
+ template<class T>
+ int inline findAndReplace(T& source, const T& find, const T& replace);
+
#include "clitkCommon.txx"
} // end namespace
}
//--------------------------------------------------------------------
+
+//--------------------------------------------------------------------
+//http://stackoverflow.com/questions/1494399/how-do-i-search-find-and-replace-in-a-standard-string
+template<class T>
+int inline findAndReplace(T& source, const T& find, const T& replace)
+{
+ int num=0;
+ int fLen = find.size();
+ int rLen = replace.size();
+ for (int pos=0; (pos=source.find(find, pos))!=T::npos; pos+=rLen) {
+ num++;
+ source.replace(pos, fLen, replace);
+ }
+ return num;
+}
+//--------------------------------------------------------------------
+
#endif /* end #define CLITKCOMMON_TXX */
//--------------------------------------------------------------------
-void clitk::FilterBase::StartNewStep(std::string s)
+void clitk::FilterBase::StartNewStep(std::string s, bool endl)
{
if (Cancelled()) {
throw clitk::ExceptionObject("Filter is canceled.");
m_CurrentStepName = "Step "+GetCurrentStepId()+" -- "+s;
if (GetVerboseStepFlag()) {
- std::cout << m_CurrentStepName << std::endl;
+ std::cout << m_CurrentStepName;
+ if (endl) std::cout << std::endl;
//"Step " << GetCurrentStepId() << " -- " << s << std::endl;
}
}
protected:
FilterBase();
virtual ~FilterBase() {}
- void StartNewStep(std::string s);
+ void StartNewStep(std::string s, bool endl=true);
template<class TInternalImageType>
void StopCurrentStep(typename TInternalImageType::Pointer p, std::string txt="");
void StopCurrentStep();
itkGetConstMacro(Radius, double);
itkSetMacro(Radius, double);
+ itkSetMacro(K1, double);
+ itkGetMacro(K1, double);
+
typename FloatImageType::Pointer GetFuzzyMap() { return m_FuzzyMap; }
// I dont want to verify inputs information
bool m_FuzzyMapOnlyFlag;
bool m_FastFlag;
double m_Radius;
+ double m_K1;
virtual void GenerateOutputInformation();
virtual void GenerateInputRequestedRegion();
FuzzyMapOnlyFlagOff();
FastFlagOff();
SetRadius(2.0);
+ SetK1(vcl_acos(-1.0)/2);
}
//--------------------------------------------------------------------
relPosFilter->SetInput(working_image);
relPosFilter->SetAlpha1(m_Angle1[i]); // xy plane
relPosFilter->SetAlpha2(m_Angle2[i]);
- relPosFilter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2
-
- // relPosFilter->SetFast(true);
- // relPosFilter->SetRadius(1); // seems sufficient in this case
-
+ relPosFilter->SetK1(GetK1());// M_PI/2.0); // Opening parameter, default = pi/2
// relPosFilter->SetVerboseProgress(true);
+
relPosFilter->Update();
relPos = relPosFilter->GetOutput();
resampleFilter->SetGaussianFilteringEnabled(false);
resampleFilter->Update();
relPos = m_FuzzyMap = resampleFilter->GetOutput();
+ StopCurrentStep<FloatImageType>(relPos);
// Need to put exactly the same size
if (relPos->GetLargestPossibleRegion() != input->GetLargestPossibleRegion()) {
m_region.SetSize(m_size);
int start = m_index[GetDirection()];
this->SetNumberOfOutputs(m_NumberOfSlices);
+ // deprecated : use SetNumberOfIndexedInputs ? FIXME
//--------------------------------------------------------------------
// loop ExtractImageFilter with region updated, push_back
extract->SetDirectionCollapseToSubmatrix();
#endif
extract->Update();
- SetNthOutput(i, extract->GetOutput());
+ this->SetNthOutput(i, extract->GetOutput());
}
return;
}
ImagePointer image_intersection = clitk::Clone<ImageType>(input1);
clitk::Or<ImageType>(image_union, input2, GetBackgroundValue());
clitk::And<ImageType>(image_intersection, input2, GetBackgroundValue());
+
+ ImagePointer image_1NotIn2 = clitk::Clone<ImageType>(input1);
+ clitk::AndNot<ImageType>(image_1NotIn2, input2, GetBackgroundValue());
+ ImagePointer image_2NotIn1 = clitk::Clone<ImageType>(input2);
+ clitk::AndNot<ImageType>(image_2NotIn1, input1, GetBackgroundValue());
+
//writeImage<ImageType>(image_union, "union.mha");
//writeImage<ImageType>(image_intersection, "intersection.mha");
+ //writeImage<ImageType>(image_1NotIn2, "image_1NotIn2.mha");
+ //writeImage<ImageType>(image_2NotIn1, "image_2NotIn1.mha");
// Compute size
typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
statFilter->Update();
int in2 = statFilter->GetCount(GetLabel1());
+ statFilter->SetInput(image_1NotIn2);
+ statFilter->SetLabelInput(image_1NotIn2);
+ statFilter->Update();
+ int l1notIn2 = statFilter->GetCount(GetLabel1());
+
+ statFilter->SetInput(image_2NotIn1);
+ statFilter->SetLabelInput(image_2NotIn1);
+ statFilter->Update();
+ int l2notIn1 = statFilter->GetCount(GetLabel1());
+
double dice = 2.0*(double)inter/(double)(in1+in2);
int width = 6;
std::cout << std::setw(width) << in1 << " "
<< std::setw(width) << in2 << " "
<< std::setw(width) << inter << " "
<< std::setw(width) << u << " "
+ << std::setw(width) << l1notIn2 << " "
+ << std::setw(width) << l2notIn1 << " "
<< std::setw(width) << dice << " "; //std::endl;
}
//--------------------------------------------------------------------
template<class ImageType>
typename ImageType::Pointer
AutoCrop(const ImageType * input,
- typename ImageType::PixelType BG) {
+ typename ImageType::PixelType BG,
+ const bool useBorderFlag=false) {
typedef clitk::AutoCropFilter<ImageType> AutoCropFilterType;
typename AutoCropFilterType::Pointer autoCropFilter = AutoCropFilterType::New();
autoCropFilter->SetInput(input);
autoCropFilter->SetBackgroundValue(BG);
+ autoCropFilter->SetUseBorder(useBorderFlag);
autoCropFilter->Update();
return autoCropFilter->GetOutput();
}
++iter;
}
if (!found) return false;
- input->TransformIndexToPhysicalPoint(max, point);
+ input->TransformIndexToPhysicalPoint(max, point); // half of the pixel
return true;
}
//--------------------------------------------------------------------
int dim, double max, bool autoCrop,
typename ImageType::PixelType BG)
{
- typename ImageType::PointType p;
+ typename ImageType::PointType p;
+
image->TransformIndexToPhysicalPoint(image->GetLargestPossibleRegion().GetIndex()+
image->GetLargestPossibleRegion().GetSize(), p);
- // Add GetSpacing because remove Lower or equal than
- // DD(max);
- // DD(p);
- // DD(max+image->GetSpacing()[dim]);
- return CropImageAlongOneAxis<ImageType>(image, dim, max+image->GetSpacing()[dim], p[dim], autoCrop, BG);
+
+ return CropImageAlongOneAxis<ImageType>(image, dim, max, p[dim], autoCrop, BG);
}
//--------------------------------------------------------------------
// Compute region size
typename ImageType::RegionType region;
typename ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
- typename ImageType::PointType p = image->GetOrigin();
+
+ // Starting index
+ typename ImageType::PointType p = image->GetOrigin(); // not at pixel center !
if (min > p[dim]) p[dim] = min; // Check if not outside the image
typename ImageType::IndexType start;
image->TransformPhysicalPointToIndex(p, start);
- double m = image->GetOrigin()[dim] + size[dim]*image->GetSpacing()[dim];
+
+ // Size of the region
+ // -1 because last point is size -1
+ double m = image->GetOrigin()[dim] + (size[dim]-1)*image->GetSpacing()[dim];
if (max > m) p[dim] = m; // Check if not outside the image
else p[dim] = max;
+
typename ImageType::IndexType end;
image->TransformPhysicalPointToIndex(p, end);
- size[dim] = abs(end[dim]-start[dim]);
+ size[dim] = abs(end[dim]-start[dim])+1;// +1 because we want to include the point.
+
+ // Set region
region.SetIndex(start);
region.SetSize(size);
#include "clitkSliceBySliceRelativePositionFilter.txx"
#endif
+typedef unsigned char PixelType_uchar;
+typedef itk::Image<PixelType_uchar, 3> ImageType_uchar;
+extern template class clitk::SliceBySliceRelativePositionFilter<ImageType_uchar>;
+
#endif
SetObjectCCLSelectionDirection(1);
ObjectCCLSelectionIgnoreSingleCCLFlagOff();
VerboseSlicesFlagOff();
+ this->SetK1(vcl_acos(-1.0)/2);
}
//--------------------------------------------------------------------
<< "ObjectCCLSelectionIgnoreSingleCCLFlag = " << this->GetObjectCCLSelectionIgnoreSingleCCLFlag() << std::endl
<< "IgnoreEmptySliceObjectFlag = " << this->GetIgnoreEmptySliceObjectFlag() << std::endl
<< "(RP) FastFlag = " << this->GetFastFlag() << std::endl
- << "(RP) Radius = " << this->GetRadius() << std::endl;
+ << "(RP) Radius = " << this->GetRadius() << std::endl
+ << "(RP) K1 = " << this->GetK1() << std::endl;
}
//--------------------------------------------------------------------
int nb=0;
mObjectSlices[i] = LabelizeAndCountNumberOfObjects<SliceType>(mObjectSlices[i], 0, true, 1, nb);
+ if (GetVerboseSlicesFlag()) {
+ std::cout << "slice " << i << " nb = " << nb << std::endl;
+ }
+
// If no object and empty slices and if we need the full fuzzy map, create a dummy one.
if ((nb==0) && (this->GetFuzzyMapOnlyFlag())) {
typename FloatSliceType::Pointer one = FloatSliceType::New();
if (GetVerboseSlicesFlag()) {
std::cout << "Slice " << i << std::endl;
relPosFilter->VerboseStepFlagOn();
+ //relPosFilter->WriteStepFlagOn();
}
relPosFilter->WriteStepFlagOff();
// relPosFilter->VerboseMemoryFlagOn();
// relPosFilter->SetComputeFuzzyMapFlag(this->GetComputeFuzzyMapFlag());
relPosFilter->SetFastFlag(this->GetFastFlag());
relPosFilter->SetRadius(this->GetRadius());
+ relPosFilter->SetK1(this->GetK1());
// Go !
relPosFilter->Update();
ADD_LIBRARY(clitkSegmentationGgoLib
clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx
clitkAnatomicalFeatureDatabase.cxx
+ clitkSliceBySliceRelativePositionFilter_uchar.cxx
${GGO_C_FILES})
#=========================================================
{
ADD_IMAGE_TYPE(Dim, uchar);
ADD_IMAGE_TYPE(Dim, short);
+ ADD_IMAGE_TYPE(Dim, ushort);
// ADD_IMAGE_TYPE(Dim, int);
ADD_IMAGE_TYPE(Dim, float);
}
for(uint i=0; i<mArgsInfoList.size(); i++) {
// clitk::PrintMemory(true, "Start");
- std::string text = "["+s+"] limits ";
+ // remove _S in station name
+ std::string sname = s;
+ clitk::findAndReplace<std::string>(sname, "_S", " ");
+ std::string text = "["+sname+"] ";
if (mArgsInfoList[i].orientation_given) text += std::string(mArgsInfoList[i].orientation_arg[0])+" ";
else text = text+"("+toString(mArgsInfoList[i].angle1_arg)+" "+
toString(mArgsInfoList[i].angle2_arg)+" "+
(mArgsInfoList[i].inverse_flag?"true":"false")+") ";
text = text+mArgsInfoList[i].object_arg+" "+toString(mArgsInfoList[i].threshold_arg);
if (mArgsInfoList[i].sliceBySlice_flag) {
- text += " slice by slice";
+ text += " SbS";
}
else text += " 3D";
- text += " spacing=" + toString(mArgsInfoList[i].spacing_arg);
+ text += " sp=" + toString(mArgsInfoList[i].spacing_arg)+" ";
- StartNewStep(text);
+ StartNewStep(text, false); // no endl
typename RelPosFilterType::Pointer relPosFilter;
// Is it slice by slice or 3D ?
f->SetUniqueConnectedComponentBySliceFlag(mArgsInfoList[i].uniqueCCL_flag);
f->SetObjectCCLSelectionFlag(mArgsInfoList[i].uniqueObjectCCL_flag);
f->IgnoreEmptySliceObjectFlagOn();
+ f->SetVerboseSlicesFlag(mArgsInfoList[i].verboseSlices_flag);
//f->SetObjectCCLSelectionDimension(0);
//f->SetObjectCCLSelectionDirection(-1);
//f->SetAutoCropFlag(false);
filter->SetInput(1, m_reference);
filter->Update();
}
-
+ std::cout << std::endl;
}
}
//--------------------------------------------------------------------
ImagePointer object = GetAFDB()->template GetImage<ImageType>(options.object_arg);
filter->SetInputObject(object);
filter->WriteStepFlagOff();
+ if (options.writeStep_flag) filter->WriteStepFlagOn();
filter->SetVerboseImageSizeFlag(GetVerboseImageSizeFlag());
filter->SetFuzzyThreshold(options.threshold_arg);
filter->SetInverseOrientationFlag(options.inverse_flag); // MUST BE BEFORE AddOrientationTypeString
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+ ======================================================================-====*/
+
+#ifndef CLITKSLICEBYSLICERELATIVEPOSITIONFILTER_UCHARH
+#define CLITKSLICEBYSLICERELATIVEPOSITIONFILTER_UCHARH
+
+// clitk
+#include "clitkSliceBySliceRelativePositionFilter.h"
+
+typedef unsigned char PixelType_uchar;
+typedef itk::Image<PixelType_uchar, 3> ImageType_uchar;
+template class clitk::SliceBySliceRelativePositionFilter<ImageType_uchar>;
+
+#endif
typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
typename itk::Vector<double,Dimension> invTrans = clitk::GetTranslationPartMatrix(invMatrix);
+ // Display warning
+ if (m_ArgsInfo.spacing_given)
+ std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl;
+ if (m_ArgsInfo.origin_given)
+ std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
+
// Spacing is influenced by affine transform matrix and input direction
typename InputImageType::SpacingType outputSpacing;
outputSpacing = invRotMatrix *
template<unsigned int Dim>
void clitk::AutoCropGenericFilter<ArgsInfoType>::InitializeImageType()
{
- // ADD_DEFAULT_IMAGE_TYPES(Dim);
+ //ADD_DEFAULT_IMAGE_TYPES(Dim);
ADD_IMAGE_TYPE(Dim, uchar);
ADD_IMAGE_TYPE(Dim, ushort);
+ ADD_IMAGE_TYPE(Dim, short);
// ADD_IMAGE_TYPE(Dim, uint);
// ADD_IMAGE_TYPE(Dim, ulong);
// ADD_IMAGE_TYPE(Dim, int);
===========================================================================**/
#ifndef CLITKIMAGEARITHM_CXX
#define CLITKIMAGEARITHM_CXX
-/**
- -------------------------------------------------
- * @file clitkImageArithm.cxx
- * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
- * @date 23 Feb 2008 08:37:53
- -------------------------------------------------*/
// clitk include
#include "clitkImageArithm_ggo.h"
option "output" o "Output image filename" string yes
option "scalar" s "Scalar value" double no
-option "operation" t "Type of operation : \n With another image : 0=add*, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff, 7=difference*, 8=relativ diff\n; For 'scalar' : 0=add*, 1=multiply*, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide* 12=normalize (divide by max); \n* operations supported with vector fields as inputs." int default="0" no
+option "operation" t "Type of operation : \n With another image : 0=add*, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff, 7=difference*, 8=relativ diff\n; For 'scalar' : 0=add*, 1=multiply*, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide* 12=normalize (divide by max) 13=-ln(I/IO)**; \n* operations supported with vector fields as inputs. \n** for fluence image, if pixel value == 0, consider value=0.5" int default="0" no
option "pixelValue" - "Default value for NaN/Inf" double default="0.0" no
option "setFloatOutput" f "Set output to float pixel type" flag off
++ito;
}
break;
+ case 13: // -ln I/I0
+ while (!it.IsAtEnd()) {
+ if (it.Get() == 0) { // special case for fluence image with 0 value in a pixel -> consider 0.5
+ ito.Set(-log(0.5 / mScalar) );
+ }
+ else {
+ ito.Set(-log(PixelTypeDownCast<double, PixelType>((double)it.Get() / mScalar)) );
+ }
+ ++it;
+ ++ito;
+ }
+ break;
default: // error ?
std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
exit(-1);
if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
UpdateWithDimAndPixelType<Dimension, unsigned char, Components>();
}
-
-// else if (PixelType == "unsigned_int"){
-// if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_int..." << std::endl;
-// UpdateWithDimAndPixelType<Dimension, unsigned int, Components>();
-// }
- // else if (PixelType == "char"){
- // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
- // UpdateWithDimAndPixelType<Dimension, signed char>();
- // }
+
else if(PixelType == "double"){
if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
UpdateWithDimAndPixelType<Dimension, double, Components>();
typename ResamplerType::Pointer resampler = ResamplerType::New();
resampler->SetInput(labelImage);
resampler->SetOutputSpacing(input->GetSpacing());
+ resampler->SetGaussianFilteringEnabled(false);
resampler->Update();
labelImage = resampler->GetOutput();
+ //writeImage<LabelImageType>(labelImage, "test1.mha");
+
+ typedef clitk::CropLikeImageFilter<LabelImageType> FilterType;
+ typename FilterType::Pointer crop = FilterType::New();
+ crop->SetInput(labelImage);
+ crop->SetCropLikeImage(input);
+ crop->Update();
+ labelImage = crop->GetOutput();
+ //writeImage<LabelImageType>(labelImage, "test2.mha");
- typename itk::ImageBase<LabelImageType::ImageDimension>::RegionType reg
- = input->GetLargestPossibleRegion();
- labelImage = ResizeImageLike<LabelImageType>(labelImage, ®, 0);
}
else {
std::cerr << "Mask image has a different size/spacing than input. Abort" << std::endl;
//----------------------------------------------------------------------------
void vvBinaryImageOverlayActor::ComputeExtent(int * inExtent, int * outExtent, vtkImageData * image, vtkImageData * overlay)
{
+ for(int i=0; i<3; i++) {
+ double a = (image->GetOrigin()[i] + inExtent[i*2]*image->GetSpacing()[i] -
+ overlay->GetOrigin()[i]) / overlay->GetSpacing()[i];
+ double b = (image->GetOrigin()[i] + inExtent[i*2+1]*image->GetSpacing()[i] -
+ overlay->GetOrigin()[i]) / overlay->GetSpacing()[i];
+ outExtent[i*2] = lrint(a);
+ outExtent[i*2+1] = lrint(b);
+ }
+
+ /* // FIXME (original)
outExtent[0] = (int)lrint(((image->GetOrigin()[0] + inExtent[0]*image->GetSpacing()[0]) - overlay->GetOrigin()[0]) / overlay->GetSpacing()[0]);
outExtent[1] = (int)lrint(((image->GetOrigin()[0] + inExtent[1]*image->GetSpacing()[0]) - overlay->GetOrigin()[0]) / overlay->GetSpacing()[0]);
outExtent[2] = (int)lrint(((image->GetOrigin()[1] + inExtent[2]*image->GetSpacing()[1]) - overlay->GetOrigin()[1]) / overlay->GetSpacing()[1]);
outExtent[3] = (int)lrint(((image->GetOrigin()[1] + inExtent[3]*image->GetSpacing()[1]) - overlay->GetOrigin()[1]) / overlay->GetSpacing()[1]);
outExtent[4] = (int)lrint(((image->GetOrigin()[2] + inExtent[4]*image->GetSpacing()[2]) - overlay->GetOrigin()[2]) / overlay->GetSpacing()[2]);
outExtent[5] = (int)lrint(((image->GetOrigin()[2] + inExtent[5]*image->GetSpacing()[2]) - overlay->GetOrigin()[2]) / overlay->GetSpacing()[2]);
+ */
}
//----------------------------------------------------------------------------