RTK  2.5.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  * https://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  using ImageType = typename TConstantImageSourceType::OutputImageType;
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 < std::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 < std::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 < std::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({});
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  using LikeReaderType = itk::ImageFileReader<ImageType>;
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 " << names->GetFileNames().size() << " file(s)..." << std::endl;
128 
129  // Check submatch in file names
130  if (args_info.submatch_given)
131  {
132  // Check that the submatch number returns something relevant
133  itksys::RegularExpression reg;
134  if (!reg.compile(args_info.regexp_arg))
135  {
136  itkGenericExceptionMacro(<< "Error compiling regular expression " << args_info.regexp_arg);
137  }
138 
139  // Store the full filename and the selected sub expression match
140  for (const std::string & name : names->GetFileNames())
141  {
142  reg.find(name);
143  if (reg.match(args_info.submatch_arg) == std::string(""))
144  {
145  itkGenericExceptionMacro(<< "Cannot find submatch " << args_info.submatch_arg << " in " << name
146  << " from regular expression " << args_info.regexp_arg);
147  }
148  }
149  }
150 
151  std::vector<std::string> fileNames = names->GetFileNames();
153  std::vector<size_t> idxtopop;
154  size_t i = 0;
155  for (const auto & fn : fileNames)
156  {
157  itk::ImageIOBase::Pointer imageio =
159 
160  if (imageio.IsNull())
161  {
162  std::cerr << "Ignoring file: " << fn << "\n";
163  idxtopop.push_back(i);
164  }
165  i++;
166  }
167  std::reverse(idxtopop.begin(), idxtopop.end());
168  for (const auto & id : idxtopop)
169  {
170  fileNames.erase(fileNames.begin() + id);
171  }
172 
173  return fileNames;
174 }
175 
176 template <class TProjectionsReaderType, class TArgsInfo>
177 void
178 SetProjectionsReaderFromGgo(typename TProjectionsReaderType::Pointer reader, const TArgsInfo & args_info)
179 {
180  const std::vector<std::string> fileNames = GetProjectionsFileNamesFromGgo(args_info);
181 
182  // Vector component extraction
183  if (args_info.component_given)
184  {
185  reader->SetVectorComponent(args_info.component_arg);
186  }
187 
188  // Change image information
189  const unsigned int Dimension = TProjectionsReaderType::OutputImageType::GetImageDimension();
190  typename TProjectionsReaderType::OutputImageDirectionType direction;
191  if (args_info.newdirection_given)
192  {
193  direction.Fill(args_info.newdirection_arg[0]);
194  for (unsigned int i = 0; i < args_info.newdirection_given; i++)
195  direction[i / Dimension][i % Dimension] = args_info.newdirection_arg[i];
196  reader->SetDirection(direction);
197  }
198  typename TProjectionsReaderType::OutputImageSpacingType spacing;
199  if (args_info.newspacing_given)
200  {
201  spacing.Fill(args_info.newspacing_arg[0]);
202  for (unsigned int i = 0; i < args_info.newspacing_given; i++)
203  spacing[i] = args_info.newspacing_arg[i];
204  reader->SetSpacing(spacing);
205  }
206  typename TProjectionsReaderType::OutputImagePointType origin;
207  if (args_info.neworigin_given)
208  {
209  origin.Fill(args_info.neworigin_arg[0]);
210  for (unsigned int i = 0; i < args_info.neworigin_given; i++)
211  origin[i] = args_info.neworigin_arg[i];
212  reader->SetOrigin(origin);
213  }
214 
215  // Crop boundaries
216  typename TProjectionsReaderType::OutputImageSizeType upperCrop, lowerCrop;
217  upperCrop.Fill(0);
218  lowerCrop.Fill(0);
219  for (unsigned int i = 0; i < args_info.lowercrop_given; i++)
220  lowerCrop[i] = args_info.lowercrop_arg[i];
221  reader->SetLowerBoundaryCropSize(lowerCrop);
222  for (unsigned int i = 0; i < args_info.uppercrop_given; i++)
223  upperCrop[i] = args_info.uppercrop_arg[i];
224  reader->SetUpperBoundaryCropSize(upperCrop);
225 
226  // Conditional median
227  typename TProjectionsReaderType::MedianRadiusType medianRadius;
228  medianRadius.Fill(0);
229  for (unsigned int i = 0; i < args_info.radius_given; i++)
230  medianRadius[i] = args_info.radius_arg[i];
231  reader->SetMedianRadius(medianRadius);
232  if (args_info.multiplier_given)
233  reader->SetConditionalMedianThresholdMultiplier(args_info.multiplier_arg);
234 
235  // Shrink / Binning
236  typename TProjectionsReaderType::ShrinkFactorsType binFactors;
237  binFactors.Fill(1);
238  for (unsigned int i = 0; i < args_info.binning_given; i++)
239  binFactors[i] = args_info.binning_arg[i];
240  reader->SetShrinkFactors(binFactors);
241 
242  // Boellaard scatter correction
243  if (args_info.spr_given)
244  reader->SetScatterToPrimaryRatio(args_info.spr_arg);
245  if (args_info.nonneg_given)
246  reader->SetNonNegativityConstraintThreshold(args_info.nonneg_arg);
247  if (args_info.airthres_given)
248  reader->SetAirThreshold(args_info.airthres_arg);
249 
250  // I0 and IDark
251  if (args_info.i0_given)
252  reader->SetI0(args_info.i0_arg);
253  reader->SetIDark(args_info.idark_arg);
254 
255  // Line integral flag
256  if (args_info.nolineint_flag)
257  reader->ComputeLineIntegralOff();
258 
259  // Water precorrection
260  if (args_info.wpc_given)
261  {
262  std::vector<double> coeffs;
263  coeffs.assign(args_info.wpc_arg, args_info.wpc_arg + args_info.wpc_given);
264  reader->SetWaterPrecorrectionCoefficients(coeffs);
265  }
266 
267  // Pass list to projections reader
268  reader->SetFileNames(fileNames);
269  TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->UpdateOutputInformation());
270 }
271 
280 template <class TArgsInfo, class TIterativeReconstructionFilter>
281 void
282 SetBackProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFilter * recon)
283 {
284  typename TIterativeReconstructionFilter::VolumeType::Pointer attenuationMap;
285  using VolumeType = typename TIterativeReconstructionFilter::VolumeType;
286  if (args_info.attenuationmap_given)
287  {
288  // Read an existing image to initialize the attenuation map
289  attenuationMap = itk::ReadImage<VolumeType>(args_info.attenuationmap_arg);
290  }
291 
292  switch (args_info.bp_arg)
293  {
294  case (-1): // bp__NULL, keep default
295  break;
296  case (0): // bp_arg_VoxelBasedBackProjection
297  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_VOXELBASED);
298  break;
299  case (1): // bp_arg_Joseph
300  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPH);
301  break;
302  case (2): // bp_arg_CudaVoxelBased
303  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDAVOXELBASED);
304  break;
305  case (3): // bp_arg_CudaRayCast
306  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDARAYCAST);
307  if (args_info.step_given)
308  recon->SetStepSize(args_info.step_arg);
309  break;
310  case (4): // bp_arg_JosephAttenuated
311  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPHATTENUATED);
312  if (args_info.attenuationmap_given)
313  recon->SetAttenuationMap(attenuationMap);
314  break;
315  case (5): // bp_arg_RotationBased
316  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_ZENG);
317  if (args_info.sigmazero_given)
318  recon->SetSigmaZero(args_info.sigmazero_arg);
319  if (args_info.alphapsf_given)
320  recon->SetAlphaPSF(args_info.alphapsf_arg);
321  if (args_info.attenuationmap_given)
322  recon->SetAttenuationMap(attenuationMap);
323  break;
324  }
325 }
326 
335 template <class TArgsInfo, class TIterativeReconstructionFilter>
336 void
337 SetForwardProjectionFromGgo(const TArgsInfo & args_info, TIterativeReconstructionFilter * recon)
338 {
339  typename TIterativeReconstructionFilter::VolumeType::Pointer attenuationMap;
340  using VolumeType = typename TIterativeReconstructionFilter::VolumeType;
341  if (args_info.attenuationmap_given)
342  {
343  // Read an existing image to initialize the attenuation map
344  attenuationMap = itk::ReadImage<VolumeType>(args_info.attenuationmap_arg);
345  }
347  typename ClipImageType::Pointer inferiorClipImage, superiorClipImage;
348  if (args_info.inferiorclipimage_given)
349  {
350  // Read an existing image to initialize the attenuation map
351  inferiorClipImage = itk::ReadImage<ClipImageType>(args_info.inferiorclipimage_arg);
352  }
353  if (args_info.superiorclipimage_given)
354  {
355  // Read an existing image to initialize the attenuation map
356  superiorClipImage = itk::ReadImage<ClipImageType>(args_info.superiorclipimage_arg);
357  }
358 
359  switch (args_info.fp_arg)
360  {
361  case (-1): // fp__NULL, keep default
362  break;
363  case (0): // fp_arg_Joseph
364  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPH);
365  if (args_info.superiorclipimage_given)
366  recon->SetSuperiorClipImage(superiorClipImage);
367  if (args_info.inferiorclipimage_given)
368  recon->SetInferiorClipImage(inferiorClipImage);
369  break;
370  case (1): // fp_arg_CudaRayCast
371  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDARAYCAST);
372  if (args_info.step_given)
373  recon->SetStepSize(args_info.step_arg);
374  break;
375  case (2): // fp_arg_JosephAttenuated
376  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPHATTENUATED);
377  if (args_info.superiorclipimage_given)
378  recon->SetSuperiorClipImage(superiorClipImage);
379  if (args_info.inferiorclipimage_given)
380  recon->SetInferiorClipImage(inferiorClipImage);
381  if (args_info.attenuationmap_given)
382  recon->SetAttenuationMap(attenuationMap);
383  break;
384  case (3): // fp_arg_RotationBased
385  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_ZENG);
386  if (args_info.sigmazero_given)
387  recon->SetSigmaZero(args_info.sigmazero_arg);
388  if (args_info.alphapsf_given)
389  recon->SetAlphaPSF(args_info.alphapsf_arg);
390  if (args_info.attenuationmap_given)
391  recon->SetAttenuationMap(attenuationMap);
392  break;
393  }
394 }
395 
396 } // namespace rtk
397 
398 #endif // rtkGgoFunctions_h
void SetConstantImageSourceFromGgo(typename TConstantImageSourceType::Pointer source, const TArgsInfo &args_info)
Create 3D image from gengetopt specifications.
static ImageIOBasePointer CreateImageIO(const char *path, IOFileModeEnum mode)
const std::vector< std::string > GetProjectionsFileNamesFromGgo(const TArgsInfo &args_info)
Read a stack of 2D projections from gengetopt specifications.
void SetBackProjectionFromGgo(const TArgsInfo &args_info, TIterativeReconstructionFilter *recon)
Set the correct RTK backprojection type from gengetopt Given a gengetopt bp_arg option from the rtkpr...
void SetProjectionsReaderFromGgo(typename TProjectionsReaderType::Pointer reader, const TArgsInfo &args_info)
void SetForwardProjectionFromGgo(const TArgsInfo &args_info, TIterativeReconstructionFilter *recon)
Set the correct RTK forward projection type from gengetopt Given a gengetopt fp_arg option from the r...
constexpr unsigned int Dimension
void RTK_EXPORT RegisterIOFactories()
#define TRY_AND_EXIT_ON_ITK_EXCEPTION(execFunc)
Update a filter and catching/displaying exceptions.
Definition: rtkMacro.h:106