RTK  1.4.0
Reconstruction Toolkit
rtkGgoFunctions.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright RTK Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 //#ifndef RTKGGOFUNCTIONS_H
20 //#define RTKGGOFUNCTIONS_H
21 
22 #ifndef rtkGgoFunctions_h
23 #define rtkGgoFunctions_h
24 
25 #include "rtkMacro.h"
26 #include "rtkConstantImageSource.h"
27 #include "rtkIOFactories.h"
28 #include "rtkProjectionsReader.h"
30 #include <itksys/RegularExpression.hxx>
31 
32 namespace rtk
33 {
34 
48 template< class TConstantImageSourceType, class TArgsInfo >
49 void
50 SetConstantImageSourceFromGgo(typename TConstantImageSourceType::Pointer source, const TArgsInfo &args_info)
51 {
52  typedef typename TConstantImageSourceType::OutputImageType ImageType;
53 
54  const unsigned int Dimension = ImageType::GetImageDimension();
55 
56  typename ImageType::SizeType imageDimension;
57  imageDimension.Fill(args_info.dimension_arg[0]);
58  for(unsigned int i=0; i<vnl_math_min(args_info.dimension_given, Dimension); i++)
59  imageDimension[i] = args_info.dimension_arg[i];
60 
61  typename ImageType::SpacingType imageSpacing;
62  imageSpacing.Fill(args_info.spacing_arg[0]);
63  for(unsigned int i=0; i<vnl_math_min(args_info.spacing_given, Dimension); i++)
64  imageSpacing[i] = args_info.spacing_arg[i];
65 
66  typename ImageType::PointType imageOrigin;
67  for(unsigned int i=0; i<Dimension; i++)
68  imageOrigin[i] = imageSpacing[i] * (imageDimension[i]-1) * -0.5;
69  for(unsigned int i=0; i<vnl_math_min(args_info.origin_given, Dimension); i++)
70  imageOrigin[i] = args_info.origin_arg[i];
71 
72  typename ImageType::DirectionType imageDirection;
73  if(args_info.direction_given)
74  for(unsigned int i=0; i<Dimension; i++)
75  for(unsigned int j=0; j<Dimension; j++)
76  imageDirection[i][j] = args_info.direction_arg[i*Dimension+j];
77  else
78  imageDirection.SetIdentity();
79 
80  source->SetOrigin( imageOrigin );
81  source->SetSpacing( imageSpacing );
82  source->SetDirection( imageDirection );
83  source->SetSize( imageDimension );
84  source->SetConstant( 0. );
85 
86  // Copy output image information from an existing file, if requested
87  // Overwrites parameters given in command line, if any
88  if (args_info.like_given)
89  {
90  typedef itk::ImageFileReader< ImageType > LikeReaderType;
91  typename LikeReaderType::Pointer likeReader = LikeReaderType::New();
92  likeReader->SetFileName( args_info.like_arg );
93  TRY_AND_EXIT_ON_ITK_EXCEPTION( likeReader->UpdateOutputInformation() );
94  source->SetInformationFromImage(likeReader->GetOutput());
95  }
96 
97  TRY_AND_EXIT_ON_ITK_EXCEPTION( source->UpdateOutputInformation() );
98 }
99 
115 template< class TArgsInfo >
116 const std::vector< std::string >
117 GetProjectionsFileNamesFromGgo(const TArgsInfo &args_info)
118 {
119  // Generate file names
121  names->SetDirectory(args_info.path_arg);
122  names->SetNumericSort(args_info.nsort_flag);
123  names->SetRegularExpression(args_info.regexp_arg);
124  names->SetSubMatch(args_info.submatch_arg);
125 
126  if(args_info.verbose_flag)
127  std::cout << "Regular expression matches "
128  << names->GetFileNames().size()
129  << " file(s)..."
130  << std::endl;
131 
132  // Check submatch in file names
133  if(args_info.submatch_given)
134  {
135  // Check that the submatch number returns something relevant
136  itksys::RegularExpression reg;
137  if ( !reg.compile( args_info.regexp_arg ) )
138  {
139  itkGenericExceptionMacro(<< "Error compiling regular expression " << args_info.regexp_arg);
140  }
141 
142  // Store the full filename and the selected sub expression match
143  for(size_t i=0; i<names->GetFileNames().size(); i++)
144  {
145  reg.find( names->GetFileNames()[i] );
146  if (reg.match(args_info.submatch_arg) == std::string(""))
147  {
148  itkGenericExceptionMacro(<< "Cannot find submatch " << args_info.submatch_arg
149  << " in " << names->GetFileNames()[i]
150  << " from regular expression " << args_info.regexp_arg);
151  }
152  }
153  }
154 
155  std::vector<std::string> fileNames = names->GetFileNames();
157  std::vector<size_t> idxtopop;
158  size_t i = 0;
159  for (auto fn : fileNames) {
160  itk::ImageIOBase::Pointer imageio = itk::ImageIOFactory::CreateImageIO(
161  fn.c_str(), itk::ImageIOFactory::ReadMode);
162 
163  if (imageio.IsNull()) {
164  std::cerr << "Ignoring file: " << fn << std::endl;
165  idxtopop.push_back(i);
166  }
167  i++;
168  }
169  std::reverse(idxtopop.begin(), idxtopop.end());
170  for (auto id : idxtopop) {
171  fileNames.erase(fileNames.begin() + id);
172  }
173 
174  return fileNames;
175 }
176 
177 template< class TProjectionsReaderType, class TArgsInfo >
178 void
179 SetProjectionsReaderFromGgo(typename TProjectionsReaderType::Pointer reader,
180  const TArgsInfo &args_info)
181 {
182  const std::vector< std::string > fileNames = GetProjectionsFileNamesFromGgo(args_info);
183 
184  // Vector component extraction
185  if(args_info.component_given)
186  {
187  reader->SetVectorComponent(args_info.component_arg);
188  }
189 
190  // Change image information
191  const unsigned int Dimension = TProjectionsReaderType::OutputImageType::GetImageDimension();
192  typename TProjectionsReaderType::OutputImageDirectionType direction;
193  if(args_info.newdirection_given)
194  {
195  direction.Fill(args_info.newdirection_arg[0]);
196  for(unsigned int i=0; i<args_info.newdirection_given; i++)
197  direction[i/Dimension][i%Dimension] = args_info.newdirection_arg[i];
198  reader->SetDirection(direction);
199  }
200  typename TProjectionsReaderType::OutputImageSpacingType spacing;
201  if(args_info.newspacing_given)
202  {
203  spacing.Fill(args_info.newspacing_arg[0]);
204  for(unsigned int i=0; i<args_info.newspacing_given; i++)
205  spacing[i] = args_info.newspacing_arg[i];
206  reader->SetSpacing(spacing);
207  }
208  typename TProjectionsReaderType::OutputImagePointType origin;
209  if(args_info.neworigin_given)
210  {
211  direction.Fill(args_info.neworigin_arg[0]);
212  for(unsigned int i=0; i<args_info.neworigin_given; i++)
213  origin[i] = args_info.neworigin_arg[i];
214  reader->SetOrigin(origin);
215  }
216 
217  // Crop boundaries
218  typename TProjectionsReaderType::OutputImageSizeType upperCrop, lowerCrop;
219  upperCrop.Fill(0);
220  lowerCrop.Fill(0);
221  for(unsigned int i=0; i<args_info.lowercrop_given; i++)
222  lowerCrop[i] = args_info.lowercrop_arg[i];
223  reader->SetLowerBoundaryCropSize(lowerCrop);
224  for(unsigned int i=0; i<args_info.uppercrop_given; i++)
225  upperCrop[i] = args_info.uppercrop_arg[i];
226  reader->SetUpperBoundaryCropSize(upperCrop);
227 
228  // Conditional median
229  typename TProjectionsReaderType::MedianRadiusType medianRadius;
230  medianRadius.Fill(0);
231  for(unsigned int i=0; i<args_info.radius_given; i++)
232  medianRadius[i] = args_info.radius_arg[i];
233  reader->SetMedianRadius(medianRadius);
234  if(args_info.multiplier_given)
235  reader->SetConditionalMedianThresholdMultiplier(args_info.multiplier_arg);
236 
237  // Shrink / Binning
238  typename TProjectionsReaderType::ShrinkFactorsType binFactors;
239  binFactors.Fill(1);
240  for(unsigned int i=0; i<args_info.binning_given; i++)
241  binFactors[i] = args_info.binning_arg[i];
242  reader->SetShrinkFactors(binFactors);
243 
244  // Boellaard scatter correction
245  if(args_info.spr_given)
246  reader->SetScatterToPrimaryRatio(args_info.spr_arg);
247  if(args_info.nonneg_given)
248  reader->SetNonNegativityConstraintThreshold(args_info.nonneg_arg);
249  if(args_info.airthres_given)
250  reader->SetAirThreshold(args_info.airthres_arg);
251 
252  // I0 and IDark
253  if(args_info.i0_given)
254  reader->SetI0(args_info.i0_arg);
255  reader->SetIDark(args_info.idark_arg);
256 
257  // Line integral flag
258  if(args_info.nolineint_flag)
259  reader->ComputeLineIntegralOff();
260 
261  // Water precorrection
262  if(args_info.wpc_given)
263  {
264  std::vector<double> coeffs;
265  coeffs.assign(args_info.wpc_arg, args_info.wpc_arg+args_info.wpc_given);
266  reader->SetWaterPrecorrectionCoefficients(coeffs);
267  }
268 
269  // Pass list to projections reader
270  reader->SetFileNames( fileNames );
271  TRY_AND_EXIT_ON_ITK_EXCEPTION( reader->UpdateOutputInformation() );
272 }
273 
282 template< class TArgsInfo, class TIterativeReconstructionFilter >
283 void
284 SetBackProjectionFromGgo(const TArgsInfo &args_info, TIterativeReconstructionFilter *recon)
285 {
286  switch(args_info.bp_arg)
287  {
288  case(-1): //bp__NULL
289  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_UNKNOWN);
290  break;
291  case(0): //bp_arg_VoxelBasedBackProjection
292  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_VOXELBASED);
293  break;
294  case(1): //bp_arg_Joseph
295  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPH);
296  break;
297  case(2): //bp_arg_CudaVoxelBased
298  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDAVOXELBASED);
299  break;
300  case(3): //bp_arg_CudaRayCast
301  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDARAYCAST);
302  break;
303  case(4): //bp_arg_JosephAttenuated
304  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPHATTENUATED);
305  break;
306  }
307 }
308 
317 template< class TArgsInfo, class TIterativeReconstructionFilter >
318 void
319 SetForwardProjectionFromGgo(const TArgsInfo &args_info, TIterativeReconstructionFilter *recon)
320 {
321  switch(args_info.fp_arg)
322  {
323  case(-1): //fp__NULL
324  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_UNKNOWN);
325  break;
326  case(0): //fp_arg_Joseph
327  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPH);
328  break;
329  case(1): //fp_arg_CudaRayCast
330  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDARAYCAST);
331  break;
332  case(2): //fp_arg_JosephAttenuated
333  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPHATTENUATED);
334  break;
335  }
336 }
337 
338 }
339 
340 #endif // rtkGgoFunctions_h
#define TRY_AND_EXIT_ON_ITK_EXCEPTION(execFunc)
Update a filter and catching/displaying exceptions.
Definition: rtkMacro.h:118
void SetProjectionsReaderFromGgo(typename TProjectionsReaderType::Pointer reader, const TArgsInfo &args_info)
constexpr unsigned int Dimension
void RTK_EXPORT RegisterIOFactories()
static ImageIOBasePointer CreateImageIO(const char *path, FileModeType mode)