RTK  2.1.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  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(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  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  switch (args_info.bp_arg)
285  {
286  case (-1): // bp__NULL, keep default
287  break;
288  case (0): // bp_arg_VoxelBasedBackProjection
289  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_VOXELBASED);
290  break;
291  case (1): // bp_arg_Joseph
292  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPH);
293  break;
294  case (2): // bp_arg_CudaVoxelBased
295  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDAVOXELBASED);
296  break;
297  case (3): // bp_arg_CudaRayCast
298  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_CUDARAYCAST);
299  break;
300  case (4): // bp_arg_JosephAttenuated
301  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_JOSEPHATTENUATED);
302  break;
303  case (5): // bp_arg_RotationBased
304  recon->SetBackProjectionFilter(TIterativeReconstructionFilter::BP_ZENG);
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, keep default
324  break;
325  case (0): // fp_arg_Joseph
326  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPH);
327  break;
328  case (1): // fp_arg_CudaRayCast
329  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_CUDARAYCAST);
330  break;
331  case (2): // fp_arg_JosephAttenuated
332  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_JOSEPHATTENUATED);
333  break;
334  case (3): // fp_arg_RotationBased
335  recon->SetForwardProjectionFilter(TIterativeReconstructionFilter::FP_ZENG);
336  break;
337  }
338 }
339 
340 } // namespace rtk
341 
342 #endif // rtkGgoFunctions_h
static ImageIOBasePointer CreateImageIO(const char *path, IOFileModeEnum mode)
void SetProjectionsReaderFromGgo(typename TProjectionsReaderType::Pointer reader, const TArgsInfo &args_info)
constexpr unsigned int Dimension
void RTK_EXPORT RegisterIOFactories()