1 #include <tclap/CmdLine.h> 2 #include <itkN4BiasFieldCorrectionImageFilter.h> 4 #include <itkConstantPadImageFilter.h> 5 #include <itkOtsuThresholdImageFilter.h> 6 #include <itkShrinkImageFilter.h> 7 #include <itkCommand.h> 8 #include <itkTimeProbe.h> 11 void eventCallback(itk::Object* caller,
const itk::EventObject& event,
void* clientData)
13 itk::ProcessObject * processObject = (itk::ProcessObject*) caller;
14 std::cout <<
"\033[K\rProgression: " << (int)(processObject->GetProgress() * 100) <<
"%" << std::flush;
22 std::istringstream(pi_rsValue)>> numRes;
30 std::vector<T> oValuesVectorRes;
32 size_t iSeparatorPosition = pi_rsVector.find(
'x', 0);
34 if (iSeparatorPosition == std::string::npos)
36 oValuesVectorRes.push_back(convertStringToNumber<T>(pi_rsVector));
40 std::string sChunk = pi_rsVector.substr(0, iSeparatorPosition);
42 oValuesVectorRes.push_back(convertStringToNumber<T>(sChunk));
43 while (iSeparatorPosition != std::string::npos)
45 std::string::size_type crossposfrom = iSeparatorPosition;
46 iSeparatorPosition = pi_rsVector.find(
'x', crossposfrom + 1);
47 if (iSeparatorPosition == std::string::npos)
49 sChunk = pi_rsVector.substr(crossposfrom + 1, pi_rsVector.length());
53 sChunk = pi_rsVector.substr(crossposfrom + 1, iSeparatorPosition);
55 oValuesVectorRes.push_back(convertStringToNumber<unsigned int>(sChunk));
59 return oValuesVectorRes;
64 char const *pchTemp = sIterations.c_str();
65 for (
int i = 0; i < sIterations.length(); ++i)
67 if (!((pchTemp[i] >=
'0' && pchTemp[i] <=
'9') || pchTemp[i] ==
'x'))
69 std::cerr <<
"error: " <<
" for argument iterations does not respect format strictly positive number separated by 'x' like 100x50x50..." << std::endl;
78 return ConvertVector<unsigned int>(pi_rsIterations);
81 int main(
int argc,
char *argv[])
83 TCLAP::CmdLine cmd(
"INRIA / IRISA - VisAGeS/Empenn Team",
' ', ANIMA_VERSION);
84 TCLAP::ValueArg<std::string> oArgInputImg(
"i",
"input",
"Input image.",
true,
"",
"input Name", cmd);
85 TCLAP::ValueArg<std::string> oArgOutputName(
"o",
"output",
"Name for output file",
true,
"",
"output Name", cmd);
86 TCLAP::ValueArg<std::string> oArgIterationsString(
"I",
"iterations",
"Table of number of iterations, default=50x40x30",
false,
"50x40x30",
"iterations table", cmd);
87 TCLAP::ValueArg<int> oArgShrinkFactors(
"S",
"shrinkFactor",
"Shrink factor, default=4",
false, 4,
"shrink Factor", cmd);
88 TCLAP::ValueArg<double> oArgWienerFilterNoise(
"W",
"wiener",
"Wiener Filter Noise, default=0.01",
false, 0.01,
"wiener noise", cmd);
89 TCLAP::ValueArg<double> oArgbfFWHM(
"B",
"bfFWHM",
"Bias field Full Width at Half Maximum, default=0.15",
false, 0.15,
"Bias field Full Width at Half Maximum", cmd);
90 TCLAP::ValueArg<double> oArgConvergenceThreshold(
"T",
"threshold",
"Convergence Threshold, default=0.0001",
false, 0.0001,
"threshold", cmd);
91 TCLAP::ValueArg<int> oArgSplineOrder(
"O",
"splineOrder",
"BSpline Order, default=3",
false, 3,
"spline Order", cmd);
92 TCLAP::ValueArg<double> oArgSplineDistance(
"D",
"splineDistance",
"B-Spline distance, default=0.0",
false, 0.0,
"spline Distance", cmd);
93 TCLAP::ValueArg<std::string> oArgInitialMeshResolutionString(
"G",
"splineGrid",
"B-Spline grid resolution. It is ignored if splineDistance>0 or if dimention of it <>3, default=1x1x1",
false,
"1x1x1",
"spline Grid", cmd);
94 TCLAP::ValueArg<unsigned int> oArgNbP(
"p",
"numberofthreads",
"Number of threads to run on (default: all cores)",
false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(),
"number of threads", cmd);
98 cmd.parse( argc, argv );
100 catch (TCLAP::ArgException & e)
102 std::cerr <<
"error: " << e.error() <<
" for arg " << e.argId() << std::endl;
106 typedef itk::Image<double, 3 > ImageType;
107 typedef itk::Image<unsigned char, 3> MaskImageType;
108 typedef itk::N4BiasFieldCorrectionImageFilter<ImageType, MaskImageType, ImageType> BiasFilter;
109 typedef itk::ConstantPadImageFilter<ImageType, ImageType> PadderType;
110 typedef itk::ConstantPadImageFilter<MaskImageType, MaskImageType> MaskPadderType;
111 typedef itk::ShrinkImageFilter<ImageType, ImageType> ShrinkerType;
112 typedef itk::ShrinkImageFilter<MaskImageType, MaskImageType> MaskShrinkerType;
113 typedef itk::BSplineControlPointImageFilter< BiasFilter::BiasFieldControlPointLatticeType, BiasFilter::ScalarImageType> BSplinerType;
114 typedef itk::ExpImageFilter<ImageType, ImageType> ExpFilterType;
115 typedef itk::DivideImageFilter<ImageType, ImageType, ImageType> DividerType;
116 typedef itk::ExtractImageFilter<ImageType, ImageType> CropperType;
119 std::vector<double> oInitialMeshResolutionVect = ConvertVector<double>(oArgInitialMeshResolutionString.getValue());
120 double fSplineDistance = oArgSplineDistance.getValue();
128 BiasFilter::Pointer filter = BiasFilter::New();
129 BiasFilter::ArrayType oNumberOfControlPointsArray;
132 ImageType::Pointer image = anima::readImage<ImageType>( oArgInputImg.getValue());
135 std::cout <<
"Creating Otsu mask." << std::endl;
136 itk::TimeProbe timer;
138 MaskImageType::Pointer maskImage = ITK_NULLPTR;
139 typedef itk::OtsuThresholdImageFilter<ImageType, MaskImageType> ThresholderType;
140 ThresholderType::Pointer otsu = ThresholderType::New();
141 otsu->SetInput(image);
142 otsu->SetNumberOfHistogramBins(200);
143 otsu->SetInsideValue(0);
144 otsu->SetOutsideValue(1);
146 otsu->SetNumberOfWorkUnits(oArgNbP.getValue());
148 maskImage = otsu->GetOutput();
151 BiasFilter::VariableSizeArrayType itkTabMaximumIterations;
152 itkTabMaximumIterations.SetSize(oMaxNumbersIterationsVector.size());
153 for (
int i=0; i<oMaxNumbersIterationsVector.size(); ++i)
155 itkTabMaximumIterations[i] = oMaxNumbersIterationsVector[i];
157 filter->SetMaximumNumberOfIterations(itkTabMaximumIterations);
160 BiasFilter::ArrayType oFittingLevelsTab;
161 oFittingLevelsTab.Fill(oMaxNumbersIterationsVector.size());
162 filter->SetNumberOfFittingLevels(oFittingLevelsTab);
165 ImageType::IndexType oImageIndex = image->GetLargestPossibleRegion().GetIndex();
166 ImageType::SizeType oImageSize = image->GetLargestPossibleRegion().GetSize();
167 ImageType::PointType newOrigin = image->GetOrigin();
169 if (fSplineDistance>0)
172 itk::SizeValueType lowerBound[3];
173 itk::SizeValueType upperBound[3];
175 for (
unsigned int i = 0; i < 3; i++)
177 double domain =
static_cast<double>(image->GetLargestPossibleRegion().GetSize()[i] - 1) * image->GetSpacing()[i];
178 unsigned int numberOfSpans =
static_cast<unsigned int>(std::ceil(domain / fSplineDistance));
179 unsigned long extraPadding =
static_cast<unsigned long>((numberOfSpans * fSplineDistance - domain) / image->GetSpacing()[i] + 0.5);
180 lowerBound[i] =
static_cast<unsigned long>(0.5 * extraPadding);
181 upperBound[i] = extraPadding - lowerBound[i];
182 newOrigin[i] -= (
static_cast<double>(lowerBound[i]) * image->GetSpacing()[i]);
183 oNumberOfControlPointsArray[i] = numberOfSpans + filter->GetSplineOrder();
187 PadderType::Pointer imagePadder = PadderType::New();
188 imagePadder->SetInput(image);
189 imagePadder->SetPadLowerBound(lowerBound);
190 imagePadder->SetPadUpperBound(upperBound);
191 imagePadder->SetConstant(0);
192 imagePadder->SetNumberOfWorkUnits(oArgNbP.getValue());
193 imagePadder->Update();
195 image = imagePadder->GetOutput();
198 MaskPadderType::Pointer maskPadder = MaskPadderType::New();
199 maskPadder->SetInput(maskImage);
200 maskPadder->SetPadLowerBound(lowerBound);
201 maskPadder->SetPadUpperBound(upperBound);
202 maskPadder->SetConstant(0);
204 maskPadder->SetNumberOfWorkUnits(oArgNbP.getValue());
205 maskPadder->Update();
207 maskImage = maskPadder->GetOutput();
210 filter->SetNumberOfControlPoints(oNumberOfControlPointsArray);
212 else if(oInitialMeshResolutionVect.size() == 3)
215 for (
unsigned i = 0; i < 3; i++)
217 oNumberOfControlPointsArray[i] =
static_cast<unsigned int>(oInitialMeshResolutionVect[i]) + filter->GetSplineOrder();
219 filter->SetNumberOfControlPoints(oNumberOfControlPointsArray);
223 std::cout <<
"No BSpline distance and Mesh Resolution is ignored because not 3 dimensions" << std::endl;
227 ShrinkerType::Pointer imageShrinker = ShrinkerType::New();
228 imageShrinker->SetInput(image);
229 imageShrinker->SetShrinkFactors(1);
232 MaskShrinkerType::Pointer maskShrinker = MaskShrinkerType::New();
233 maskShrinker->SetInput(maskImage);
234 maskShrinker->SetShrinkFactors(1);
237 imageShrinker->SetShrinkFactors(oArgShrinkFactors.getValue());
238 maskShrinker->SetShrinkFactors(oArgShrinkFactors.getValue());
239 imageShrinker->SetNumberOfWorkUnits(oArgNbP.getValue());
240 maskShrinker->SetNumberOfWorkUnits(oArgNbP.getValue());
241 imageShrinker->Update();
242 maskShrinker->Update();
245 filter->SetSplineOrder(oArgSplineOrder.getValue());
246 filter->SetWienerFilterNoise(oArgWienerFilterNoise.getValue());
247 filter->SetBiasFieldFullWidthAtHalfMaximum(oArgbfFWHM.getValue());
248 filter->SetConvergenceThreshold(oArgConvergenceThreshold.getValue());
249 filter->SetInput(imageShrinker->GetOutput());
250 filter->SetMaskImage(maskShrinker->GetOutput());
253 itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New();
255 filter->AddObserver(itk::ProgressEvent(), callback);
258 filter->SetNumberOfWorkUnits(oArgNbP.getValue());
261 catch (itk::ExceptionObject & err)
263 std::cerr <<
"ExceptionObject caught !" << std::endl;
264 std::cerr << err << std::endl;
273 BSplinerType::Pointer bspliner = BSplinerType::New();
274 bspliner->SetInput(filter->GetLogBiasFieldControlPointLattice());
275 bspliner->SetSplineOrder(filter->GetSplineOrder());
276 bspliner->SetSize(image->GetLargestPossibleRegion().GetSize());
277 bspliner->SetOrigin(newOrigin);
278 bspliner->SetDirection(image->GetDirection());
279 bspliner->SetSpacing(image->GetSpacing());
280 bspliner->SetNumberOfWorkUnits(oArgNbP.getValue());
283 ImageType::Pointer logField = ImageType::New();
284 logField->SetOrigin(image->GetOrigin());
285 logField->SetSpacing(image->GetSpacing());
286 logField->SetRegions(image->GetLargestPossibleRegion());
287 logField->SetDirection(image->GetDirection());
288 logField->Allocate();
290 itk::ImageRegionIterator<BiasFilter::ScalarImageType> IB(bspliner->GetOutput(), bspliner->GetOutput()->GetLargestPossibleRegion());
292 itk::ImageRegionIterator<ImageType> IF(logField, logField->GetLargestPossibleRegion());
294 for (IB.GoToBegin(), IF.GoToBegin(); !IB.IsAtEnd(); ++IB, ++IF)
299 ExpFilterType::Pointer expFilter = ExpFilterType::New();
300 expFilter->SetInput(logField);
301 expFilter->SetNumberOfWorkUnits(oArgNbP.getValue());
304 DividerType::Pointer divider = DividerType::New();
305 divider->SetInput1(image);
306 divider->SetInput2(expFilter->GetOutput());
307 divider->SetNumberOfWorkUnits(oArgNbP.getValue());
310 ImageType::RegionType inputRegion;
311 inputRegion.SetIndex(oImageIndex);
312 inputRegion.SetSize(oImageSize);
314 CropperType::Pointer cropper = CropperType::New();
315 cropper->SetInput(divider->GetOutput());
316 cropper->SetExtractionRegion(inputRegion);
317 cropper->SetDirectionCollapseToSubmatrix();
318 cropper->SetNumberOfWorkUnits(oArgNbP.getValue());
322 std::cout <<
"\nComputation time : " << timer.GetTotal() << std::endl;
327 anima::writeImage<ImageType>(oArgOutputName.getValue(), cropper->GetOutput());
329 catch (itk::ExceptionObject & err)
331 std::cerr <<
"ExceptionObject caught !" << std::endl;
332 std::cerr << err << std::endl;
int main(int argc, char *argv[])
void checkIterationArg(std::string &sIterations)
T convertStringToNumber(std::string const &pi_rsValue)
std::vector< unsigned int > extractMaxNumberOfIterationsVector(std::string pi_rsIterations)
std::vector< T > ConvertVector(std::string &pi_rsVector)
void eventCallback(itk::Object *caller, const itk::EventObject &event, void *clientData)