diff --git a/.clang-format b/.clang-format
index 1606cc0cd..19d80912e 100644
--- a/.clang-format
+++ b/.clang-format
@@ -25,6 +25,8 @@ BraceWrapping:
SplitEmptyFunction: false
SplitEmptyRecord: true
SplitEmptyNamespace: true
+BinPackParameters: false
+BinPackArguments: false
ColumnLimit: 130
IndentPPDirectives: AfterHash
PointerAlignment: Left
@@ -32,5 +34,6 @@ SortIncludes: false
SortUsingDeclarations: false
SpaceBeforeParens: ControlStatements
PackConstructorInitializers: Never
+AllowShortFunctionsOnASingleLine: InlineOnly
Standard: Cpp11
...
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index fe361a407..3b29a13e0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,5 +1,5 @@
repos:
-- repo: git://github.com/doublify/pre-commit-clang-format
+- repo: https://github.com/doublify/pre-commit-clang-format
rev: '62302476'
hooks:
- id: clang-format
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 3c8ea7411..c8ce45824 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -29,23 +29,25 @@ This is our recommended process. If it sounds too daunting, ask for help.
3. Create a branch in your fork with a descriptive name and put your fixes there. If your fix is
simple you could do it on github by editing a file, otherwise clone your project (or add a remote
to your current git clone) and work as usual.
-4. If your change is important, add it to the release notesfor the upcoming version, [see](https://github.com/UCL/STIR/blob/master/documentation/)
+4. Configure your editor and potentially even [pre-commit](https://pre-commit.com/), see
+[documentation/devel/README.md](documentation/devel/README.md).
+5. If your change is important, add it to the release notes for the upcoming version in the [documentation folder](https://github.com/UCL/STIR/tree/master/documentation/)
and even the [User's Guide](https://github.com/UCL/STIR/blob/master/documentation/STIR-UsersGuide.tex) or other documentation files.
-5. Use [well-formed commit messages](http://tbaggery.com/2008/04/19/a-note-about-git-commit-messages.html)
+6. Use [well-formed commit messages](http://tbaggery.com/2008/04/19/a-note-about-git-commit-messages.html)
for each change (in particular with a single "subject" line
-followed by an empty line and then more details). If the change affects comments only, it is recommended to put `[ci skip]` in your subject line. This avoids unnecessary computation, and clogging our Travis/Appveyor queues.
-6. Push the commits to your fork and submit a [pull request (PR)](https://help.github.com/articles/creating-a-pull-request)
-(enable changes by project admins.) Give your pull request a descriptive name (i.e. don't call if *Fix #issuenumber*. Be prepared to add further commits to your branch after discussion.
-In the description of the PR, add a statement about which Issue this applies to
-using [a phrase such that github auto-closes the issue when merged to master](https://help.github.com/articles/closing-issues-using-keywords/).
-7. Be prepared to add further commits to your branch after discussion.
+followed by an empty line and then more details).
Please by mindful about the resources used by our Continuous Integration (CI) workflows:
- Group your commits and only push once your code compiles and tests succeed on your machine (ideally you have sensible commit messages at every stage)
- Use specific keywords in the first line of the last commit that you push to prevent CI being run:
- `[ci skip]` skips all CI runs (e.g. when you only change documentation, or when your update isn't ready yet)
- `[actions skip]` does not run GitHub Actions, see [here](https://github.blog/changelog/2021-02-08-github-actions-skip-pull-request-and-push-workflows-with-skip-ci/). Note: this can be in the main commit message.
- `[skip appveyor]` does not run Appveyor, see [here](https://www.appveyor.com/docs/how-to/filtering-commits/#skip-directive-in-commit-message)
-8. After acceptance of your PR, go home with a nice warm feeling.
+7. Push the commits to your fork and submit a [pull request (PR)](https://help.github.com/articles/creating-a-pull-request)
+(enable changes by project admins.) Give your pull request a descriptive name (i.e. don't call if *Fix #issuenumber*. Be prepared to add further commits to your branch after discussion.
+In the description of the PR, add a statement about which Issue this applies to
+using [a phrase such that github auto-closes the issue when merged to master](https://help.github.com/articles/closing-issues-using-keywords/).
+8. Be prepared to add further commits to your branch after discussion.
+9. After acceptance of your PR, go home with a nice warm feeling.
Suggested reading:
https://help.github.com/articles/fork-a-repo/, https://git-scm.com/book/en/v2/GitHub-Contributing-to-a-Project or https://guides.github.com/activities/forking/.
diff --git a/documentation/STIR-developers-overview.tex b/documentation/STIR-developers-overview.tex
index 08f4d7121..5746a292f 100644
--- a/documentation/STIR-developers-overview.tex
+++ b/documentation/STIR-developers-overview.tex
@@ -1331,7 +1331,10 @@ \section{
\section{Contributing to STIR}
See
-\R2Lurl{https://github.com/UCL/STIR/blob/master/CONTRIBUTING.md}{STIR/CONTRIBUTING.md}.
+\R2Lurl{https://github.com/UCL/STIR/blob/master/CONTRIBUTING.md}{STIR/CONTRIBUTING.md}
+and
+\R2Lurl{https://github.com/UCL/STIR/blob/master/documentation/devel/README.md}{documentation/devel/README.md}
+for more information on editor settings etc.
\subsection{Continuous integration testing}
We use GitHub Actions and Appveyor for testing of every pull-request that was submitted on GitHub. PRs
diff --git a/documentation/release_6.0.htm b/documentation/release_6.0.htm
index 9c4e249ce..38cf9921e 100644
--- a/documentation/release_6.0.htm
+++ b/documentation/release_6.0.htm
@@ -9,7 +9,7 @@
Summary of changes in STIR release 6.0
This version is 99% backwards compatible with STIR 5.x for the user (see below).
Developers might need to make code changes as
- detailed below. Note though that the location of installed files have changed.
+ detailed below. Note though that the locations of installed files have changed.
Developers of other software that uses STIR via CMake will need to adapt (see below).
Overall summary
@@ -285,6 +285,18 @@ New deprecations for future versions
What's new for developers (aside from what should be obvious
from the above):
+ White-space and style enforcement
+
+ - We now use clang-format to enforce C++-style, including white-space settings, line-breaks
+ etc. This uses the .clang-format file in the root directory of STIR. Developers should
+ configure their editor encordingly, and ideally use pre-commit. It also has
+ consequences for existing branches as you might experience more conflicts than usual
+ during a merge. More detail is in
+ documentation/devel/README.md.
+ PR #1368.
+
+
+
Backward incompatibities
-
diff --git a/examples/C++/General_Reconstruction/General_Reconstruction.cxx b/examples/C++/General_Reconstruction/General_Reconstruction.cxx
index 277dc792f..1926455b8 100644
--- a/examples/C++/General_Reconstruction/General_Reconstruction.cxx
+++ b/examples/C++/General_Reconstruction/General_Reconstruction.cxx
@@ -5,53 +5,50 @@
#include
START_NAMESPACE_STIR
-General_Reconstruction::
-General_Reconstruction()
+General_Reconstruction::General_Reconstruction()
{
- this->set_defaults();
+ this->set_defaults();
}
void
General_Reconstruction::set_defaults()
-{
-
-}
+{}
void
General_Reconstruction::initialise_keymap()
{
- this->parser.add_start_key("General reconstruction");
- this->parser.add_stop_key("End General reconstruction");
+ this->parser.add_start_key("General reconstruction");
+ this->parser.add_stop_key("End General reconstruction");
- this->parser.add_parsing_key("reconstruction method", &this->reconstruction_method_sptr);
+ this->parser.add_parsing_key("reconstruction method", &this->reconstruction_method_sptr);
}
bool
General_Reconstruction::post_processing()
{
- return false;
+ return false;
}
Succeeded
General_Reconstruction::process_data()
{
- HighResWallClockTimer t;
- t.reset();
- t.start();
-
- //return reconstruction_object.reconstruct() == Succeeded::yes ?
- // EXIT_SUCCESS : EXIT_FAILURE;
- if (reconstruction_method_sptr->reconstruct() == Succeeded::yes)
- {
- t.stop();
- std::cout << "Total Wall clock time: " << t.value() << " seconds" << std::endl;
- return Succeeded::yes;
- }
- else
- {
- t.stop();
- return Succeeded::no;
- }
+ HighResWallClockTimer t;
+ t.reset();
+ t.start();
+
+ // return reconstruction_object.reconstruct() == Succeeded::yes ?
+ // EXIT_SUCCESS : EXIT_FAILURE;
+ if (reconstruction_method_sptr->reconstruct() == Succeeded::yes)
+ {
+ t.stop();
+ std::cout << "Total Wall clock time: " << t.value() << " seconds" << std::endl;
+ return Succeeded::yes;
+ }
+ else
+ {
+ t.stop();
+ return Succeeded::no;
+ }
}
END_NAMESPACE_STIR
diff --git a/examples/C++/General_Reconstruction/General_Reconstruction.h b/examples/C++/General_Reconstruction/General_Reconstruction.h
index a97321c60..b49fa29ab 100644
--- a/examples/C++/General_Reconstruction/General_Reconstruction.h
+++ b/examples/C++/General_Reconstruction/General_Reconstruction.h
@@ -16,7 +16,6 @@
#include "stir/CartesianCoordinate3D.h"
#include "Reconstruction.h"
-
START_NAMESPACE_STIR
class Succeeded;
@@ -24,23 +23,20 @@ class Succeeded;
class General_Reconstruction : public ParsingObject
{
public:
- //!
- //! \brief General_Reconstuction
- //! \details Default constructor
- General_Reconstruction();
+ //!
+ //! \brief General_Reconstuction
+ //! \details Default constructor
+ General_Reconstruction();
- virtual Succeeded process_data();
-protected:
+ virtual Succeeded process_data();
- void set_defaults();
- void initialise_keymap();
- bool post_processing();
+protected:
+ void set_defaults();
+ void initialise_keymap();
+ bool post_processing();
private:
-
- shared_ptr < Reconstruction < DiscretisedDensity < 3, float > > >
- reconstruction_method_sptr;
-
+ shared_ptr>> reconstruction_method_sptr;
};
END_NAMESPACE_STIR
diff --git a/examples/C++/using_STIR_LOCAL/demo1.cxx b/examples/C++/using_STIR_LOCAL/demo1.cxx
index 8162c8ded..42f279132 100644
--- a/examples/C++/using_STIR_LOCAL/demo1.cxx
+++ b/examples/C++/using_STIR_LOCAL/demo1.cxx
@@ -6,20 +6,20 @@
\brief A simple program that backprojects some projection data.
It illustrates
- - basic interaction with the user,
- - reading of images and projection data
- - construction of a specified type of back-projector,
- - how to use back-project all projection data
- - output of images
+ - basic interaction with the user,
+ - reading of images and projection data
+ - construction of a specified type of back-projector,
+ - how to use back-project all projection data
+ - output of images
See README.txt in the directory where this file is located.
- \author Kris Thielemans
+ \author Kris Thielemans
*/
/*
Copyright (C) 2004- 2011, Hammersmith Imanet Ltd
- This software is distributed under the terms
+ This software is distributed under the terms
of the GNU General Public Licence (GPL)
See STIR/LICENSE.txt for details
*/
@@ -35,25 +35,21 @@
#include "stir/utilities.h"
#include "stir/Succeeded.h"
-int main()
+int
+main()
{
using namespace stir;
-
+
/////////////// input sinogram
- const std::string input_filename =
- ask_filename_with_extension("Input file",".hs");
+ const std::string input_filename = ask_filename_with_extension("Input file", ".hs");
- shared_ptr
- proj_data_sptr(ProjData::read_from_file(input_filename));
- shared_ptr
- proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
+ shared_ptr proj_data_sptr(ProjData::read_from_file(input_filename));
+ shared_ptr proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
/////////////// template image (for sizes etc)
- const std::string template_filename =
- ask_filename_with_extension("Template image file",".hv");
+ const std::string template_filename = ask_filename_with_extension("Template image file", ".hv");
- shared_ptr >
- density_sptr(read_from_file >(template_filename));
+ shared_ptr> density_sptr(read_from_file>(template_filename));
density_sptr->fill(0);
diff --git a/examples/C++/using_STIR_LOCAL/demo2.cxx b/examples/C++/using_STIR_LOCAL/demo2.cxx
index 411c3fbc4..68fd828b4 100644
--- a/examples/C++/using_STIR_LOCAL/demo2.cxx
+++ b/examples/C++/using_STIR_LOCAL/demo2.cxx
@@ -4,23 +4,23 @@
\file
\ingroup examples
\brief A small modification of demo1.cxx to ask the user for the
- back projector she wants to use.
+ back projector she wants to use.
It illustrates
- - how to ask the user for objects for which different types
- exist (e.g. back-projector, forward-projectors, image processors
- etc), anything based on the RegisteredObject hierarchy.
- - that STIR is able to select basic processing units at run-time
- - how to use the (very) basic display facilities in STIR
+ - how to ask the user for objects for which different types
+ exist (e.g. back-projector, forward-projectors, image processors
+ etc), anything based on the RegisteredObject hierarchy.
+ - that STIR is able to select basic processing units at run-time
+ - how to use the (very) basic display facilities in STIR
See README.txt in the directory where this file is located.
- \author Kris Thielemans
+ \author Kris Thielemans
*/
/*
Copyright (C) 2004- 2012, Hammersmith Imanet Ltd
- This software is distributed under the terms
+ This software is distributed under the terms
of the GNU General Public Licence (GPL)
See STIR/LICENSE.txt for details
*/
@@ -34,31 +34,26 @@
#include "stir/Succeeded.h"
#include "stir/display.h"
-int main()
+int
+main()
{
using namespace stir;
-
+
/////////////// input sinogram
- const std::string input_filename =
- ask_filename_with_extension("Input file",".hs");
+ const std::string input_filename = ask_filename_with_extension("Input file", ".hs");
- shared_ptr
- proj_data_sptr(ProjData::read_from_file(input_filename));
- shared_ptr
- proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
+ shared_ptr proj_data_sptr(ProjData::read_from_file(input_filename));
+ shared_ptr proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
/////////////// template image (for sizes etc)
- const std::string template_filename =
- ask_filename_with_extension("Template image file",".hv");
+ const std::string template_filename = ask_filename_with_extension("Template image file", ".hv");
- shared_ptr >
- density_sptr(read_from_file >(template_filename));
+ shared_ptr> density_sptr(read_from_file>(template_filename));
density_sptr->fill(0);
/////////////// back project
- shared_ptr back_projector_sptr
- (BackProjectorByBin::ask_type_and_parameters());
+ shared_ptr back_projector_sptr(BackProjectorByBin::ask_type_and_parameters());
back_projector_sptr->set_up(proj_data_info_sptr, density_sptr);
diff --git a/examples/C++/using_STIR_LOCAL/demo3.cxx b/examples/C++/using_STIR_LOCAL/demo3.cxx
index ec2ecba3d..4df712be0 100644
--- a/examples/C++/using_STIR_LOCAL/demo3.cxx
+++ b/examples/C++/using_STIR_LOCAL/demo3.cxx
@@ -6,10 +6,10 @@
\brief A modification of demo2.cxx that parses all parameters from a parameter file.
It illustrates
- - basic class derivation principles
- - how to use ParsingObject to have automatic capabilities of parsing
- parameters files (and interactive questions to the user)
- - how most STIR programs parse the parameter files.
+ - basic class derivation principles
+ - how to use ParsingObject to have automatic capabilities of parsing
+ parameters files (and interactive questions to the user)
+ - how most STIR programs parse the parameter files.
Note that the same functionality could be provided without deriving
a new class from ParsingObject. One could have a KeyParser object
@@ -17,12 +17,12 @@
See README.txt in the directory where this file is located.
- \author Kris Thielemans
+ \author Kris Thielemans
*/
/*
Copyright (C) 2004- 2012, Hammersmith Imanet Ltd
- This software is distributed under the terms
+ This software is distributed under the terms
of the GNU General Public Licence (GPL)
See STIR/LICENSE.txt for details
*/
@@ -39,9 +39,10 @@
#include "stir/display.h"
#include
-namespace stir {
+namespace stir
+{
-class MyStuff: public ParsingObject
+class MyStuff : public ParsingObject
{
public:
void set_defaults();
@@ -52,7 +53,7 @@ class MyStuff: public ParsingObject
std::string input_filename;
std::string template_filename;
shared_ptr back_projector_sptr;
- shared_ptr > > output_file_format_sptr;
+ shared_ptr>> output_file_format_sptr;
};
void
@@ -60,10 +61,10 @@ MyStuff::set_defaults()
{
auto projection_matrix_sptr = std::make_shared();
back_projector_sptr = std::make_shared(projection_matrix_sptr);
- output_file_format_sptr = OutputFileFormat >::default_sptr();
+ output_file_format_sptr = OutputFileFormat>::default_sptr();
}
-void
+void
MyStuff::initialise_keymap()
{
parser.add_start_key("MyStuff parameters");
@@ -78,13 +79,10 @@ void
MyStuff::run(const bool display_off)
{
- shared_ptr
- proj_data_sptr(ProjData::read_from_file(input_filename));
- shared_ptr
- proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
+ shared_ptr proj_data_sptr(ProjData::read_from_file(input_filename));
+ shared_ptr proj_data_info_sptr(proj_data_sptr->get_proj_data_info_sptr()->clone());
- shared_ptr >
- density_sptr(read_from_file >(template_filename));
+ shared_ptr> density_sptr(read_from_file>(template_filename));
density_sptr->fill(0);
@@ -100,28 +98,29 @@ MyStuff::run(const bool display_off)
display(*density_sptr, density_sptr->find_max(), "Output");
}
-}// end of namespace stir
+} // end of namespace stir
-int main(int argc, char **argv)
+int
+main(int argc, char** argv)
{
using namespace stir;
MyStuff my_stuff;
my_stuff.set_defaults();
bool display_off = false;
- if (argc<2)
- {
- std::cerr << "Normal usage: " << argv[0] << " parameter-file [--display_off]\n";
- std::cerr << "I will now ask you the questions interactively\n";
- my_stuff.ask_parameters();
- }
+ if (argc < 2)
+ {
+ std::cerr << "Normal usage: " << argv[0] << " parameter-file [--display_off]\n";
+ std::cerr << "I will now ask you the questions interactively\n";
+ my_stuff.ask_parameters();
+ }
else
- {
- my_stuff.parse(argv[1]);
- if (argc>=3)
- // Set the display_off to true if the second argument is "--display_off"
- display_off = (strcmp(argv[2], "--display_off")==0);
- }
+ {
+ my_stuff.parse(argv[1]);
+ if (argc >= 3)
+ // Set the display_off to true if the second argument is "--display_off"
+ display_off = (strcmp(argv[2], "--display_off") == 0);
+ }
my_stuff.run(display_off);
return EXIT_SUCCESS;
}
diff --git a/examples/C++/using_STIR_LOCAL/demo4_obj_fun.cxx b/examples/C++/using_STIR_LOCAL/demo4_obj_fun.cxx
index 041fc9a5e..791b6de94 100644
--- a/examples/C++/using_STIR_LOCAL/demo4_obj_fun.cxx
+++ b/examples/C++/using_STIR_LOCAL/demo4_obj_fun.cxx
@@ -21,12 +21,12 @@
See README.txt in the directory where this file is located.
- \author Kris Thielemans and Robert Twyman
+ \author Kris Thielemans and Robert Twyman
*/
/*
Copyright (C) 2020 University College London
- This software is distributed under the terms
+ This software is distributed under the terms
of the GNU General Public Licence (GPL)
See STIR/LICENSE.txt for details
*/
@@ -37,25 +37,26 @@
#include "stir/is_null_ptr.h"
#include "stir/error.h"
-namespace stir {
+namespace stir
+{
-class MyStuff: public ParsingObject
+class MyStuff : public ParsingObject
{
public:
MyStuff();
void set_defaults();
void run();
- typedef DiscretisedDensity<3,float> target_type;
+ typedef DiscretisedDensity<3, float> target_type;
protected:
- shared_ptr > objective_function_sptr;
+ shared_ptr> objective_function_sptr;
private:
std::string input_filename;
std::string image_filename;
int num_iterations;
float step_size;
- shared_ptr > > output_file_format_sptr;
+ shared_ptr>> output_file_format_sptr;
void initialise_keymap();
bool post_processing();
};
@@ -69,12 +70,12 @@ void
MyStuff::set_defaults()
{
objective_function_sptr.reset(new PoissonLogLikelihoodWithLinearModelForMeanAndProjData);
- output_file_format_sptr = OutputFileFormat >::default_sptr();
+ output_file_format_sptr = OutputFileFormat>::default_sptr();
num_iterations = 10;
step_size = 0.001;
}
-void
+void
MyStuff::initialise_keymap()
{
parser.add_start_key("MyStuff parameters");
@@ -86,14 +87,14 @@ MyStuff::initialise_keymap()
parser.add_stop_key("End");
}
-bool MyStuff::
-post_processing()
+bool
+MyStuff::post_processing()
{
if (is_null_ptr(this->objective_function_sptr))
- {
+ {
error("objective_function_sptr is null");
return true;
- }
+ }
return false;
}
@@ -101,12 +102,10 @@ void
MyStuff::run()
{
/////// load initial density from file
- shared_ptr >
- density_sptr(read_from_file >(image_filename));
+ shared_ptr> density_sptr(read_from_file>(image_filename));
//////// gradient it copied Density filled with 0's
- shared_ptr >
- gradient_sptr(density_sptr->get_empty_copy());
+ shared_ptr> gradient_sptr(density_sptr->get_empty_copy());
/////// setup objective function object
objective_function_sptr->set_up(density_sptr);
@@ -117,41 +116,43 @@ MyStuff::run()
/////// Iteratively compute and add objective function gradients to density_sptr
/////// Update formula: x_{k+1} = x_k + step_size * gradient
/////// Note the lack of positivity constraint on the images.
- for (int k = 0; k < num_iterations; k++) {
- std::cout << "Iteration number: " << k << "\n";
- objective_function_sptr->compute_sub_gradient(*gradient_sptr, *density_sptr, 0);
- *density_sptr += (*gradient_sptr) * step_size;
-
- /////// save density and estimates
- std::string density_filename = "demo4_density_" + std::to_string(k);
- std::string gradient_filename = "demo4_gradient_" + std::to_string(k);
- output_file_format_sptr->write_to_file(density_filename, *density_sptr);
- output_file_format_sptr->write_to_file(gradient_filename, *gradient_sptr);
- }
+ for (int k = 0; k < num_iterations; k++)
+ {
+ std::cout << "Iteration number: " << k << "\n";
+ objective_function_sptr->compute_sub_gradient(*gradient_sptr, *density_sptr, 0);
+ *density_sptr += (*gradient_sptr) * step_size;
+
+ /////// save density and estimates
+ std::string density_filename = "demo4_density_" + std::to_string(k);
+ std::string gradient_filename = "demo4_gradient_" + std::to_string(k);
+ output_file_format_sptr->write_to_file(density_filename, *density_sptr);
+ output_file_format_sptr->write_to_file(gradient_filename, *gradient_sptr);
+ }
/////// Compute objective function value of image + gradient
const double my_objective_function_value2 = objective_function_sptr->compute_objective_function(*density_sptr);
/////// Return the objective function values and improvement
std::cout << "The initial Objective Function Value = " << my_objective_function_value1 << "\n";
- std::cout << "The Objective Function Value of after " << num_iterations << " iteration(s) ="
- << my_objective_function_value2 << "\n";
+ std::cout << "The Objective Function Value of after " << num_iterations << " iteration(s) =" << my_objective_function_value2
+ << "\n";
std::cout << "A change of " << my_objective_function_value2 - my_objective_function_value1 << "\n";
}
-}// end of namespace stir
+} // end of namespace stir
-int main(int argc, char **argv)
+int
+main(int argc, char** argv)
{
using namespace stir;
- if (argc!=2)
+ if (argc != 2)
{
std::cerr << "Normal usage: " << argv[0] << " parameter-file\n";
std::cerr << "I will now ask you the questions interactively\n";
}
MyStuff my_stuff;
my_stuff.set_defaults();
- if (argc!=2)
+ if (argc != 2)
my_stuff.ask_parameters();
else
my_stuff.parse(argv[1]);
diff --git a/examples/C++/using_STIR_LOCAL/demo5_line_search.cxx b/examples/C++/using_STIR_LOCAL/demo5_line_search.cxx
index b40a3ea8d..66a9315f7 100644
--- a/examples/C++/using_STIR_LOCAL/demo5_line_search.cxx
+++ b/examples/C++/using_STIR_LOCAL/demo5_line_search.cxx
@@ -48,9 +48,8 @@ compute_linear_alphas(const float alpha_min, const float alpha_max, const float
float d_alpha = (alpha_max - alpha_min) / num_evaluations;
std::cout << "\nComputing linear alphas:"
- "\n alpha_min = " << alpha_min <<
- "\n alpha_max = " << alpha_max <<
- "\n delta_alpha = " << d_alpha << "\n";
+ "\n alpha_min = "
+ << alpha_min << "\n alpha_max = " << alpha_max << "\n delta_alpha = " << d_alpha << "\n";
/// Explicitly add alpha = 0.0 and/or alpha_min
alphas.push_back(0.0);
@@ -64,7 +63,6 @@ compute_linear_alphas(const float alpha_min, const float alpha_max, const float
return alphas;
}
-
std::vector
compute_exponential_alphas(const float alpha_min, const float alpha_max, const float num_evaluations)
{
@@ -73,9 +71,8 @@ compute_exponential_alphas(const float alpha_min, const float alpha_max, const f
float d_alpha = (alpha_max - alpha_min) / num_evaluations;
std::cout << "\nComputing exponential alphas:"
- "\n exponential min = " << alpha_min <<
- "\n exponential max = " << alpha_max <<
- "\n exponential delta = " << d_alpha << "\n";
+ "\n exponential min = "
+ << alpha_min << "\n exponential max = " << alpha_max << "\n exponential delta = " << d_alpha << "\n";
/// Explicitly add alpha = 0.0 and/or alpha_min
alphas.push_back(0.0);
@@ -87,65 +84,65 @@ compute_exponential_alphas(const float alpha_min, const float alpha_max, const f
return alphas;
}
-
void
save_doubles_vector_to_file(std::string filename, std::vector vector)
{
/// This function is used to save the line search results (alpha and Phi values) to separate files.
- std::ofstream myfile (filename);
+ std::ofstream myfile(filename);
int precision = 40;
- if (myfile.is_open()){
- for (double v : vector){
- myfile << std::fixed << std::setprecision(precision) << v << std::endl;
+ if (myfile.is_open())
+ {
+ for (double v : vector)
+ {
+ myfile << std::fixed << std::setprecision(precision) << v << std::endl;
+ }
+ myfile.close();
}
- myfile.close();
- }
}
-
using namespace stir;
-class LineSearcher: public ParsingObject
+class LineSearcher : public ParsingObject
{
public:
- LineSearcher();
- /// Methods
- void set_defaults();
- void setup();
- double compute_line_search_value(const double alpha);
- void apply_update_step(const double alpha);
- void perform_line_search();
- void save_data();
- void save_max_line_search_image();
-
- typedef DiscretisedDensity<3,float> target_type;
-
- /// Class variables
- int num_evaluations;
- float alpha_min;
- float alpha_max;
- bool use_exponential_alphas;
- double image_lower_bound;
-
- /// Measurements
- std::vector alphas;
- std::vector Phis;
-
- /// Image volumes
- shared_ptr > image_sptr;
- shared_ptr > gradient_sptr;
- shared_ptr > eval_image_sptr;
+ LineSearcher();
+ /// Methods
+ void set_defaults();
+ void setup();
+ double compute_line_search_value(const double alpha);
+ void apply_update_step(const double alpha);
+ void perform_line_search();
+ void save_data();
+ void save_max_line_search_image();
+
+ typedef DiscretisedDensity<3, float> target_type;
+
+ /// Class variables
+ int num_evaluations;
+ float alpha_min;
+ float alpha_max;
+ bool use_exponential_alphas;
+ double image_lower_bound;
+
+ /// Measurements
+ std::vector alphas;
+ std::vector Phis;
+
+ /// Image volumes
+ shared_ptr> image_sptr;
+ shared_ptr> gradient_sptr;
+ shared_ptr> eval_image_sptr;
protected:
- shared_ptr > objective_function_sptr;
+ shared_ptr> objective_function_sptr;
private:
- void initialise_keymap();
- bool post_processing();
+ void initialise_keymap();
+ bool post_processing();
- std::string image_filename;
- bool is_setup;
- shared_ptr > > output_file_format_sptr;
+ std::string image_filename;
+ bool is_setup;
+ shared_ptr>> output_file_format_sptr;
};
LineSearcher::LineSearcher()
@@ -153,7 +150,6 @@ LineSearcher::LineSearcher()
set_defaults();
}
-
void
LineSearcher::set_defaults()
{
@@ -164,10 +160,9 @@ LineSearcher::set_defaults()
use_exponential_alphas = false;
is_setup = false;
image_lower_bound = 0.0;
- output_file_format_sptr = OutputFileFormat >::default_sptr();
+ output_file_format_sptr = OutputFileFormat>::default_sptr();
}
-
void
LineSearcher::initialise_keymap()
{
@@ -182,21 +177,19 @@ LineSearcher::initialise_keymap()
parser.add_stop_key("End");
}
-
-bool LineSearcher::
-post_processing()
+bool
+LineSearcher::post_processing()
{
if (is_null_ptr(this->objective_function_sptr))
- {
- error("objective_function_sptr is null");
- return true;
- }
+ {
+ error("objective_function_sptr is null");
+ return true;
+ }
return false;
}
-
-void LineSearcher::
-setup()
+void
+LineSearcher::setup()
{
/// Setup LineSearcher
this->is_setup = false;
@@ -206,7 +199,7 @@ setup()
error("LineSearcher setup. No image filename has been given.");
std::cout << "Loading image: \n " << image_filename << "\n";
- this->image_sptr = read_from_file >(image_filename);
+ this->image_sptr = read_from_file>(image_filename);
/// Gradient it copied density filled with 0's
this->gradient_sptr.reset(this->image_sptr->get_empty_copy());
@@ -222,9 +215,9 @@ setup()
this->is_setup = true;
}
-
void
-LineSearcher::perform_line_search() {
+LineSearcher::perform_line_search()
+{
/// Performs the line search
/// Gets the step sizes
/// Computes the objective function value of the image at every step size
@@ -234,12 +227,12 @@ LineSearcher::perform_line_search() {
error("LineSearcher is not setup, please run setup()");
double Phi; // Used to store objective function values of each iterate
- std::cout << "Computing objective function values of alphas from " << this->alpha_min << " to "
- << this->alpha_max << " in increments of " << this->num_evaluations << "\n";
+ std::cout << "Computing objective function values of alphas from " << this->alpha_min << " to " << this->alpha_max
+ << " in increments of " << this->num_evaluations << "\n";
/// get alpha values as a vector
{
- if ( this->use_exponential_alphas )
+ if (this->use_exponential_alphas)
alphas = compute_exponential_alphas(this->alpha_min, this->alpha_max, this->num_evaluations);
else
alphas = compute_linear_alphas(this->alpha_min, this->alpha_max, this->num_evaluations);
@@ -247,29 +240,26 @@ LineSearcher::perform_line_search() {
/// Iterate over each of the alphas and compute Phi
for (auto a = alphas.begin(); a != alphas.end(); ++a)
- {
- Phi = this->compute_line_search_value(*a);
- Phis.push_back(Phi);
- std::cout << "alpha = " << *a << ". Phi = " << Phi << "\n";
- }
+ {
+ Phi = this->compute_line_search_value(*a);
+ Phis.push_back(Phi);
+ std::cout << "alpha = " << *a << ". Phi = " << Phi << "\n";
+ }
/// Output alpha and Phi values to console
std::cout << "\n\n====================================\n"
"Alpha and Phi values: \n";
- for (int i = 0 ; i < alphas.size() ; ++i)
+ for (int i = 0; i < alphas.size(); ++i)
std::cout << std::setprecision(20) << " alpha = " << alphas[i] << ". Phis = " << Phis[i] << "\n";
-
}
-
double
LineSearcher::compute_line_search_value(const double alpha)
{
/// For a given alpha, computes the objective function value at the update step
apply_update_step(alpha);
- std::cout << "\nimage_min = " << image_sptr->find_min()
- << "\ngrad_min = " << gradient_sptr->find_min()
+ std::cout << "\nimage_min = " << image_sptr->find_min() << "\ngrad_min = " << gradient_sptr->find_min()
<< "\neval_min = " << eval_image_sptr->find_min() << "\n";
return objective_function_sptr->compute_objective_function(*eval_image_sptr);
}
@@ -283,7 +273,6 @@ LineSearcher::apply_update_step(const double alpha)
this->eval_image_sptr->apply_lower_threshold(this->image_lower_bound);
}
-
void
LineSearcher::save_data()
{
@@ -299,7 +288,6 @@ LineSearcher::save_data()
output_file_format_sptr->write_to_file("LineSearchGradient.hv", *this->gradient_sptr);
}
-
void
LineSearcher::save_max_line_search_image()
{
@@ -308,37 +296,41 @@ LineSearcher::save_max_line_search_image()
double max_alpha = alphas[0];
/// Find the alpha Phi combination that is max Phi and save that image.
- for (int i = 0 ; i < alphas.size() ; ++i){
- if (Phis[i] > max_Phi){
- max_Phi = Phis[i];
- max_alpha = alphas[i];
+ for (int i = 0; i < alphas.size(); ++i)
+ {
+ if (Phis[i] > max_Phi)
+ {
+ max_Phi = Phis[i];
+ max_alpha = alphas[i];
+ }
}
- }
/// Check if alpha == 0 is optimal, otherwise save the image at the maximum evaluation.
if (max_alpha != alphas[0])
- {
- apply_update_step(max_alpha);
- std::cout << "Saving max line search value image, computed at alpha = " << max_alpha << "\n"
- << "and Phi = " << max_Phi << '\n';
- output_file_format_sptr->write_to_file("MaxObjectiveFunctionImage.hv", *this->eval_image_sptr);
- } else {
- std::cout << "Max line search value image is at alpha = 0.0\n";
- }
+ {
+ apply_update_step(max_alpha);
+ std::cout << "Saving max line search value image, computed at alpha = " << max_alpha << "\n"
+ << "and Phi = " << max_Phi << '\n';
+ output_file_format_sptr->write_to_file("MaxObjectiveFunctionImage.hv", *this->eval_image_sptr);
+ }
+ else
+ {
+ std::cout << "Max line search value image is at alpha = 0.0\n";
+ }
}
-
-int main(int argc, char **argv)
+int
+main(int argc, char** argv)
{
using namespace stir;
- if (argc!=2)
- {
- std::cerr << "Normal usage: " << argv[0] << " parameter-file\n";
- std::cerr << "I will now ask you the questions interactively\n";
- }
+ if (argc != 2)
+ {
+ std::cerr << "Normal usage: " << argv[0] << " parameter-file\n";
+ std::cerr << "I will now ask you the questions interactively\n";
+ }
LineSearcher my_stuff;
- if (argc!=2)
+ if (argc != 2)
my_stuff.ask_parameters();
else
my_stuff.parse(argv[1]);
diff --git a/examples/C++/using_installed_STIR/demo_create_image.cxx b/examples/C++/using_installed_STIR/demo_create_image.cxx
index f3d391c9b..ffc3c53ea 100644
--- a/examples/C++/using_installed_STIR/demo_create_image.cxx
+++ b/examples/C++/using_installed_STIR/demo_create_image.cxx
@@ -25,7 +25,8 @@
#include "stir/VoxelsOnCartesianGrid.h"
#include
-int main()
+int
+main()
{
std::string output_filename("test.hv");
@@ -36,7 +37,7 @@ int main()
float output_voxel_size_y = 4.F;
float output_voxel_size_z = 3.F;
double image_duration = 100.; // secs
- double rel_start_time = 1; //secs
+ double rel_start_time = 1; // secs
auto exam_info_sptr = std::make_shared(stir::ImagingModality::PT);
{
@@ -45,24 +46,22 @@ int main()
frame_defs.set_time_frame(1, rel_start_time, image_duration);
exam_info_sptr->set_time_frame_definitions(frame_defs);
}
- stir::VoxelsOnCartesianGrid
- image(exam_info_sptr,
- stir::IndexRange3D(0,output_image_size_z-1,
- -(output_image_size_y/2),
- -(output_image_size_y/2)+output_image_size_y-1,
- -(output_image_size_x/2),
- -(output_image_size_x/2)+output_image_size_x-1),
- stir::CartesianCoordinate3D(0,0,0),
- stir::CartesianCoordinate3D(output_voxel_size_z,
- output_voxel_size_y,
- output_voxel_size_x));
+ stir::VoxelsOnCartesianGrid image(
+ exam_info_sptr,
+ stir::IndexRange3D(0,
+ output_image_size_z - 1,
+ -(output_image_size_y / 2),
+ -(output_image_size_y / 2) + output_image_size_y - 1,
+ -(output_image_size_x / 2),
+ -(output_image_size_x / 2) + output_image_size_x - 1),
+ stir::CartesianCoordinate3D(0, 0, 0),
+ stir::CartesianCoordinate3D(output_voxel_size_z, output_voxel_size_y, output_voxel_size_x));
// add shape to image
{
- const auto centre =
- image.get_physical_coordinates_for_indices(stir::make_coordinate(output_image_size_z/2,0,0));
+ const auto centre = image.get_physical_coordinates_for_indices(stir::make_coordinate(output_image_size_z / 2, 0, 0));
stir::EllipsoidalCylinder shape(40.F, 30.F, 20.F, centre);
- shape.construct_volume(image, stir::make_coordinate(2,2,2));
+ shape.construct_volume(image, stir::make_coordinate(2, 2, 2));
}
// write output for checking
@@ -75,5 +74,5 @@ int main()
{
return EXIT_FAILURE;
}
- return EXIT_SUCCESS;
+ return EXIT_SUCCESS;
}
diff --git a/src/IO/ECAT6OutputFileFormat.cxx b/src/IO/ECAT6OutputFileFormat.cxx
index 68bbd0762..89bc80537 100644
--- a/src/IO/ECAT6OutputFileFormat.cxx
+++ b/src/IO/ECAT6OutputFileFormat.cxx
@@ -29,22 +29,17 @@ START_NAMESPACE_STIR
START_NAMESPACE_ECAT
START_NAMESPACE_ECAT6
+const char* const ECAT6OutputFileFormat::registered_name = "ECAT6";
-const char * const
-ECAT6OutputFileFormat::registered_name = "ECAT6";
-
-ECAT6OutputFileFormat::
-ECAT6OutputFileFormat(const NumericType& type,
- const ByteOrder& byte_order)
+ECAT6OutputFileFormat::ECAT6OutputFileFormat(const NumericType& type, const ByteOrder& byte_order)
{
base_type::set_defaults();
set_type_of_numbers(type);
set_byte_order(byte_order);
}
-void
-ECAT6OutputFileFormat::
-initialise_keymap()
+void
+ECAT6OutputFileFormat::initialise_keymap()
{
parser.add_start_key("ECAT6 Output File Format Parameters");
parser.add_stop_key("End ECAT6 Output File Format Parameters");
@@ -52,103 +47,84 @@ initialise_keymap()
base_type::initialise_keymap();
}
-void
-ECAT6OutputFileFormat::
-set_defaults()
+void
+ECAT6OutputFileFormat::set_defaults()
{
default_scanner_name = "ECAT 953";
base_type::set_defaults();
file_byte_order = ByteOrder::little_endian;
type_of_numbers = NumericType::SHORT;
-
}
bool
-ECAT6OutputFileFormat::
-post_processing()
+ECAT6OutputFileFormat::post_processing()
{
if (base_type::post_processing())
return true;
- shared_ptr scanner_ptr (
- Scanner::get_scanner_from_name(default_scanner_name));
+ shared_ptr scanner_ptr(Scanner::get_scanner_from_name(default_scanner_name));
- if (find_ECAT_system_type(*scanner_ptr)==0)
+ if (find_ECAT_system_type(*scanner_ptr) == 0)
{
- warning("ECAT6OutputFileFormat: default_scanner_name %s is not supported\n",
- default_scanner_name.c_str());
+ warning("ECAT6OutputFileFormat: default_scanner_name %s is not supported\n", default_scanner_name.c_str());
return true;
}
return false;
}
-NumericType
-ECAT6OutputFileFormat::
-set_type_of_numbers(const NumericType& new_type, const bool warn)
+NumericType
+ECAT6OutputFileFormat::set_type_of_numbers(const NumericType& new_type, const bool warn)
{
- const NumericType supported_type_of_numbers =
- NumericType("signed integer", 2);
+ const NumericType supported_type_of_numbers = NumericType("signed integer", 2);
if (new_type != supported_type_of_numbers)
- {
- if (warn)
- warning("ECAT6OutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
- type_of_numbers = supported_type_of_numbers;
- }
+ {
+ if (warn)
+ warning("ECAT6OutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
+ type_of_numbers = supported_type_of_numbers;
+ }
else
type_of_numbers = new_type;
return type_of_numbers;
-
}
-ByteOrder
-ECAT6OutputFileFormat::
-set_byte_order(const ByteOrder& new_byte_order, const bool warn)
+ByteOrder
+ECAT6OutputFileFormat::set_byte_order(const ByteOrder& new_byte_order, const bool warn)
{
if (new_byte_order != ByteOrder::little_endian)
- {
- if (warn)
- warning("ECAT6OutputFileFormat: byte_order is currently fixed to little-endian\n");
- file_byte_order = ByteOrder::little_endian;
- }
+ {
+ if (warn)
+ warning("ECAT6OutputFileFormat: byte_order is currently fixed to little-endian\n");
+ file_byte_order = ByteOrder::little_endian;
+ }
else
file_byte_order = new_byte_order;
return file_byte_order;
}
-Succeeded
-ECAT6OutputFileFormat::
-actual_write_to_file(std::string& filename,
- const DiscretisedDensity<3,float>& density) const
+Succeeded
+ECAT6OutputFileFormat::actual_write_to_file(std::string& filename, const DiscretisedDensity<3, float>& density) const
{
- shared_ptr scanner_ptr(
- Scanner::get_scanner_from_name(default_scanner_name));
-
+ shared_ptr scanner_ptr(Scanner::get_scanner_from_name(default_scanner_name));
+
add_extension(filename, ".img");
ECAT6_Main_header mhead;
make_ECAT6_Main_header(mhead, *scanner_ptr, "", density);
mhead.num_frames = 1;
-
- FILE *fptr= cti_create (filename.c_str(), &mhead);
+
+ FILE* fptr = cti_create(filename.c_str(), &mhead);
if (fptr == NULL)
- {
- warning("ECAT6OutputFileFormat::write_to_file: error opening output file %s\n",
- filename.c_str());
- return Succeeded::no;
- }
-
- Succeeded success =
- DiscretisedDensity_to_ECAT6(fptr,
- density,
- mhead,
- 1 /*frame_num*/);
+ {
+ warning("ECAT6OutputFileFormat::write_to_file: error opening output file %s\n", filename.c_str());
+ return Succeeded::no;
+ }
+
+ Succeeded success = DiscretisedDensity_to_ECAT6(fptr, density, mhead, 1 /*frame_num*/);
fclose(fptr);
-
+
return success;
}
END_NAMESPACE_ECAT6
END_NAMESPACE_ECAT
END_NAMESPACE_STIR
-
-
diff --git a/src/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.cxx b/src/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.cxx
index 6009369b6..e019db573 100644
--- a/src/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.cxx
+++ b/src/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.cxx
@@ -27,9 +27,8 @@
#include
#include
-
#ifndef HAVE_LLN_MATRIX
-#error HAVE_LLN_MATRIX not defined: you need the lln ecat library.
+# error HAVE_LLN_MATRIX not defined: you need the lln ecat library.
#endif
#include "stir/IO/stir_ecat7.h"
@@ -42,10 +41,8 @@ START_NAMESPACE_ECAT7
\preliminary
*/
-bool
-ECAT7DynamicDiscretisedDensityInputFileFormat::
-actual_can_read(const FileSignature& signature,
- std::istream& input) const
+bool
+ECAT7DynamicDiscretisedDensityInputFileFormat::actual_can_read(const FileSignature& signature, std::istream& input) const
{
if (strncmp(signature.get_signature(), "MATRIX", 6) != 0)
return false;
@@ -56,38 +53,32 @@ actual_can_read(const FileSignature& signature,
}
unique_ptr
-ECAT7DynamicDiscretisedDensityInputFileFormat::
-read_from_file(std::istream& input) const
+ECAT7DynamicDiscretisedDensityInputFileFormat::read_from_file(std::istream& input) const
{
- //TODO
- error("read_from_file for ECAT7 with istream not implemented %s:%s. Sorry",
- __FILE__, __LINE__);
- return
- unique_ptr();
+ // TODO
+ error("read_from_file for ECAT7 with istream not implemented %s:%s. Sorry", __FILE__, __LINE__);
+ return unique_ptr();
}
unique_ptr
-ECAT7DynamicDiscretisedDensityInputFileFormat::
-read_from_file(const std::string& filename) const
+ECAT7DynamicDiscretisedDensityInputFileFormat::read_from_file(const std::string& filename) const
{
if (is_ECAT7_image_file(filename))
{
Main_header mhead;
if (read_ECAT7_main_header(mhead, filename) == Succeeded::no)
- {
- error("ECAT7DynamicDiscretisedDensityInputFileFormat::read_from_file cannot read %s as ECAT7 (failed to read main header)", filename.c_str());
- return unique_ptr();
- }
+ {
+ error("ECAT7DynamicDiscretisedDensityInputFileFormat::read_from_file cannot read %s as ECAT7 (failed to read main "
+ "header)",
+ filename.c_str());
+ return unique_ptr();
+ }
TimeFrameDefinitions time_frame_definitions(filename);
shared_ptr scanner_sptr(find_scanner_from_ECAT_system_type(mhead.system_type));
- unique_ptr
- dynamic_image_ptr
- (new DynamicDiscretisedDensity(time_frame_definitions,
- static_cast(mhead.scan_start_time),
- scanner_sptr)
- );
+ unique_ptr dynamic_image_ptr(
+ new DynamicDiscretisedDensity(time_frame_definitions, static_cast(mhead.scan_start_time), scanner_sptr));
dynamic_image_ptr->set_calibration_factor(mhead.calibration_factor);
@@ -98,17 +89,18 @@ read_from_file(const std::string& filename) const
// shead.processing_code & DecayPrc
dynamic_image_ptr->set_if_decay_corrected(false);
- for (unsigned int frame_num=1; frame_num <= dynamic_image_ptr->get_num_time_frames(); ++ frame_num)
- {
- shared_ptr dens_sptr
- (ECAT7_to_VoxelsOnCartesianGrid(filename,
- frame_num,
- /* gate_num, data_num, bed_num */ 1,0,0)
- );
- if (is_null_ptr(dens_sptr))
- error("read_from_file for DynamicDiscretisedDensity: No frame %d available", frame_num);
- dynamic_image_ptr->set_density(*dens_sptr, frame_num );
- }
+ for (unsigned int frame_num = 1; frame_num <= dynamic_image_ptr->get_num_time_frames(); ++frame_num)
+ {
+ shared_ptr dens_sptr(
+ ECAT7_to_VoxelsOnCartesianGrid(filename,
+ frame_num,
+ /* gate_num, data_num, bed_num */ 1,
+ 0,
+ 0));
+ if (is_null_ptr(dens_sptr))
+ error("read_from_file for DynamicDiscretisedDensity: No frame %d available", frame_num);
+ dynamic_image_ptr->set_density(*dens_sptr, frame_num);
+ }
return dynamic_image_ptr;
}
else
@@ -122,4 +114,3 @@ read_from_file(const std::string& filename) const
END_NAMESPACE_ECAT
END_NAMESPACE_ECAT7
END_NAMESPACE_STIR
-
diff --git a/src/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.cxx b/src/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.cxx
index 8e3c592e4..3ed0b1343 100644
--- a/src/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.cxx
+++ b/src/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.cxx
@@ -28,21 +28,18 @@ START_NAMESPACE_STIR
START_NAMESPACE_ECAT
START_NAMESPACE_ECAT7
-const char * const
-ECAT7DynamicDiscretisedDensityOutputFileFormat::registered_name = "ECAT7";
+const char* const ECAT7DynamicDiscretisedDensityOutputFileFormat::registered_name = "ECAT7";
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-ECAT7DynamicDiscretisedDensityOutputFileFormat(const NumericType& type,
- const ByteOrder& byte_order)
+ECAT7DynamicDiscretisedDensityOutputFileFormat::ECAT7DynamicDiscretisedDensityOutputFileFormat(const NumericType& type,
+ const ByteOrder& byte_order)
{
this->set_defaults();
this->set_type_of_numbers(type);
this->set_byte_order(byte_order);
}
-void
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-initialise_keymap()
+void
+ECAT7DynamicDiscretisedDensityOutputFileFormat::initialise_keymap()
{
this->parser.add_start_key("ECAT7 Output File Format Parameters");
this->parser.add_stop_key("End ECAT7 Output File Format Parameters");
@@ -50,9 +47,8 @@ initialise_keymap()
base_type::initialise_keymap();
}
-void
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-set_defaults()
+void
+ECAT7DynamicDiscretisedDensityOutputFileFormat::set_defaults()
{
this->default_scanner_name = "ECAT 962";
base_type::set_defaults();
@@ -63,69 +59,62 @@ set_defaults()
}
bool
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-post_processing()
+ECAT7DynamicDiscretisedDensityOutputFileFormat::post_processing()
{
if (base_type::post_processing())
return true;
- shared_ptr scanner_ptr(
- Scanner::get_scanner_from_name(this->default_scanner_name));
+ shared_ptr scanner_ptr(Scanner::get_scanner_from_name(this->default_scanner_name));
- if (find_ECAT_system_type(*scanner_ptr)==0)
+ if (find_ECAT_system_type(*scanner_ptr) == 0)
{
warning("ECAT7DynamicDiscretisedDensityOutputFileFormat: default_scanner_name %s is not supported\n",
- this->default_scanner_name.c_str());
+ this->default_scanner_name.c_str());
return true;
}
return false;
}
-NumericType
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-set_type_of_numbers(const NumericType& new_type, const bool warn)
+NumericType
+ECAT7DynamicDiscretisedDensityOutputFileFormat::set_type_of_numbers(const NumericType& new_type, const bool warn)
{
- const NumericType supported_type_of_numbers =
- NumericType("signed integer", 2);
+ const NumericType supported_type_of_numbers = NumericType("signed integer", 2);
if (new_type != supported_type_of_numbers)
- {
- if (warn)
- warning("ECAT7DynamicDiscretisedDensityOutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
- this->type_of_numbers = supported_type_of_numbers;
- }
+ {
+ if (warn)
+ warning("ECAT7DynamicDiscretisedDensityOutputFileFormat: output type of numbers is currently fixed to short (2 byte "
+ "signed integers)\n");
+ this->type_of_numbers = supported_type_of_numbers;
+ }
else
this->type_of_numbers = new_type;
return this->type_of_numbers;
}
-ByteOrder
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-set_byte_order(const ByteOrder& new_byte_order, const bool warn)
+ByteOrder
+ECAT7DynamicDiscretisedDensityOutputFileFormat::set_byte_order(const ByteOrder& new_byte_order, const bool warn)
{
if (new_byte_order != ByteOrder::big_endian)
- {
- if (warn)
- warning("ECAT7DynamicDiscretisedDensityOutputFileFormat: byte_order is currently fixed to big-endian\n");
- this->file_byte_order = ByteOrder::big_endian;
- }
+ {
+ if (warn)
+ warning("ECAT7DynamicDiscretisedDensityOutputFileFormat: byte_order is currently fixed to big-endian\n");
+ this->file_byte_order = ByteOrder::big_endian;
+ }
else
this->file_byte_order = new_byte_order;
return this->file_byte_order;
}
-Succeeded
-ECAT7DynamicDiscretisedDensityOutputFileFormat::
-actual_write_to_file(std::string& filename,
- const DynamicDiscretisedDensity& dynamic_density) const
+Succeeded
+ECAT7DynamicDiscretisedDensityOutputFileFormat::actual_write_to_file(std::string& filename,
+ const DynamicDiscretisedDensity& dynamic_density) const
{
add_extension(filename, ".img");
- return
- dynamic_density.write_to_ecat7(filename);
+ return dynamic_density.write_to_ecat7(filename);
}
-
-//template class ECAT7DynamicDiscretisedDensityOutputFileFormat;
+// template class ECAT7DynamicDiscretisedDensityOutputFileFormat;
END_NAMESPACE_ECAT7
END_NAMESPACE_ECAT
diff --git a/src/IO/ECAT7OutputFileFormat.cxx b/src/IO/ECAT7OutputFileFormat.cxx
index defdc3238..28e13e3c0 100644
--- a/src/IO/ECAT7OutputFileFormat.cxx
+++ b/src/IO/ECAT7OutputFileFormat.cxx
@@ -27,21 +27,17 @@ START_NAMESPACE_STIR
START_NAMESPACE_ECAT
START_NAMESPACE_ECAT7
-const char * const
-ECAT7OutputFileFormat::registered_name = "ECAT7";
+const char* const ECAT7OutputFileFormat::registered_name = "ECAT7";
-ECAT7OutputFileFormat::
-ECAT7OutputFileFormat(const NumericType& type,
- const ByteOrder& byte_order)
+ECAT7OutputFileFormat::ECAT7OutputFileFormat(const NumericType& type, const ByteOrder& byte_order)
{
base_type::set_defaults();
set_type_of_numbers(type);
set_byte_order(byte_order);
}
-void
-ECAT7OutputFileFormat::
-initialise_keymap()
+void
+ECAT7OutputFileFormat::initialise_keymap()
{
parser.add_start_key("ECAT7 Output File Format Parameters");
parser.add_stop_key("End ECAT7 Output File Format Parameters");
@@ -49,9 +45,8 @@ initialise_keymap()
base_type::initialise_keymap();
}
-void
-ECAT7OutputFileFormat::
-set_defaults()
+void
+ECAT7OutputFileFormat::set_defaults()
{
default_scanner_name = "ECAT 962";
base_type::set_defaults();
@@ -59,88 +54,76 @@ set_defaults()
type_of_numbers = NumericType::SHORT;
set_key_values();
-
}
bool
-ECAT7OutputFileFormat::
-post_processing()
+ECAT7OutputFileFormat::post_processing()
{
if (base_type::post_processing())
return true;
shared_ptr scanner_ptr(Scanner::get_scanner_from_name(default_scanner_name));
- if (find_ECAT_system_type(*scanner_ptr)==0)
+ if (find_ECAT_system_type(*scanner_ptr) == 0)
{
- warning("ECAT7OutputFileFormat: default_scanner_name %s is not supported\n",
- default_scanner_name.c_str());
+ warning("ECAT7OutputFileFormat: default_scanner_name %s is not supported\n", default_scanner_name.c_str());
return true;
}
return false;
}
-NumericType
-ECAT7OutputFileFormat::
-set_type_of_numbers(const NumericType& new_type, const bool warn)
+NumericType
+ECAT7OutputFileFormat::set_type_of_numbers(const NumericType& new_type, const bool warn)
{
- const NumericType supported_type_of_numbers =
- NumericType("signed integer", 2);
+ const NumericType supported_type_of_numbers = NumericType("signed integer", 2);
if (new_type != supported_type_of_numbers)
- {
- if (warn)
- warning("ECAT7OutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
- type_of_numbers = supported_type_of_numbers;
- }
+ {
+ if (warn)
+ warning("ECAT7OutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
+ type_of_numbers = supported_type_of_numbers;
+ }
else
type_of_numbers = new_type;
return type_of_numbers;
-
}
-ByteOrder
-ECAT7OutputFileFormat::
-set_byte_order(const ByteOrder& new_byte_order, const bool warn)
+ByteOrder
+ECAT7OutputFileFormat::set_byte_order(const ByteOrder& new_byte_order, const bool warn)
{
if (new_byte_order != ByteOrder::big_endian)
- {
- if (warn)
- warning("ECAT7OutputFileFormat: byte_order is currently fixed to big-endian\n");
- file_byte_order = ByteOrder::big_endian;
- }
+ {
+ if (warn)
+ warning("ECAT7OutputFileFormat: byte_order is currently fixed to big-endian\n");
+ file_byte_order = ByteOrder::big_endian;
+ }
else
file_byte_order = new_byte_order;
return file_byte_order;
}
-Succeeded
-ECAT7OutputFileFormat::
-actual_write_to_file(std::string& filename,
- const DiscretisedDensity<3,float>& density) const
+Succeeded
+ECAT7OutputFileFormat::actual_write_to_file(std::string& filename, const DiscretisedDensity<3, float>& density) const
{
shared_ptr scanner_ptr(Scanner::get_scanner_from_name(default_scanner_name));
-
+
add_extension(filename, ".img");
Main_header mhead;
make_ECAT7_main_header(mhead, *scanner_ptr, "", density);
mhead.num_frames = 1;
- mhead.acquisition_type =
- mhead.num_frames>1 ? DynamicEmission : StaticEmission;
-
- MatrixFile* mptr= matrix_create (filename.c_str(), MAT_CREATE, &mhead);
+ mhead.acquisition_type = mhead.num_frames > 1 ? DynamicEmission : StaticEmission;
+
+ MatrixFile* mptr = matrix_create(filename.c_str(), MAT_CREATE, &mhead);
if (mptr == 0)
- {
- warning("ECAT7OutputFileFormat::write_to_file: error opening output file %s\n",
- filename.c_str());
- return Succeeded::no;
- }
-
- Succeeded success =
- DiscretisedDensity_to_ECAT7(mptr, density, 1 /*frame_num*/);
+ {
+ warning("ECAT7OutputFileFormat::write_to_file: error opening output file %s\n", filename.c_str());
+ return Succeeded::no;
+ }
+
+ Succeeded success = DiscretisedDensity_to_ECAT7(mptr, density, 1 /*frame_num*/);
matrix_close(mptr);
-
+
return success;
}
diff --git a/src/IO/ECAT7ParametricDensityOutputFileFormat.cxx b/src/IO/ECAT7ParametricDensityOutputFileFormat.cxx
index 9da8ca8a8..b731d5d9c 100644
--- a/src/IO/ECAT7ParametricDensityOutputFileFormat.cxx
+++ b/src/IO/ECAT7ParametricDensityOutputFileFormat.cxx
@@ -28,13 +28,11 @@ START_NAMESPACE_ECAT
START_NAMESPACE_ECAT7
template
-const char * const
-ECAT7ParametricDensityOutputFileFormat::registered_name = "ECAT7";
+const char* const ECAT7ParametricDensityOutputFileFormat::registered_name = "ECAT7";
template
-ECAT7ParametricDensityOutputFileFormat::
-ECAT7ParametricDensityOutputFileFormat(const NumericType& type,
- const ByteOrder& byte_order)
+ECAT7ParametricDensityOutputFileFormat::ECAT7ParametricDensityOutputFileFormat(const NumericType& type,
+ const ByteOrder& byte_order)
{
this->set_defaults();
this->set_type_of_numbers(type);
@@ -42,9 +40,8 @@ ECAT7ParametricDensityOutputFileFormat(const NumericType& type,
}
template
-void
-ECAT7ParametricDensityOutputFileFormat::
-initialise_keymap()
+void
+ECAT7ParametricDensityOutputFileFormat::initialise_keymap()
{
this->parser.add_start_key("ECAT7 Output File Format Parameters");
this->parser.add_stop_key("End ECAT7 Output File Format Parameters");
@@ -53,9 +50,8 @@ initialise_keymap()
}
template
-void
-ECAT7ParametricDensityOutputFileFormat::
-set_defaults()
+void
+ECAT7ParametricDensityOutputFileFormat::set_defaults()
{
this->default_scanner_name = "ECAT 962";
base_type::set_defaults();
@@ -63,24 +59,21 @@ set_defaults()
this->type_of_numbers = NumericType::SHORT;
this->set_key_values();
-
}
template
bool
-ECAT7ParametricDensityOutputFileFormat::
-post_processing()
+ECAT7ParametricDensityOutputFileFormat::post_processing()
{
if (base_type::post_processing())
return true;
- shared_ptr scanner_ptr(
- Scanner::get_scanner_from_name(this->default_scanner_name));
+ shared_ptr scanner_ptr(Scanner::get_scanner_from_name(this->default_scanner_name));
- if (find_ECAT_system_type(*scanner_ptr)==0)
+ if (find_ECAT_system_type(*scanner_ptr) == 0)
{
warning("ECAT7ParametricDensityOutputFileFormat: default_scanner_name %s is not supported\n",
- this->default_scanner_name.c_str());
+ this->default_scanner_name.c_str());
return true;
}
@@ -88,54 +81,47 @@ post_processing()
}
template
-NumericType
-ECAT7ParametricDensityOutputFileFormat::
-set_type_of_numbers(const NumericType& new_type, const bool warn)
+NumericType
+ECAT7ParametricDensityOutputFileFormat::set_type_of_numbers(const NumericType& new_type, const bool warn)
{
- const NumericType supported_type_of_numbers =
- NumericType("signed integer", 2);
+ const NumericType supported_type_of_numbers = NumericType("signed integer", 2);
if (new_type != supported_type_of_numbers)
- {
- if (warn)
- warning("ECAT7ParametricDensityOutputFileFormat: output type of numbers is currently fixed to short (2 byte signed integers)\n");
- this->type_of_numbers = supported_type_of_numbers;
- }
+ {
+ if (warn)
+ warning("ECAT7ParametricDensityOutputFileFormat: output type of numbers is currently fixed to short (2 byte signed "
+ "integers)\n");
+ this->type_of_numbers = supported_type_of_numbers;
+ }
else
this->type_of_numbers = new_type;
return this->type_of_numbers;
-
}
template
-ByteOrder
-ECAT7ParametricDensityOutputFileFormat::
-set_byte_order(const ByteOrder& new_byte_order, const bool warn)
+ByteOrder
+ECAT7ParametricDensityOutputFileFormat::set_byte_order(const ByteOrder& new_byte_order, const bool warn)
{
if (new_byte_order != ByteOrder::big_endian)
- {
- if (warn)
- warning("ECAT7ParametricDensityOutputFileFormat: byte_order is currently fixed to big-endian\n");
- this->file_byte_order = ByteOrder::big_endian;
- }
+ {
+ if (warn)
+ warning("ECAT7ParametricDensityOutputFileFormat: byte_order is currently fixed to big-endian\n");
+ this->file_byte_order = ByteOrder::big_endian;
+ }
else
this->file_byte_order = new_byte_order;
return this->file_byte_order;
}
template
-Succeeded
-ECAT7ParametricDensityOutputFileFormat::
-actual_write_to_file(std::string& filename,
- const ParametricDiscretisedDensity& parametric_density) const
+Succeeded
+ECAT7ParametricDensityOutputFileFormat::actual_write_to_file(
+ std::string& filename, const ParametricDiscretisedDensity& parametric_density) const
{
- shared_ptr scanner_ptr(
- Scanner::get_scanner_from_name(this->default_scanner_name));
-
+ shared_ptr scanner_ptr(Scanner::get_scanner_from_name(this->default_scanner_name));
+
add_extension(filename, ".img");
- typedef
- typename ParametricDiscretisedDensity::SingleDiscretisedDensityType
- SingleDensityType;
+ typedef typename ParametricDiscretisedDensity::SingleDiscretisedDensityType SingleDensityType;
// somewhat naughty trick to get elemT of DiscretisedDensityT
typedef typename DiscretisedDensityT::full_value_type KinParsT;
@@ -146,34 +132,31 @@ actual_write_to_file(std::string& filename,
mhead.num_frames = KinParsT::size();
mhead.acquisition_type = DynamicEmission;
-
- MatrixFile* mptr= matrix_create (filename.c_str(), MAT_CREATE, &mhead);
+
+ MatrixFile* mptr = matrix_create(filename.c_str(), MAT_CREATE, &mhead);
if (mptr == 0)
- {
- warning("ECAT7ParametricDensityOutputFileFormat::write_to_file: error opening output file %s\n",
- filename.c_str());
- return Succeeded::no;
- }
-
- for (unsigned int f=1; f<=KinParsT::size(); ++f)
- {
- parametric_density.construct_single_density(single_density,f);
- Succeeded success =
- DiscretisedDensity_to_ECAT7(mptr, single_density, f);
+ {
+ warning("ECAT7ParametricDensityOutputFileFormat::write_to_file: error opening output file %s\n", filename.c_str());
+ return Succeeded::no;
+ }
+
+ for (unsigned int f = 1; f <= KinParsT::size(); ++f)
+ {
+ parametric_density.construct_single_density(single_density, f);
+ Succeeded success = DiscretisedDensity_to_ECAT7(mptr, single_density, f);
if (success == Succeeded::no)
- {
- matrix_close(mptr);
- return success;
- }
+ {
+ matrix_close(mptr);
+ return success;
+ }
}
matrix_close(mptr);
-
+
return Succeeded::yes;
}
-template class ECAT7ParametricDensityOutputFileFormat;
+template class ECAT7ParametricDensityOutputFileFormat;
END_NAMESPACE_ECAT7
END_NAMESPACE_ECAT
END_NAMESPACE_STIR
-
diff --git a/src/IO/GEHDF5Wrapper.cxx b/src/IO/GEHDF5Wrapper.cxx
index 130ed98ec..45978ba2e 100644
--- a/src/IO/GEHDF5Wrapper.cxx
+++ b/src/IO/GEHDF5Wrapper.cxx
@@ -32,15 +32,16 @@
#include "stir/error.h"
#include
-
START_NAMESPACE_STIR
-namespace GE {
-namespace RDF_HDF5 {
+namespace GE
+{
+namespace RDF_HDF5
+{
std::uint32_t
GEHDF5Wrapper::read_dataset_uint32(const std::string& dataset_name)
{
- std::uint32_t tmp=0U;
+ std::uint32_t tmp = 0U;
H5::DataSet dataset = file.openDataSet(dataset_name);
dataset.read(&tmp, H5::PredType::NATIVE_UINT32);
return tmp;
@@ -49,14 +50,14 @@ GEHDF5Wrapper::read_dataset_uint32(const std::string& dataset_name)
std::int32_t
GEHDF5Wrapper::read_dataset_int32(const std::string& dataset_name)
{
- std::int32_t tmp=0;
+ std::int32_t tmp = 0;
H5::DataSet dataset = file.openDataSet(dataset_name);
dataset.read(&tmp, H5::PredType::NATIVE_INT32);
return tmp;
}
-
-static float read_float(const H5::H5File& file, const std::string& dataset)
+static float
+read_float(const H5::H5File& file, const std::string& dataset)
{
float tmp = 0.F;
H5::DataSet ds = file.openDataSet(dataset.c_str());
@@ -64,532 +65,563 @@ static float read_float(const H5::H5File& file, const std::string& dataset)
return tmp;
}
-bool GEHDF5Wrapper::check_GE_signature(const std::string& filename)
+bool
+GEHDF5Wrapper::check_GE_signature(const std::string& filename)
{
- try
+ try
{
- if(!H5::H5File::isHdf5(filename))
+ if (!H5::H5File::isHdf5(filename))
return false;
H5::H5File file;
- file.openFile( filename, H5F_ACC_RDONLY );
+ file.openFile(filename, H5F_ACC_RDONLY);
return (check_GE_signature(file));
}
- catch(...)
+ catch (...)
{
- return false;
- }
+ return false;
+ }
}
-bool GEHDF5Wrapper::check_GE_signature(H5::H5File& file)
+bool
+GEHDF5Wrapper::check_GE_signature(H5::H5File& file)
{
- if(file.getId() == -1)
- error("File is not open. Aborting");
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+
+ H5::StrType vlst(
+ 0,
+ 37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
+ std::string read_str_manufacturer;
- H5::StrType vlst(0,37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
- std::string read_str_manufacturer;
+ H5::DataSet dataset2 = file.openDataSet("/HeaderData/ExamData/manufacturer");
+ dataset2.read(read_str_manufacturer, vlst);
- H5::DataSet dataset2= file.openDataSet("/HeaderData/ExamData/manufacturer");
- dataset2.read(read_str_manufacturer,vlst);
-
- if(read_str_manufacturer == "GE MEDICAL SYSTEMS")
+ if (read_str_manufacturer == "GE MEDICAL SYSTEMS")
{
- return true;
+ return true;
}
- return false;
+ return false;
}
bool
GEHDF5Wrapper::is_list_file() const
{
- // have we already checked this? (note: initially set to false in check_file())
- if(is_list)
- return true;
-
- if(file.getId() == -1)
- error("File is not open. Aborting");
-
- // All RDF files shoudl have this DataSet
- H5::DataSet dataset = file.openDataSet("/HeaderData/RDFConfiguration/isListFile");
- std::uint32_t is_list_file;
- dataset.read(&is_list_file, H5::PredType::NATIVE_UINT32);
- return is_list_file;
-
+ // have we already checked this? (note: initially set to false in check_file())
+ if (is_list)
+ return true;
+
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+
+ // All RDF files shoudl have this DataSet
+ H5::DataSet dataset = file.openDataSet("/HeaderData/RDFConfiguration/isListFile");
+ std::uint32_t is_list_file;
+ dataset.read(&is_list_file, H5::PredType::NATIVE_UINT32);
+ return is_list_file;
}
-// Checks if file is a sino file.
+// Checks if file is a sino file.
// AB todo: only valid for RDF9 (until they tell us otherwise)
bool
GEHDF5Wrapper::is_sino_file() const
{
- if(is_sino)
- return true;
- if(file.getId() == -1)
- error("File is not open. Aborting");
-
- // If this Group exists, its a sino file.
- // huh, no C++ way to do this without try catch. https://stackoverflow.com/questions/35668056/test-group-existence-in-hdf5-c
- return H5Lexists( file.getId(), "/SegmentData/Segment2", H5P_DEFAULT ) > 0;
+ if (is_sino)
+ return true;
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+
+ // If this Group exists, its a sino file.
+ // huh, no C++ way to do this without try catch. https://stackoverflow.com/questions/35668056/test-group-existence-in-hdf5-c
+ return H5Lexists(file.getId(), "/SegmentData/Segment2", H5P_DEFAULT) > 0;
}
bool
GEHDF5Wrapper::is_geo_file() const
{
- // Apparently the norm file has the geo file inside. This means that the geo file is superfluous.
- // For now, lets just make this fucntion is_geo_or_norm_file(). Perhaps later versions will no be like this and this we need to change this function.
-
- if(is_geo || is_norm)
- return true;
- if(file.getId() == -1)
- error("File is not open. Aborting");
- return H5Lexists( file.getId(), "/SegmentData/Segment4/3D_Norm_Correction/slice1", H5P_DEFAULT ) > 0;
- // AB if you want to make sure geo is definetly geo, and not geo or norm, add to the avobe line: && !H5Lexists( file.getId(), "/3DCrystalEfficiency/crystalEfficiency", H5P_DEFAULT );
+ // Apparently the norm file has the geo file inside. This means that the geo file is superfluous.
+ // For now, lets just make this fucntion is_geo_or_norm_file(). Perhaps later versions will no be like this and this we need to
+ // change this function.
+
+ if (is_geo || is_norm)
+ return true;
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+ return H5Lexists(file.getId(), "/SegmentData/Segment4/3D_Norm_Correction/slice1", H5P_DEFAULT) > 0;
+ // AB if you want to make sure geo is definetly geo, and not geo or norm, add to the avobe line: && !H5Lexists( file.getId(),
+ // "/3DCrystalEfficiency/crystalEfficiency", H5P_DEFAULT );
}
bool
GEHDF5Wrapper::is_norm_file() const
{
- if(is_norm)
- return true;
- if(file.getId() == -1)
- error("File is not open. Aborting");
+ if (is_norm)
+ return true;
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
- return H5Lexists( file.getId(), "/3DCrystalEfficiency/crystalEfficiency", H5P_DEFAULT ) > 0;
+ return H5Lexists(file.getId(), "/3DCrystalEfficiency/crystalEfficiency", H5P_DEFAULT) > 0;
}
GEHDF5Wrapper::GEHDF5Wrapper()
{
- // Not much.
+ // Not much.
}
GEHDF5Wrapper::GEHDF5Wrapper(const std::string& filename)
{
- if(open(filename) == Succeeded::no)
- error("GEHDF5Wrapper: Error opening HDF5 file. Abort.");
+ if (open(filename) == Succeeded::no)
+ error("GEHDF5Wrapper: Error opening HDF5 file. Abort.");
}
unsigned int
GEHDF5Wrapper::check_geo_type()
{
- if(!is_geo )
- error("Not a geo file. Aborting");
- if(file.getId() == -1)
- error("File is not open. Aborting");
-
- unsigned int geo_type=0;
- H5::DataSet str_geo_size = file.openDataSet("/HeaderData/Sorter/Segment4/dimension3Size");
- str_geo_size.read(&geo_type, H5::PredType::NATIVE_UINT32);
- if (geo_type>1)
- geo_type=3;
- else
- geo_type=2;
- return geo_type;
+ if (!is_geo)
+ error("Not a geo file. Aborting");
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+
+ unsigned int geo_type = 0;
+ H5::DataSet str_geo_size = file.openDataSet("/HeaderData/Sorter/Segment4/dimension3Size");
+ str_geo_size.read(&geo_type, H5::PredType::NATIVE_UINT32);
+ if (geo_type > 1)
+ geo_type = 3;
+ else
+ geo_type = 2;
+ return geo_type;
}
// Function that error checks the input file and sets flags for the correct formats. GEHDF5Wrapper.file must be already open.
-// AB todo: this file is only valid for RDF 9.
-Succeeded
+// AB todo: this file is only valid for RDF 9.
+Succeeded
GEHDF5Wrapper::check_file()
{
- if(file.getId()==-1)
- error("File is not open. Aborting");
- if(!check_GE_signature(file))
- error("File is HDF5 but not GE data. Aborting");
-
- //AB: We are openign a new file. The same class should not be used twice, but lets make sure that if it happens, we reseted the file identifiers.
- is_list=false; is_norm=false;is_geo=false;is_sino=false;
-
- //AB Find out the RDF version of the file.
- H5::DataSet str_file_version = file.openDataSet("/HeaderData/RDFConfiguration/fileVersion/majorVersion");
- str_file_version.read(&rdf_ver, H5::PredType::NATIVE_UINT32);
-
- //AB Lets just error for now.
- if (rdf_ver!=9)
- error("Only RDF version 9 supported. Aborting");
-
- if(is_list_file())
+ if (file.getId() == -1)
+ error("File is not open. Aborting");
+ if (!check_GE_signature(file))
+ error("File is HDF5 but not GE data. Aborting");
+
+ // AB: We are openign a new file. The same class should not be used twice, but lets make sure that if it happens, we reseted the
+ // file identifiers.
+ is_list = false;
+ is_norm = false;
+ is_geo = false;
+ is_sino = false;
+
+ // AB Find out the RDF version of the file.
+ H5::DataSet str_file_version = file.openDataSet("/HeaderData/RDFConfiguration/fileVersion/majorVersion");
+ str_file_version.read(&rdf_ver, H5::PredType::NATIVE_UINT32);
+
+ // AB Lets just error for now.
+ if (rdf_ver != 9)
+ error("Only RDF version 9 supported. Aborting");
+
+ if (is_list_file())
{
- is_list = true;
- // AB Now lets check all things that are required
- if (rdf_ver == 9) //AB todo: is this valid for 10?
+ is_list = true;
+ // AB Now lets check all things that are required
+ if (rdf_ver == 9) // AB todo: is this valid for 10?
{
- // Check 1: Is the file compressed?
- unsigned int is_compressed;
- H5::DataSet str_file_version = file.openDataSet("/HeaderData/ListHeader/isListCompressed");
- str_file_version.read(&is_compressed, H5::PredType::NATIVE_UINT32);
- if (is_compressed)
- error("The RDF9 Listmode file is compressed, we won't be able to read it. Please uncompress it and retry. Aborting");
+ // Check 1: Is the file compressed?
+ unsigned int is_compressed;
+ H5::DataSet str_file_version = file.openDataSet("/HeaderData/ListHeader/isListCompressed");
+ str_file_version.read(&is_compressed, H5::PredType::NATIVE_UINT32);
+ if (is_compressed)
+ error("The RDF9 Listmode file is compressed, we won't be able to read it. Please uncompress it and retry. Aborting");
}
- return Succeeded::yes;
+ return Succeeded::yes;
}
- if(is_sino_file())
+ if (is_sino_file())
{
- is_sino = true;
- return Succeeded::yes;
+ is_sino = true;
+ return Succeeded::yes;
}
- if(is_norm_file())
+ if (is_norm_file())
{
- is_norm=true;
- is_geo =true; // in RFD9, if its norm, it is also geo (it contains it)
- // Check the type of geo file it contains.
- geo_dims = check_geo_type();
-
- return Succeeded::yes;
+ is_norm = true;
+ is_geo = true; // in RFD9, if its norm, it is also geo (it contains it)
+ // Check the type of geo file it contains.
+ geo_dims = check_geo_type();
+
+ return Succeeded::yes;
}
- if(is_geo_file())
+ if (is_geo_file())
{
- is_geo=true;
- geo_dims = check_geo_type();
-
- return Succeeded::yes;
+ is_geo = true;
+ geo_dims = check_geo_type();
+
+ return Succeeded::yes;
}
- // should not get here.
- return Succeeded::no;
+ // should not get here.
+ return Succeeded::no;
}
Succeeded
GEHDF5Wrapper::open(const std::string& filename)
{
- if(!file.isHdf5(filename))
- error("GEHDF5Wrapper: The input file is not HDF5! Abort.");
+ if (!file.isHdf5(filename))
+ error("GEHDF5Wrapper: The input file is not HDF5! Abort.");
- file.openFile(filename, H5F_ACC_RDONLY);
+ file.openFile(filename, H5F_ACC_RDONLY);
- //AB: check if the input file is valid, not only as a HDF5, also a valid PET file.
- check_file();
+ // AB: check if the input file is valid, not only as a HDF5, also a valid PET file.
+ check_file();
- initialise_exam_info();
- initialise_proj_data_info_from_HDF5();
+ initialise_exam_info();
+ initialise_proj_data_info_from_HDF5();
- // functions above will throw actually, so if we get here, it worked.
- return Succeeded::yes;
+ // functions above will throw actually, so if we get here, it worked.
+ return Succeeded::yes;
}
-shared_ptr GEHDF5Wrapper::get_scanner_from_HDF5()
+shared_ptr
+GEHDF5Wrapper::get_scanner_from_HDF5()
{
- std::string read_str_scanner;
- H5::StrType vlst(0,37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
-
- H5::DataSet dataset = file.openDataSet("/HeaderData/ExamData/scannerDesc");
- dataset.read(read_str_scanner,vlst);
-
- float effective_ring_diameter;
- int num_transaxial_blocks_per_bucket = 0;
- int num_axial_blocks_per_bucket = 0;
- int axial_blocks_per_unit = 0;
- int radial_blocks_per_unit = 0;
- int axial_units_per_module = 0;
- int radial_units_per_module = 0;
- int axial_modules_per_system = 0;
- int radial_modules_per_system = 0;
- int max_num_non_arccorrected_bins = 0;
- int num_transaxial_crystals_per_block = 0;
- int num_axial_crystals_per_block = 0 ;
- float detector_axial_size = 0.0;
- float intrinsic_tilt = 0.0;
- int num_detector_layers = 1;
-
- H5::DataSet str_effective_ring_diameter = file.openDataSet("/HeaderData/SystemGeometry/effectiveRingDiameter");
- H5::DataSet str_axial_blocks_per_module = file.openDataSet("/HeaderData/SystemGeometry/axialBlocksPerModule");
- H5::DataSet str_radial_blocks_per_module = file.openDataSet("/HeaderData/SystemGeometry/radialBlocksPerModule");
- H5::DataSet str_axial_blocks_per_unit = file.openDataSet("/HeaderData/SystemGeometry/axialBlocksPerUnit");
- H5::DataSet str_radial_blocks_per_unit = file.openDataSet("/HeaderData/SystemGeometry/radialBlocksPerUnit");
- H5::DataSet str_axial_units_per_module = file.openDataSet("/HeaderData/SystemGeometry/axialUnitsPerModule");
- H5::DataSet str_radial_units_per_module = file.openDataSet("/HeaderData/SystemGeometry/radialUnitsPerModule");
- H5::DataSet str_axial_modules_per_system = file.openDataSet("/HeaderData/SystemGeometry/axialModulesPerSystem");
- H5::DataSet str_radial_modules_per_system = file.openDataSet("/HeaderData/SystemGeometry/radialModulesPerSystem");
- //! \todo P.W: Find the crystal gaps and other info missing.
- H5::DataSet str_detector_axial_size = file.openDataSet("/HeaderData/SystemGeometry/detectorAxialSize");
- H5::DataSet str_intrinsic_tilt = file.openDataSet("/HeaderData/SystemGeometry/transaxial_crystal_0_offset");
-
- H5::DataSet str_max_number_of_non_arc_corrected_bins;
- // TODO RDF10, what happens here?
- if (rdf_ver == 9)
- { // Bug in RDF9 makes this dimension2Size instead of the expected dimension1Size
- if(is_sino_file())
- str_max_number_of_non_arc_corrected_bins = file.openDataSet("/HeaderData/Sorter/dimension2Size");
- else
- str_max_number_of_non_arc_corrected_bins = file.openDataSet("/HeaderData/Sorter/dimension1Size");
+ std::string read_str_scanner;
+ H5::StrType vlst(
+ 0,
+ 37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF
+
+ H5::DataSet dataset = file.openDataSet("/HeaderData/ExamData/scannerDesc");
+ dataset.read(read_str_scanner, vlst);
+
+ float effective_ring_diameter;
+ int num_transaxial_blocks_per_bucket = 0;
+ int num_axial_blocks_per_bucket = 0;
+ int axial_blocks_per_unit = 0;
+ int radial_blocks_per_unit = 0;
+ int axial_units_per_module = 0;
+ int radial_units_per_module = 0;
+ int axial_modules_per_system = 0;
+ int radial_modules_per_system = 0;
+ int max_num_non_arccorrected_bins = 0;
+ int num_transaxial_crystals_per_block = 0;
+ int num_axial_crystals_per_block = 0;
+ float detector_axial_size = 0.0;
+ float intrinsic_tilt = 0.0;
+ int num_detector_layers = 1;
+
+ H5::DataSet str_effective_ring_diameter = file.openDataSet("/HeaderData/SystemGeometry/effectiveRingDiameter");
+ H5::DataSet str_axial_blocks_per_module = file.openDataSet("/HeaderData/SystemGeometry/axialBlocksPerModule");
+ H5::DataSet str_radial_blocks_per_module = file.openDataSet("/HeaderData/SystemGeometry/radialBlocksPerModule");
+ H5::DataSet str_axial_blocks_per_unit = file.openDataSet("/HeaderData/SystemGeometry/axialBlocksPerUnit");
+ H5::DataSet str_radial_blocks_per_unit = file.openDataSet("/HeaderData/SystemGeometry/radialBlocksPerUnit");
+ H5::DataSet str_axial_units_per_module = file.openDataSet("/HeaderData/SystemGeometry/axialUnitsPerModule");
+ H5::DataSet str_radial_units_per_module = file.openDataSet("/HeaderData/SystemGeometry/radialUnitsPerModule");
+ H5::DataSet str_axial_modules_per_system = file.openDataSet("/HeaderData/SystemGeometry/axialModulesPerSystem");
+ H5::DataSet str_radial_modules_per_system = file.openDataSet("/HeaderData/SystemGeometry/radialModulesPerSystem");
+ //! \todo P.W: Find the crystal gaps and other info missing.
+ H5::DataSet str_detector_axial_size = file.openDataSet("/HeaderData/SystemGeometry/detectorAxialSize");
+ H5::DataSet str_intrinsic_tilt = file.openDataSet("/HeaderData/SystemGeometry/transaxial_crystal_0_offset");
+
+ H5::DataSet str_max_number_of_non_arc_corrected_bins;
+ // TODO RDF10, what happens here?
+ if (rdf_ver == 9)
+ { // Bug in RDF9 makes this dimension2Size instead of the expected dimension1Size
+ if (is_sino_file())
+ str_max_number_of_non_arc_corrected_bins = file.openDataSet("/HeaderData/Sorter/dimension2Size");
+ else
+ str_max_number_of_non_arc_corrected_bins = file.openDataSet("/HeaderData/Sorter/dimension1Size");
+ }
+ H5::DataSet str_axial_crystals_per_block = file.openDataSet("/HeaderData/SystemGeometry/axialCrystalsPerBlock");
+ H5::DataSet str_radial_crystals_per_block = file.openDataSet("/HeaderData/SystemGeometry/radialCrystalsPerBlock");
+ // Convert to numbers.
+
+ str_radial_blocks_per_module.read(&num_transaxial_blocks_per_bucket, H5::PredType::NATIVE_UINT32);
+ str_axial_blocks_per_module.read(&num_axial_blocks_per_bucket, H5::PredType::NATIVE_UINT32);
+ str_axial_blocks_per_unit.read(&axial_blocks_per_unit, H5::PredType::NATIVE_UINT32);
+ str_radial_blocks_per_unit.read(&radial_blocks_per_unit, H5::PredType::NATIVE_UINT32);
+ str_axial_units_per_module.read(&axial_units_per_module, H5::PredType::NATIVE_UINT32);
+ str_radial_units_per_module.read(&radial_units_per_module, H5::PredType::NATIVE_UINT32);
+ str_axial_modules_per_system.read(&axial_modules_per_system, H5::PredType::NATIVE_UINT32);
+ str_radial_modules_per_system.read(&radial_modules_per_system, H5::PredType::NATIVE_UINT32);
+ str_detector_axial_size.read(&detector_axial_size, H5::PredType::NATIVE_FLOAT);
+ str_intrinsic_tilt.read(&intrinsic_tilt, H5::PredType::NATIVE_FLOAT);
+ str_effective_ring_diameter.read(&effective_ring_diameter, H5::PredType::NATIVE_FLOAT);
+ str_max_number_of_non_arc_corrected_bins.read(&max_num_non_arccorrected_bins, H5::PredType::NATIVE_UINT32);
+ str_radial_crystals_per_block.read(&num_transaxial_crystals_per_block, H5::PredType::NATIVE_UINT32);
+ str_axial_crystals_per_block.read(&num_axial_crystals_per_block, H5::PredType::NATIVE_UINT32);
+
+ // TOF related
+ const float timingResolutionInPico = read_float(file, "/HeaderData/SystemGeometry/timingResolutionInPico");
+ const int posCoincidenceWindow = read_dataset_int32("/HeaderData/AcqParameters/EDCATParameters/posCoincidenceWindow");
+ const int negCoincidenceWindow = read_dataset_int32("/HeaderData/AcqParameters/EDCATParameters/negCoincidenceWindow");
+ const float coincTimingPrecisionInPico
+ = read_float(file, "/HeaderData/AcqParameters/EDCATParameters/coincTimingPrecision") * 1000; // in nanoSecs in file
+ const int num_tof_bins = posCoincidenceWindow + negCoincidenceWindow + 1;
+
+ int num_rings = num_axial_blocks_per_bucket * num_axial_crystals_per_block * axial_modules_per_system;
+ int num_detectors_per_ring = num_transaxial_blocks_per_bucket * num_transaxial_crystals_per_block * radial_modules_per_system;
+ float ring_spacing = detector_axial_size / num_rings;
+ // PW Bin Size, default num of arc corrected bins and inner ring radius not found in RDF header.
+ // They will be set from the default STIR values
+ shared_ptr scanner_sptr(Scanner::get_scanner_from_name(read_str_scanner));
+ if (is_null_ptr(scanner_sptr))
+ error("Scanner read from RDF file is " + read_str_scanner + ", but this is not supported yet");
+
+ scanner_sptr->set_num_detectors_per_ring(num_detectors_per_ring);
+ scanner_sptr->set_num_rings(num_rings);
+ if (!is_list_file())
+ scanner_sptr->set_max_num_non_arccorrected_bins(max_num_non_arccorrected_bins);
+ scanner_sptr->set_ring_spacing(ring_spacing);
+ scanner_sptr->set_intrinsic_azimuthal_tilt(intrinsic_tilt * _PI / 180);
+ scanner_sptr->set_num_axial_blocks_per_bucket(num_axial_blocks_per_bucket);
+ scanner_sptr->set_num_transaxial_blocks_per_bucket(num_transaxial_blocks_per_bucket);
+ scanner_sptr->set_num_axial_crystals_per_block(num_axial_crystals_per_block);
+ scanner_sptr->set_num_transaxial_crystals_per_block(num_transaxial_crystals_per_block);
+ scanner_sptr->set_num_detector_layers(num_detector_layers);
+ scanner_sptr->set_reference_energy(511.F);
+
+ if (fabs(scanner_sptr->get_effective_ring_radius() - effective_ring_diameter / 2) > .1F)
+ {
+ const float def_DOI = .0F;
+ warning("GEHDF5Wrapper: default STIR effective ring radius is " + std::to_string(scanner_sptr->get_effective_ring_radius())
+ + ", while RDF says " + std::to_string(effective_ring_diameter / 2)
+ + "\nWill adjust scanner info to fit with the RDF file using default average DOI of " + std::to_string(def_DOI)
+ + "mm");
+ scanner_sptr->set_inner_ring_radius((effective_ring_diameter / 2) - def_DOI);
+ scanner_sptr->set_average_depth_of_interaction(def_DOI);
+ }
+ if (timingResolutionInPico > 0 // Signa files seem to have zero in this field
+ && (fabs(scanner_sptr->get_timing_resolution() - timingResolutionInPico) > .1F))
+ {
+ warning("GEHDF5Wrapper: default STIR timing resolution is " + std::to_string(scanner_sptr->get_timing_resolution())
+ + ", while RDF says " + std::to_string(timingResolutionInPico)
+ + "\nWill adjust scanner info to fit with the RDF file");
+ scanner_sptr->set_timing_resolution(timingResolutionInPico);
+ }
+ if (fabs(scanner_sptr->get_size_of_timing_pos() - coincTimingPrecisionInPico) > .1F)
+ {
+ warning("GEHDF5Wrapper: default STIR size of (unmashed) TOF bins is "
+ + std::to_string(scanner_sptr->get_size_of_timing_pos()) + ", while RDF says "
+ + std::to_string(coincTimingPrecisionInPico) + "\nWill adjust scanner info to fit with the RDF file");
+ scanner_sptr->set_size_of_timing_poss(coincTimingPrecisionInPico);
+ }
+ if (std::abs(scanner_sptr->get_max_num_timing_poss() - num_tof_bins) > 0)
+ {
+ warning("GEHDF5Wrapper: default STIR number of (unmashed) TOF bins is "
+ + std::to_string(scanner_sptr->get_max_num_timing_poss()) + ", while RDF says " + std::to_string(num_tof_bins)
+ + "\nWill adjust scanner info to fit with the RDF file");
+ scanner_sptr->set_max_num_timing_poss(num_tof_bins);
+ }
+ if (scanner_sptr->get_default_bin_size() <= 0.F)
+ {
+ warning("GEHDF5Wrapper: default bin-size is not set. This will create trouble for FBP etc");
+ }
+ if (scanner_sptr->get_default_num_arccorrected_bins() <= 0)
+ {
+ warning("GEHDF5Wrapper: default num_arccorrected bins is not set. This will create trouble for FBP etc");
+ }
+ if (scanner_sptr->get_energy_resolution() <= 0)
+ {
+ warning("GEHDF5Wrapper: energy resolution is not set. This will create trouble for scatter estimation");
}
- H5::DataSet str_axial_crystals_per_block = file.openDataSet("/HeaderData/SystemGeometry/axialCrystalsPerBlock");
- H5::DataSet str_radial_crystals_per_block = file.openDataSet("/HeaderData/SystemGeometry/radialCrystalsPerBlock");
- // Convert to numbers.
-
- str_radial_blocks_per_module.read(&num_transaxial_blocks_per_bucket, H5::PredType::NATIVE_UINT32);
- str_axial_blocks_per_module.read(&num_axial_blocks_per_bucket, H5::PredType::NATIVE_UINT32);
- str_axial_blocks_per_unit.read(&axial_blocks_per_unit, H5::PredType::NATIVE_UINT32);
- str_radial_blocks_per_unit.read(&radial_blocks_per_unit, H5::PredType::NATIVE_UINT32);
- str_axial_units_per_module.read(&axial_units_per_module, H5::PredType::NATIVE_UINT32);
- str_radial_units_per_module.read(&radial_units_per_module, H5::PredType::NATIVE_UINT32);
- str_axial_modules_per_system.read(&axial_modules_per_system, H5::PredType::NATIVE_UINT32);
- str_radial_modules_per_system.read(&radial_modules_per_system, H5::PredType::NATIVE_UINT32);
- str_detector_axial_size.read(&detector_axial_size, H5::PredType::NATIVE_FLOAT);
- str_intrinsic_tilt.read(&intrinsic_tilt, H5::PredType::NATIVE_FLOAT);
- str_effective_ring_diameter.read(&effective_ring_diameter, H5::PredType::NATIVE_FLOAT);
- str_max_number_of_non_arc_corrected_bins.read(&max_num_non_arccorrected_bins, H5::PredType::NATIVE_UINT32);
- str_radial_crystals_per_block.read(&num_transaxial_crystals_per_block, H5::PredType::NATIVE_UINT32);
- str_axial_crystals_per_block.read(&num_axial_crystals_per_block, H5::PredType::NATIVE_UINT32);
-
- // TOF related
- const float timingResolutionInPico = read_float(file, "/HeaderData/SystemGeometry/timingResolutionInPico");
- const int posCoincidenceWindow = read_dataset_int32("/HeaderData/AcqParameters/EDCATParameters/posCoincidenceWindow");
- const int negCoincidenceWindow = read_dataset_int32("/HeaderData/AcqParameters/EDCATParameters/negCoincidenceWindow");
- const float coincTimingPrecisionInPico = read_float(file, "/HeaderData/AcqParameters/EDCATParameters/coincTimingPrecision") * 1000; // in nanoSecs in file
- const int num_tof_bins = posCoincidenceWindow + negCoincidenceWindow + 1;
-
- int num_rings = num_axial_blocks_per_bucket*num_axial_crystals_per_block*axial_modules_per_system;
- int num_detectors_per_ring = num_transaxial_blocks_per_bucket*num_transaxial_crystals_per_block*radial_modules_per_system;
- float ring_spacing = detector_axial_size/num_rings;
- //PW Bin Size, default num of arc corrected bins and inner ring radius not found in RDF header.
- // They will be set from the default STIR values
- shared_ptr scanner_sptr(Scanner::get_scanner_from_name(read_str_scanner));
- if (is_null_ptr(scanner_sptr))
- error("Scanner read from RDF file is " + read_str_scanner + ", but this is not supported yet");
-
- scanner_sptr->set_num_detectors_per_ring(num_detectors_per_ring);
- scanner_sptr->set_num_rings(num_rings);
- if (!is_list_file())
- scanner_sptr->set_max_num_non_arccorrected_bins(max_num_non_arccorrected_bins);
- scanner_sptr->set_ring_spacing(ring_spacing);
- scanner_sptr->set_intrinsic_azimuthal_tilt(intrinsic_tilt*_PI/180);
- scanner_sptr->set_num_axial_blocks_per_bucket(num_axial_blocks_per_bucket);
- scanner_sptr->set_num_transaxial_blocks_per_bucket(num_transaxial_blocks_per_bucket);
- scanner_sptr->set_num_axial_crystals_per_block(num_axial_crystals_per_block);
- scanner_sptr->set_num_transaxial_crystals_per_block(num_transaxial_crystals_per_block);
- scanner_sptr->set_num_detector_layers(num_detector_layers);
- scanner_sptr->set_reference_energy(511.F);
-
- if (fabs(scanner_sptr->get_effective_ring_radius() - effective_ring_diameter/2) > .1F)
- {
- const float def_DOI = .0F;
- warning("GEHDF5Wrapper: default STIR effective ring radius is "
- + std::to_string(scanner_sptr->get_effective_ring_radius())
- + ", while RDF says " + std::to_string(effective_ring_diameter/2)
- + "\nWill adjust scanner info to fit with the RDF file using default average DOI of "
- + std::to_string(def_DOI) + "mm");
- scanner_sptr->set_inner_ring_radius((effective_ring_diameter/2) - def_DOI);
- scanner_sptr->set_average_depth_of_interaction(def_DOI);
- }
- if (timingResolutionInPico > 0 // Signa files seem to have zero in this field
- && (fabs(scanner_sptr->get_timing_resolution() - timingResolutionInPico) > .1F))
- {
- warning("GEHDF5Wrapper: default STIR timing resolution is "
- + std::to_string(scanner_sptr->get_timing_resolution())
- + ", while RDF says " + std::to_string(timingResolutionInPico)
- + "\nWill adjust scanner info to fit with the RDF file");
- scanner_sptr->set_timing_resolution(timingResolutionInPico);
- }
- if (fabs(scanner_sptr->get_size_of_timing_pos() - coincTimingPrecisionInPico)> .1F)
- {
- warning("GEHDF5Wrapper: default STIR size of (unmashed) TOF bins is "
- + std::to_string(scanner_sptr->get_size_of_timing_pos())
- + ", while RDF says " + std::to_string(coincTimingPrecisionInPico)
- + "\nWill adjust scanner info to fit with the RDF file");
- scanner_sptr->set_size_of_timing_poss(coincTimingPrecisionInPico);
- }
- if (std::abs(scanner_sptr->get_max_num_timing_poss() - num_tof_bins) > 0)
- {
- warning("GEHDF5Wrapper: default STIR number of (unmashed) TOF bins is "
- + std::to_string(scanner_sptr->get_max_num_timing_poss())
- + ", while RDF says " + std::to_string(num_tof_bins)
- + "\nWill adjust scanner info to fit with the RDF file");
- scanner_sptr->set_max_num_timing_poss(num_tof_bins);
- }
- if (scanner_sptr->get_default_bin_size() <= 0.F)
- {
- warning("GEHDF5Wrapper: default bin-size is not set. This will create trouble for FBP etc");
- }
- if (scanner_sptr->get_default_num_arccorrected_bins() <= 0)
- {
- warning("GEHDF5Wrapper: default num_arccorrected bins is not set. This will create trouble for FBP etc");
- }
- if (scanner_sptr->get_energy_resolution() <= 0)
- {
- warning("GEHDF5Wrapper: energy resolution is not set. This will create trouble for scatter estimation");
- }
- return scanner_sptr;
+ return scanner_sptr;
}
-void GEHDF5Wrapper::initialise_proj_data_info_from_HDF5()
+void
+GEHDF5Wrapper::initialise_proj_data_info_from_HDF5()
{
shared_ptr scanner_sptr = get_scanner_from_HDF5();
- this->proj_data_info_sptr =
- ProjDataInfo::construct_proj_data_info(scanner_sptr,
- /*span*/ 2,
- /* max_delta*/ scanner_sptr->get_num_rings()-1,
- /* num_views */ scanner_sptr->get_num_detectors_per_ring()/2,
- /* num_tangential_poss */ scanner_sptr->get_max_num_non_arccorrected_bins(),
- /* arc_corrected */ false,
- this->is_list_file() ? 1 : 0 // TODO change when reading sinos as TOF
- );
- this->proj_data_info_sptr->
- set_bed_position_horizontal(this->read_dataset_int32("/HeaderData/AcqParameters/LandmarkParameters/absTableLongitude")/10.F); /* units in RDF are 0.1 mm */
- //this->proj_data_info_sptr->set_gantry_tilt(this->read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/gantryTilt")); /* units in RDF are 0.25 degrees, patient relative */
- this->proj_data_info_sptr->
- set_bed_position_vertical(this->read_dataset_int32("/HeaderData/AcqParameters/LandmarkParameters/tableElevation")/10.F); /* units in RDF are 0.1 mm */
+ this->proj_data_info_sptr
+ = ProjDataInfo::construct_proj_data_info(scanner_sptr,
+ /*span*/ 2,
+ /* max_delta*/ scanner_sptr->get_num_rings() - 1,
+ /* num_views */ scanner_sptr->get_num_detectors_per_ring() / 2,
+ /* num_tangential_poss */ scanner_sptr->get_max_num_non_arccorrected_bins(),
+ /* arc_corrected */ false,
+ this->is_list_file() ? 1 : 0 // TODO change when reading sinos as TOF
+ );
+ this->proj_data_info_sptr->set_bed_position_horizontal(
+ this->read_dataset_int32("/HeaderData/AcqParameters/LandmarkParameters/absTableLongitude")
+ / 10.F); /* units in RDF are 0.1 mm */
+ // this->proj_data_info_sptr->set_gantry_tilt(this->read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/gantryTilt"));
+ // /* units in RDF are 0.25 degrees, patient relative */
+ this->proj_data_info_sptr->set_bed_position_vertical(
+ this->read_dataset_int32("/HeaderData/AcqParameters/LandmarkParameters/tableElevation")
+ / 10.F); /* units in RDF are 0.1 mm */
}
-unsigned int GEHDF5Wrapper::get_num_singles_samples()
+unsigned int
+GEHDF5Wrapper::get_num_singles_samples()
{
- return m_num_singles_samples;
+ return m_num_singles_samples;
}
-
-void GEHDF5Wrapper::initialise_exam_info()
+void
+GEHDF5Wrapper::initialise_exam_info()
{
- this->exam_info_sptr.reset(new ExamInfo());
- this->exam_info_sptr->imaging_modality = ImagingModality(ImagingModality::PT);
- {
- const std::uint32_t patientEntry = read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/patientEntry");
- const std::uint32_t patientPosition = read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/patientPosition");
- PatientPosition::OrientationValue orientation;
- PatientPosition::RotationValue rotation;
- switch (patientEntry)
- {
- case AcqPatientEntries::ACQ_HEAD_FIRST: orientation = PatientPosition::OrientationValue::head_in; break;
- case AcqPatientEntries::ACQ_FEET_FIRST: orientation = PatientPosition::OrientationValue::feet_in; break;
- default: orientation = PatientPosition::OrientationValue::unknown_orientation;
- }
- switch (patientPosition)
- {
- case AcqPatientPositions::ACQ_SUPINE: rotation = PatientPosition::RotationValue::supine; break;
- case AcqPatientPositions::ACQ_PRONE: rotation = PatientPosition::RotationValue::prone; break;
- case AcqPatientPositions::ACQ_LEFT_DECUB: rotation = PatientPosition::RotationValue::left; break;
- case AcqPatientPositions::ACQ_RIGHT_DECUB: rotation = PatientPosition::RotationValue::right; break;
- default: rotation = PatientPosition::RotationValue::unknown_rotation; break;
- }
- exam_info_sptr->patient_position = PatientPosition(orientation, rotation);
- }
- // PW Get the high and low energy threshold values from HDF5 header.
- unsigned int low_energy_thres = 0;
- unsigned int high_energy_thres = 0;
-
- H5::DataSet str_low_energy_thres = file.openDataSet("/HeaderData/AcqParameters/EDCATParameters/lower_energy_limit");
- H5::DataSet str_high_energy_thres = file.openDataSet("/HeaderData/AcqParameters/EDCATParameters/upper_energy_limit");
-
- str_low_energy_thres.read(&low_energy_thres, H5::PredType::NATIVE_UINT32);
- str_high_energy_thres.read(&high_energy_thres, H5::PredType::NATIVE_UINT32);
-
- // PW Set these values in exam_info_sptr.
- exam_info_sptr->set_high_energy_thres(static_cast(high_energy_thres));
- exam_info_sptr->set_low_energy_thres(static_cast(low_energy_thres));
-
- // read time since 1970
- std::uint32_t scanStartTime = 0;
- H5::DataSet scanStartTime_dataset = file.openDataSet("/HeaderData/AcqStats/scanStartTime");
- scanStartTime_dataset.read(&scanStartTime, H5::PredType::NATIVE_UINT32);
- exam_info_sptr->start_time_in_secs_since_1970=double(scanStartTime);
-
- // get time frame
- std::uint32_t frameStartTime = 0;
- std::uint32_t frameDuration = 0;
- H5::DataSet frameStartTime_dataset = file.openDataSet("/HeaderData/AcqStats/frameStartTime");
- H5::DataSet frameDuration_dataset = file.openDataSet("/HeaderData/AcqStats/frameDuration");
-
- frameStartTime_dataset.read(&frameStartTime, H5::PredType::NATIVE_UINT32);
- frameDuration_dataset.read(&frameDuration, H5::PredType::NATIVE_UINT32);
- const double frame_start_time = double(frameStartTime - scanStartTime);
-
- std::vector >tf{{frame_start_time,frame_start_time+frameDuration/1000}};
-
- TimeFrameDefinitions tm(tf);
- exam_info_sptr->set_time_frame_definitions(tm);
-
- // radionuclide
- {
- auto rn_name_ds = file.openDataSet("/HeaderData/ExamData/radionuclideName");
- H5::StrType str_type(rn_name_ds);
- std::string rn_name;
- rn_name_ds.read(rn_name, str_type);
- RadionuclideDB radionuclide_db;
- Radionuclide radionuclide = radionuclide_db.get_radionuclide(exam_info_sptr->imaging_modality,rn_name);
-
- const float positron_fraction = read_float(file, "/HeaderData/ExamData/positronFraction");
- const float half_life = read_float(file, "/HeaderData/ExamData/halfLife");
- if (radionuclide.get_half_life(false) < 0)
- radionuclide = Radionuclide(rn_name, 511.F, positron_fraction, half_life, exam_info_sptr->imaging_modality);
- exam_info_sptr->set_radionuclide(radionuclide);
- }
+ this->exam_info_sptr.reset(new ExamInfo());
+ this->exam_info_sptr->imaging_modality = ImagingModality(ImagingModality::PT);
+ {
+ const std::uint32_t patientEntry = read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/patientEntry");
+ const std::uint32_t patientPosition = read_dataset_uint32("/HeaderData/AcqParameters/LandmarkParameters/patientPosition");
+ PatientPosition::OrientationValue orientation;
+ PatientPosition::RotationValue rotation;
+ switch (patientEntry)
+ {
+ case AcqPatientEntries::ACQ_HEAD_FIRST:
+ orientation = PatientPosition::OrientationValue::head_in;
+ break;
+ case AcqPatientEntries::ACQ_FEET_FIRST:
+ orientation = PatientPosition::OrientationValue::feet_in;
+ break;
+ default:
+ orientation = PatientPosition::OrientationValue::unknown_orientation;
+ }
+ switch (patientPosition)
+ {
+ case AcqPatientPositions::ACQ_SUPINE:
+ rotation = PatientPosition::RotationValue::supine;
+ break;
+ case AcqPatientPositions::ACQ_PRONE:
+ rotation = PatientPosition::RotationValue::prone;
+ break;
+ case AcqPatientPositions::ACQ_LEFT_DECUB:
+ rotation = PatientPosition::RotationValue::left;
+ break;
+ case AcqPatientPositions::ACQ_RIGHT_DECUB:
+ rotation = PatientPosition::RotationValue::right;
+ break;
+ default:
+ rotation = PatientPosition::RotationValue::unknown_rotation;
+ break;
+ }
+ exam_info_sptr->patient_position = PatientPosition(orientation, rotation);
+ }
+ // PW Get the high and low energy threshold values from HDF5 header.
+ unsigned int low_energy_thres = 0;
+ unsigned int high_energy_thres = 0;
+
+ H5::DataSet str_low_energy_thres = file.openDataSet("/HeaderData/AcqParameters/EDCATParameters/lower_energy_limit");
+ H5::DataSet str_high_energy_thres = file.openDataSet("/HeaderData/AcqParameters/EDCATParameters/upper_energy_limit");
+
+ str_low_energy_thres.read(&low_energy_thres, H5::PredType::NATIVE_UINT32);
+ str_high_energy_thres.read(&high_energy_thres, H5::PredType::NATIVE_UINT32);
+
+ // PW Set these values in exam_info_sptr.
+ exam_info_sptr->set_high_energy_thres(static_cast(high_energy_thres));
+ exam_info_sptr->set_low_energy_thres(static_cast(low_energy_thres));
+
+ // read time since 1970
+ std::uint32_t scanStartTime = 0;
+ H5::DataSet scanStartTime_dataset = file.openDataSet("/HeaderData/AcqStats/scanStartTime");
+ scanStartTime_dataset.read(&scanStartTime, H5::PredType::NATIVE_UINT32);
+ exam_info_sptr->start_time_in_secs_since_1970 = double(scanStartTime);
+
+ // get time frame
+ std::uint32_t frameStartTime = 0;
+ std::uint32_t frameDuration = 0;
+ H5::DataSet frameStartTime_dataset = file.openDataSet("/HeaderData/AcqStats/frameStartTime");
+ H5::DataSet frameDuration_dataset = file.openDataSet("/HeaderData/AcqStats/frameDuration");
+
+ frameStartTime_dataset.read(&frameStartTime, H5::PredType::NATIVE_UINT32);
+ frameDuration_dataset.read(&frameDuration, H5::PredType::NATIVE_UINT32);
+ const double frame_start_time = double(frameStartTime - scanStartTime);
+
+ std::vector> tf{ { frame_start_time, frame_start_time + frameDuration / 1000 } };
+
+ TimeFrameDefinitions tm(tf);
+ exam_info_sptr->set_time_frame_definitions(tm);
+
+ // radionuclide
+ {
+ auto rn_name_ds = file.openDataSet("/HeaderData/ExamData/radionuclideName");
+ H5::StrType str_type(rn_name_ds);
+ std::string rn_name;
+ rn_name_ds.read(rn_name, str_type);
+ RadionuclideDB radionuclide_db;
+ Radionuclide radionuclide = radionuclide_db.get_radionuclide(exam_info_sptr->imaging_modality, rn_name);
+
+ const float positron_fraction = read_float(file, "/HeaderData/ExamData/positronFraction");
+ const float half_life = read_float(file, "/HeaderData/ExamData/halfLife");
+ if (radionuclide.get_half_life(false) < 0)
+ radionuclide = Radionuclide(rn_name, 511.F, positron_fraction, half_life, exam_info_sptr->imaging_modality);
+ exam_info_sptr->set_radionuclide(radionuclide);
+ }
}
-Succeeded GEHDF5Wrapper::initialise_listmode_data()
+Succeeded
+GEHDF5Wrapper::initialise_listmode_data()
{
- if(!is_list_file())
- error("The file provided is not listmode. Aborting");
+ if (!is_list_file())
+ error("The file provided is not listmode. Aborting");
- if(rdf_ver==9)
+ if (rdf_ver == 9)
{
- m_address = "/ListData/listData";
- // These values are not in the file are come from info shared by GE.
- m_size_of_record_signature = 6;
- m_max_size_of_record = 16;
-
- unsigned int num_time_slices = 0;
- H5::DataSet timeframe_dataspace = file.openDataSet("/HeaderData/SinglesHeader/numValidSamples");
- timeframe_dataspace.read(&num_time_slices, H5::PredType::NATIVE_UINT32);
- if(num_time_slices==0)
+ m_address = "/ListData/listData";
+ // These values are not in the file are come from info shared by GE.
+ m_size_of_record_signature = 6;
+ m_max_size_of_record = 16;
+
+ unsigned int num_time_slices = 0;
+ H5::DataSet timeframe_dataspace = file.openDataSet("/HeaderData/SinglesHeader/numValidSamples");
+ timeframe_dataspace.read(&num_time_slices, H5::PredType::NATIVE_UINT32);
+ if (num_time_slices == 0)
{
- error("Zero number of valid samples singles samples in data. Aborting");
+ error("Zero number of valid samples singles samples in data. Aborting");
}
- m_num_singles_samples=num_time_slices;
-
+ m_num_singles_samples = num_time_slices;
}
- else
- return Succeeded::no;
+ else
+ return Succeeded::no;
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
- m_dataspace = m_dataset_sptr->getSpace();
- m_dataset_list_Ndims = m_dataspace.getSimpleExtentNdims();
+ m_dataspace = m_dataset_sptr->getSpace();
+ m_dataset_list_Ndims = m_dataspace.getSimpleExtentNdims();
- // We allocate dims_out in the stack for efficiecy and safety but we need an error check just in case then
- if (m_dataset_list_Ndims>m_max_dataset_dims)
- error("Dataset dimensions ("+ std::to_string(m_dataset_list_Ndims) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
- hsize_t dims_out[m_max_dataset_dims];
+ // We allocate dims_out in the stack for efficiecy and safety but we need an error check just in case then
+ if (m_dataset_list_Ndims > m_max_dataset_dims)
+ error("Dataset dimensions (" + std::to_string(m_dataset_list_Ndims) + ") bigger than maximum of"
+ + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
+ hsize_t dims_out[m_max_dataset_dims];
- m_dataspace.getSimpleExtentDims( dims_out, NULL);
- m_list_size=dims_out[0];
+ m_dataspace.getSimpleExtentDims(dims_out, NULL);
+ m_list_size = dims_out[0];
- return Succeeded::yes;
+ return Succeeded::yes;
}
-Succeeded GEHDF5Wrapper::initialise_singles_data()
+Succeeded
+GEHDF5Wrapper::initialise_singles_data()
{
-
- if(!is_list_file() && !is_sino_file())
- error("The file provided is not listmode or sinogram data. Aborting");
-
- if(rdf_ver==9)
- {
- m_address = "/Singles/CrystalSingles/sample";
- // Get the DataSpace (metadata) corresponding to the DataSet that we want to read
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + "1")));
- m_dataspace = m_dataset_sptr->getSpace();
- // Create an array to host the size of the dimensions
- const int rank = m_dataspace.getSimpleExtentNdims();
- // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
- if (rank>m_max_dataset_dims)
- error("Dataset dimensions ("+ std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
- hsize_t dims[m_max_dataset_dims];
- // Read size of dimensions
- m_dataspace.getSimpleExtentDims( dims, NULL);
-
- m_NX_SUB = dims[0]; // hyperslab dimensions
- m_NY_SUB = dims[1];
- m_NZ_SUB = (rank==2)? 1 : dims[2]; // Signa has rank==2, but this stay shere just in case...
+ if (!is_list_file() && !is_sino_file())
+ error("The file provided is not listmode or sinogram data. Aborting");
- unsigned int num_time_slices = 0;
- H5::DataSet timeframe_dataspace = file.openDataSet("/HeaderData/SinglesHeader/numValidSamples");
- timeframe_dataspace.read(&num_time_slices, H5::PredType::NATIVE_UINT32);
- if(num_time_slices==0)
+ if (rdf_ver == 9)
+ {
+ m_address = "/Singles/CrystalSingles/sample";
+ // Get the DataSpace (metadata) corresponding to the DataSet that we want to read
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + "1")));
+ m_dataspace = m_dataset_sptr->getSpace();
+ // Create an array to host the size of the dimensions
+ const int rank = m_dataspace.getSimpleExtentNdims();
+ // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
+ if (rank > m_max_dataset_dims)
+ error("Dataset dimensions (" + std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims)
+ + ". This is unexpected, Aborting.");
+ hsize_t dims[m_max_dataset_dims];
+ // Read size of dimensions
+ m_dataspace.getSimpleExtentDims(dims, NULL);
+
+ m_NX_SUB = dims[0]; // hyperslab dimensions
+ m_NY_SUB = dims[1];
+ m_NZ_SUB = (rank == 2) ? 1 : dims[2]; // Signa has rank==2, but this stay shere just in case...
+
+ unsigned int num_time_slices = 0;
+ H5::DataSet timeframe_dataspace = file.openDataSet("/HeaderData/SinglesHeader/numValidSamples");
+ timeframe_dataspace.read(&num_time_slices, H5::PredType::NATIVE_UINT32);
+ if (num_time_slices == 0)
{
- error("Zero number of valid samples singles samples in data. Aborting");
+ error("Zero number of valid samples singles samples in data. Aborting");
}
- m_num_singles_samples=num_time_slices;
+ m_num_singles_samples = num_time_slices;
#if 0
m_NX = dims[0]; // output buffer dimensions
@@ -597,18 +629,19 @@ Succeeded GEHDF5Wrapper::initialise_singles_data()
m_NZ = (rank==2)? 1 : dims[2];
#endif
}
- else
- return Succeeded::no;
+ else
+ return Succeeded::no;
- return Succeeded::yes;
+ return Succeeded::yes;
}
-Succeeded GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
+Succeeded
+GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
{
- if(!is_sino_file())
- error("The file provided is not sinogram data. Aborting");
+ if (!is_sino_file())
+ error("The file provided is not sinogram data. Aborting");
- if(rdf_ver==9)
+ if (rdf_ver == 9)
{
// Is the file compressed?
unsigned int is_compressed;
@@ -618,28 +651,29 @@ Succeeded GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
error("The RDF9 file sinogram is compressed, we won't be able to read it. Please uncompress it and retry. Aborting");
}
- if(view_num == 0 || view_num > static_cast(this->get_scanner_sptr()->get_num_detectors_per_ring()/2))
- error("internal error in GE HDF5 code: view number "+ std::to_string(view_num) +" is incorrect");
+ if (view_num == 0 || view_num > static_cast(this->get_scanner_sptr()->get_num_detectors_per_ring() / 2))
+ error("internal error in GE HDF5 code: view number " + std::to_string(view_num) + " is incorrect");
- if(rdf_ver==9)
+ if (rdf_ver == 9)
{
- m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num);
-
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
- m_dataspace = m_dataset_sptr->getSpace();
- // Create an array to host the size of the dimensions
- const int rank = m_dataspace.getSimpleExtentNdims();
- // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
- if (rank>m_max_dataset_dims)
- error("Dataset dimensions ("+ std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
- hsize_t dims[m_max_dataset_dims];
- // Read size of dimensions
- m_dataspace.getSimpleExtentDims( dims, NULL);
-
- // AB for signa, these where [1981,27,357] and [45,448,357]
- m_NX_SUB = dims[0] ; // hyperslab dimensions
- m_NY_SUB = dims[1];
- m_NZ_SUB = dims[2];
+ m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num);
+
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+ m_dataspace = m_dataset_sptr->getSpace();
+ // Create an array to host the size of the dimensions
+ const int rank = m_dataspace.getSimpleExtentNdims();
+ // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
+ if (rank > m_max_dataset_dims)
+ error("Dataset dimensions (" + std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims)
+ + ". This is unexpected, Aborting.");
+ hsize_t dims[m_max_dataset_dims];
+ // Read size of dimensions
+ m_dataspace.getSimpleExtentDims(dims, NULL);
+
+ // AB for signa, these where [1981,27,357] and [45,448,357]
+ m_NX_SUB = dims[0]; // hyperslab dimensions
+ m_NY_SUB = dims[1];
+ m_NZ_SUB = dims[2];
#if 0
// AB todo: ??? why are these different?
m_NX = 45; // output buffer dimensions
@@ -647,255 +681,257 @@ Succeeded GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num)
m_NZ = 357;
#endif
}
- else
- return Succeeded::no;
+ else
+ return Succeeded::no;
- return Succeeded::yes;
+ return Succeeded::yes;
}
// PW The geo factors are stored in geo3d file under the file path called /SegmentData/Segment4/3D_Norm_correction/slice%d where
// slice numbers go from 1 to 16. Here this path is initialised, along with the output buffer and hyperslab.
//
-Succeeded GEHDF5Wrapper::initialise_geo_factors_data(const unsigned int slice_num)
+Succeeded
+GEHDF5Wrapper::initialise_geo_factors_data(const unsigned int slice_num)
{
- if(!is_geo_file())
- error("The file provided is not geometry data. Aborting");
+ if (!is_geo_file())
+ error("The file provided is not geometry data. Aborting");
- if(slice_num == 0)
- error("internal error in GE HDF5 geo code: slice number "+ std::to_string(slice_num) +" is incorrect");
+ if (slice_num == 0)
+ error("internal error in GE HDF5 geo code: slice number " + std::to_string(slice_num) + " is incorrect");
- if(rdf_ver==9)
+ if (rdf_ver == 9)
{
- m_address = "/SegmentData/Segment4/3D_Norm_Correction/slice";
- {
- // Open Dataset and get Dataspace(metadata)
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + std::to_string(slice_num))));
- m_dataspace = m_dataset_sptr->getSpace();
-
- // Read dimensions
- const int rank = m_dataspace.getSimpleExtentNdims();
- // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
- if (rank>m_max_dataset_dims)
- error("Dataset dimensions ("+ std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
- hsize_t dims[m_max_dataset_dims];
-
- m_dataspace.getSimpleExtentDims( dims, NULL);
-
- m_NX_SUB = dims[0]; // hyperslab dimensions
- m_NY_SUB = dims[1];
- m_NZ_SUB = (rank==2)? 1 : dims[2]; // Signa has rank==2, but this stay shere just in case...
+ m_address = "/SegmentData/Segment4/3D_Norm_Correction/slice";
+ {
+ // Open Dataset and get Dataspace(metadata)
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + std::to_string(slice_num))));
+ m_dataspace = m_dataset_sptr->getSpace();
+
+ // Read dimensions
+ const int rank = m_dataspace.getSimpleExtentNdims();
+ // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
+ if (rank > m_max_dataset_dims)
+ error("Dataset dimensions (" + std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims)
+ + ". This is unexpected, Aborting.");
+ hsize_t dims[m_max_dataset_dims];
+
+ m_dataspace.getSimpleExtentDims(dims, NULL);
+
+ m_NX_SUB = dims[0]; // hyperslab dimensions
+ m_NY_SUB = dims[1];
+ m_NZ_SUB = (rank == 2) ? 1 : dims[2]; // Signa has rank==2, but this stay shere just in case...
#if 0
m_NX = dims[0]; // output buffer dimensions
m_NY = dims[1];
m_NZ = (rank==2)? 1 : dims[2]; // Signa has rank==2, but this stay shere just in case...
#endif
- }
+ }
}
- else
- return Succeeded::no;
+ else
+ return Succeeded::no;
- return Succeeded::yes;
+ return Succeeded::yes;
}
// This info is in norm3d file
-Succeeded GEHDF5Wrapper::initialise_efficiency_factors()
+Succeeded
+GEHDF5Wrapper::initialise_efficiency_factors()
{
- if(!is_norm_file())
- error("The file provided is not norm data. Aborting");
+ if (!is_norm_file())
+ error("The file provided is not norm data. Aborting");
- if(rdf_ver==9)
+ if (rdf_ver == 9)
{
-
- m_address = "/3DCrystalEfficiency/crystalEfficiency";
- // Get the DataSpace (metadata) corresponding to the DataSet that we want to read
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
- m_dataspace = m_dataset_sptr->getSpace();
- // Create an array to host the size of the dimensions
- const int rank = m_dataspace.getSimpleExtentNdims();
- // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
- if (rank>m_max_dataset_dims)
- error("Dataset dimensions ("+ std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims) + ". This is unexpected, Aborting.");
- hsize_t dims[m_max_dataset_dims];
- // Read size of dimensions
- m_dataspace.getSimpleExtentDims( dims, NULL);
-
- m_NX_SUB = dims[0]; // hyperslab dimensions
- // TODO: Why is this divided by 2??
- m_NY_SUB = dims[1]/2; // should be equal to scanner_sptr->get_num_detectors_per_ring();
- m_NZ_SUB = (rank==2)? 1 : dims[2];
+ m_address = "/3DCrystalEfficiency/crystalEfficiency";
+ // Get the DataSpace (metadata) corresponding to the DataSet that we want to read
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address)));
+ m_dataspace = m_dataset_sptr->getSpace();
+ // Create an array to host the size of the dimensions
+ const int rank = m_dataspace.getSimpleExtentNdims();
+ // We allocate dims in the stack for efficiecy and safety but we need an error check just in case then
+ if (rank > m_max_dataset_dims)
+ error("Dataset dimensions (" + std::to_string(rank) + ") bigger than maximum of" + std::to_string(m_max_dataset_dims)
+ + ". This is unexpected, Aborting.");
+ hsize_t dims[m_max_dataset_dims];
+ // Read size of dimensions
+ m_dataspace.getSimpleExtentDims(dims, NULL);
+
+ m_NX_SUB = dims[0]; // hyperslab dimensions
+ // TODO: Why is this divided by 2??
+ m_NY_SUB = dims[1] / 2; // should be equal to scanner_sptr->get_num_detectors_per_ring();
+ m_NZ_SUB = (rank == 2) ? 1 : dims[2];
#if 0
m_NX = dims[0]; // output buffer dimensions
m_NY = dims[1]/2; // should be equal to scanner_sptr->get_num_detectors_per_ring();
m_NZ = (rank==2)? 1 : dims[2];
#endif
}
- else
- return Succeeded::no;
-
+ else
+ return Succeeded::no;
- return Succeeded::yes;
+ return Succeeded::yes;
}
// Developed for listmode access
-Succeeded GEHDF5Wrapper::read_list_data( char* output,
- const std::streampos offset, const hsize_t size) const
+Succeeded
+GEHDF5Wrapper::read_list_data(char* output, const std::streampos offset, const hsize_t size) const
{
- if(!is_list_file())
- error("The file provided is not list data. Aborting");
- const hsize_t pos = static_cast(offset);
- m_dataspace.selectHyperslab( H5S_SELECT_SET, &size, &pos );
- const H5::DataSpace memspace( m_dataset_list_Ndims, &size);
- m_dataset_sptr->read( output, H5::PredType::STD_U8LE, memspace, m_dataspace );
-
- return Succeeded::yes;
+ if (!is_list_file())
+ error("The file provided is not list data. Aborting");
+ const hsize_t pos = static_cast(offset);
+ m_dataspace.selectHyperslab(H5S_SELECT_SET, &size, &pos);
+ const H5::DataSpace memspace(m_dataset_list_Ndims, &size);
+ m_dataset_sptr->read(output, H5::PredType::STD_U8LE, memspace, m_dataspace);
+
+ return Succeeded::yes;
}
// Developed for ProjData
-Succeeded GEHDF5Wrapper::read_sinogram(Array<3, unsigned char> &output,
- const std::array& offset,
- const std::array& stride)
+Succeeded
+GEHDF5Wrapper::read_sinogram(Array<3, unsigned char>& output,
+ const std::array& offset,
+ const std::array& stride)
{
- // AB: this is only used for proj data, so for now lets ensure the file is sino. If its reused, we can change this.
- if(!is_sino_file())
- error("File is not sinogram. Aborting");
+ // AB: this is only used for proj data, so for now lets ensure the file is sino. If its reused, we can change this.
+ if (!is_sino_file())
+ error("File is not sinogram. Aborting");
- if(offset[0] != 0 || offset[1] != 0 || offset[2] != 0) //AB there are other C++ ways of doing this, but this is the most readable really.
- error("Only {0,0,0} offset supported. Aborting");
- if(stride[0] != 1 || stride[1] != 1 || stride[2] != 1)
- error("Only {1,1,1} stride supported. Aborting");
+ if (offset[0] != 0 || offset[1] != 0
+ || offset[2] != 0) // AB there are other C++ ways of doing this, but this is the most readable really.
+ error("Only {0,0,0} offset supported. Aborting");
+ if (stride[0] != 1 || stride[1] != 1 || stride[2] != 1)
+ error("Only {1,1,1} stride supported. Aborting");
- // We know the size of the DataSpace
- hsize_t str_dimsf[3] {m_NX_SUB, m_NY_SUB, m_NZ_SUB} ;
+ // We know the size of the DataSpace
+ hsize_t str_dimsf[3]{ m_NX_SUB, m_NY_SUB, m_NZ_SUB };
+ // get a temporary buffer here,
+ std::vector aux_buffer(m_NX_SUB * m_NY_SUB * m_NZ_SUB);
- // get a temporary buffer here,
- std::vector aux_buffer(m_NX_SUB*m_NY_SUB*m_NZ_SUB);
-
- m_dataspace.selectHyperslab(H5S_SELECT_SET, str_dimsf, offset.data());
- H5::DataSpace memspace(3, str_dimsf);
- m_dataset_sptr->read(static_cast(aux_buffer.data()), H5::PredType::STD_U8LE, memspace, m_dataspace);
+ m_dataspace.selectHyperslab(H5S_SELECT_SET, str_dimsf, offset.data());
+ H5::DataSpace memspace(3, str_dimsf);
+ m_dataset_sptr->read(static_cast(aux_buffer.data()), H5::PredType::STD_U8LE, memspace, m_dataspace);
- // the data is not in the correct size if its RDF9, so we will need to transpose the output of the data read.
- if(rdf_ver==9)
+ // the data is not in the correct size if its RDF9, so we will need to transpose the output of the data read.
+ if (rdf_ver == 9)
{
- output.resize(IndexRange3D(m_NZ_SUB,m_NY_SUB,m_NX_SUB));
- // now transpose the output.
- // AB: todo
- // todo 1: test
- for(unsigned int i=0;i& output,
+ const std::array& offset,
+ const std::array& count,
+ const std::array& stride)
-//PW Developed for Geometric Correction Factors
-Succeeded GEHDF5Wrapper::read_geometric_factors(Array<1, unsigned int> &output,
- const std::array& offset,
- const std::array& count,
- const std::array& stride)
-
{
- // AB: this is only used for geo data, so for now lets ensure the file is sino. If its reused, we can change this.
- if(!is_geo_file())
- error("File is Geometry. Aborting");
-
- if(count[0] == 0 || count[1] == 0) //AB there are other C++ ways of doing this, but this is the most readable really.
- error("Requested zero data to read. Aborting");
- if(stride[0] != 1 || stride[1] != 1)
- error("Only {1,1} stride supported. Aborting");
-
- Array<1, unsigned int> aux_reader;
- output.resize(count[0]*count[1]);
- aux_reader.resize(count[0]*count[1]);
-
- m_dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data());
- H5::DataSpace memspace(2, count.data());
- m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_UINT32, memspace, m_dataspace);
- aux_reader.release_data_ptr();
-
- // GE/RDF9 stores the tangetial axis reversed to STIR. Flip.
- for(unsigned int ax=0;ax aux_reader;
+ output.resize(count[0] * count[1]);
+ aux_reader.resize(count[0] * count[1]);
+
+ m_dataspace.selectHyperslab(H5S_SELECT_SET, count.data(), offset.data());
+ H5::DataSpace memspace(2, count.data());
+ m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_UINT32, memspace, m_dataspace);
+ aux_reader.release_data_ptr();
+
+ // GE/RDF9 stores the tangetial axis reversed to STIR. Flip.
+ for (unsigned int ax = 0; ax < count[0]; ++ax)
+ for (unsigned int tan = 0; tan < count[1]; ++tan)
+ output[ax * count[1] + ((count[1] - 1) - tan)] = aux_reader[ax * count[1] + tan];
+ return Succeeded::yes;
}
-//PW Developed for Efficiency Factors
-Succeeded GEHDF5Wrapper::read_efficiency_factors(Array<1, float> &output,
- const std::array& offset,
- const std::array& stride)
-
+// PW Developed for Efficiency Factors
+Succeeded
+GEHDF5Wrapper::read_efficiency_factors(Array<1, float>& output,
+ const std::array& offset,
+ const std::array& stride)
+
{
- if(!is_norm_file())
- error("The file provided is not norm data. Aborting");
-
- if(offset[0] != 0 || offset[1] != 0) //AB there are other C++ ways of doing this, but this is the most readable really.
- error("Only {0,0} offset supported. Aborting");
- if(stride[0] != 1 || stride[1] != 1)
- error("Only {1,1} stride supported. Aborting");
-
- // We know the size of the DataSpace
- hsize_t str_dimsf[2] {m_NX_SUB, m_NY_SUB} ;
- Array<1, float> aux_reader;
- output.resize(m_NX_SUB*m_NY_SUB);
- aux_reader.resize(m_NX_SUB*m_NY_SUB);
-
- m_dataspace.selectHyperslab(H5S_SELECT_SET, str_dimsf, offset.data());
- H5::DataSpace memspace(2, str_dimsf);
- m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_FLOAT, memspace, m_dataspace);
- aux_reader.release_data_ptr();
-
- // GE/RDF9 stores the tangetial axis reversed to STIR. Flip.
- for(unsigned int ax=0;ax aux_reader;
+ output.resize(m_NX_SUB * m_NY_SUB);
+ aux_reader.resize(m_NX_SUB * m_NY_SUB);
+
+ m_dataspace.selectHyperslab(H5S_SELECT_SET, str_dimsf, offset.data());
+ H5::DataSpace memspace(2, str_dimsf);
+ m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_FLOAT, memspace, m_dataspace);
+ aux_reader.release_data_ptr();
+
+ // GE/RDF9 stores the tangetial axis reversed to STIR. Flip.
+ for (unsigned int ax = 0; ax < m_NX_SUB; ++ax)
+ for (unsigned int tan = 0; tan < m_NY_SUB; ++tan)
+ output[ax * m_NY_SUB + ((m_NY_SUB - 1) - tan)] = aux_reader[ax * m_NY_SUB + tan];
+
+ return Succeeded::yes;
}
// Developed for Singles
-Succeeded GEHDF5Wrapper::read_singles(Array<1, unsigned int>& output, const unsigned int current_id)
+Succeeded
+GEHDF5Wrapper::read_singles(Array<1, unsigned int>& output, const unsigned int current_id)
{
- BOOST_STATIC_ASSERT(sizeof(unsigned int)==4); // Compilation time assert.
- if(!is_list_file()&& !is_sino_file())
- error("The file provided is not listmode or sinogram data. Aborting");
+ BOOST_STATIC_ASSERT(sizeof(unsigned int) == 4); // Compilation time assert.
+ if (!is_list_file() && !is_sino_file())
+ error("The file provided is not listmode or sinogram data. Aborting");
- if(current_id == 0 || current_id > this->m_num_singles_samples)
- error("internal error in GE HDF5 code: singles slice_id "+ std::to_string(current_id) +" is incorrect");
+ if (current_id == 0 || current_id > this->m_num_singles_samples)
+ error("internal error in GE HDF5 code: singles slice_id " + std::to_string(current_id) + " is incorrect");
- Array<1, unsigned int> aux_reader;
- output.resize(m_NX_SUB*m_NY_SUB);
- aux_reader.resize(m_NX_SUB*m_NY_SUB);
+ Array<1, unsigned int> aux_reader;
+ output.resize(m_NX_SUB * m_NY_SUB);
+ aux_reader.resize(m_NX_SUB * m_NY_SUB);
- // AB: todo check if output allocated data size is correct.
- m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + std::to_string(current_id))));
- m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_UINT32);
- aux_reader.release_data_ptr();
+ // AB: todo check if output allocated data size is correct.
+ m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address + std::to_string(current_id))));
+ m_dataset_sptr->read(aux_reader.get_data_ptr(), H5::PredType::NATIVE_UINT32);
+ aux_reader.release_data_ptr();
- // GE/RDF9 stores the tangetial axis reversed to STIR. Flip.
- for(unsigned int ax=0;axm_image_type=data_type_case;
- this->MaxLength=num_voxels;
- if(this->m_image_type == 64)
+ this->m_image_type = data_type_case;
+ this->MaxLength = num_voxels;
+ if (this->m_image_type == 64)
{
vData_f = new float[this->MaxLength];
for (int i = 0; i < this->MaxLength; i++)
- this->vData_f[i] = 0.F;
+ this->vData_f[i] = 0.F;
vData = NULL;
}
- else if(this->m_image_type == 15)
+ else if (this->m_image_type == 15)
{
vData = new short[this->MaxLength];
for (int i = 0; i < this->MaxLength; i++)
- this->vData[i] = 0;
+ this->vData[i] = 0;
vData_f = NULL;
}
else
{
- vData = NULL;
- vData_f = NULL;
+ vData = NULL;
+ vData_f = NULL;
}
- vDownsample[0] = 1;
- vDownsample[1] = 1;
- vDownsample[2] = 1;
-
- vCenter[0] = 0;
- vCenter[1] = 0;
- vCenter[2] = 0;
- vCenter[3] = 0;
+ vDownsample[0] = 1;
+ vDownsample[1] = 1;
+ vDownsample[2] = 1;
+
+ vCenter[0] = 0;
+ vCenter[1] = 0;
+ vCenter[2] = 0;
+ vCenter[3] = 0;
}
// -------------------------------------------------------------------------
@@ -109,10 +109,10 @@ Image::Image(const int num_voxels, const short data_type_case)
Image::~Image()
{
- if (vData)
- delete[] vData;
- if (vData_f)
- delete[] vData_f;
+ if (vData)
+ delete[] vData;
+ if (vData_f)
+ delete[] vData_f;
}
// -------------------------------------------------------------------------
@@ -121,555 +121,572 @@ Image::~Image()
// -------------------------------------------------------------------------
/**
-* \brief Read data from GIPL filename.
-*
-* \param filename Input filename
-*/
+ * \brief Read data from GIPL filename.
+ *
+ * \param filename Input filename
+ */
// -------------------------------------------------------------------------
-void Image::GiplRead(char* filename)
+void
+Image::GiplRead(char* filename)
{
- printf("Read %s\n",filename);
-
- // Open input binary file
- std::fstream myFile(filename, std::ios_base::in | std::ios_base::binary);
-
- // Initialize image dimensions by zero
- m_dim[0] = 0;
- m_dim[1] = 0;
- m_dim[2] = 0;
-
- // open file for reading
- /*char *buf;
- try {
- buf = new char[512];
- if( buf == 0 )
- throw "Memory allocation failure!";
- }
- catch( char * str ) {
- cout << "Exception raised: " << str << '\n';
- }*/
-
- // Read gipl header
- ReadGiplHeader(&myFile);
-
- // Require big endian format
- ByteSwapHeader();
-
- // Length of the data to be read in
- MaxLength = m_dim[0]*m_dim[1]*m_dim[2];
-
- // File cannot be read, wrong path or filename
- if (MaxLength == 0)
- {
- printf("Error: File %s cannot be found\n",filename);
- exit(1);
- }
-
- // Update image dimension
- ImageDimension = 2;
- if(m_dim[2]>1) ImageDimension = 3;
-
- // Update parameters dimension
- ParametersDimension = ImageDimension*ImageDimension+ImageDimension;
-
- // Set offset vector
- ImageOffset[0] = m_dim[0]; // XLen
- ImageOffset[1] = m_dim[0]*m_dim[1]; // XLen*YLen
-
- // Set center of rotation of the image
- for (int i = 0; i < ImageDimension; i++)
- vCenter[i] = 0;//m_dim[i]/2.0*m_pixdim[i] - m_origin[i];
-
- // Delete possible old data vectors
- if (vData)
- delete[] vData;
- if (vData_f)
- delete[] vData_f;
-
-
-
- // Check data type
- switch(m_image_type)
- {
- case 15: // shorts (default)
-
- // Initialize image vector
- vData = new short[MaxLength];
-
- // Read image data
- myFile.read((char*)vData, sizeof(short)*MaxLength);
-
- // Swap hi lo bytes
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData[i]);
-
- // Get min and max gray value
- GetMinMaxValue();
-
- break;
-
- case 64: // floats
-
- // Initialize image vector
- vData_f = new float[MaxLength];
-
- // Read image data
- myFile.read((char*)vData_f, sizeof(float)*MaxLength);
-
- // Swap hi lo bytes
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData_f[i]);
-
- break;
-
- default:
- break;
- }
-
- //int SourceIndex = ((int)(5 + 10*ImageOffset[0] + 20*ImageOffset[1]));
- //float test = vData_f[SourceIndex];
-
- // Initially no downscaling
- vDownsample[0] = 1;
- vDownsample[1] = 1;
- vDownsample[2] = 1;
-
- // Close file
- myFile.close();
+ printf("Read %s\n", filename);
+
+ // Open input binary file
+ std::fstream myFile(filename, std::ios_base::in | std::ios_base::binary);
+
+ // Initialize image dimensions by zero
+ m_dim[0] = 0;
+ m_dim[1] = 0;
+ m_dim[2] = 0;
+
+ // open file for reading
+ /*char *buf;
+try {
+buf = new char[512];
+if( buf == 0 )
+ throw "Memory allocation failure!";
+}
+catch( char * str ) {
+cout << "Exception raised: " << str << '\n';
+}*/
+
+ // Read gipl header
+ ReadGiplHeader(&myFile);
+
+ // Require big endian format
+ ByteSwapHeader();
+
+ // Length of the data to be read in
+ MaxLength = m_dim[0] * m_dim[1] * m_dim[2];
+
+ // File cannot be read, wrong path or filename
+ if (MaxLength == 0)
+ {
+ printf("Error: File %s cannot be found\n", filename);
+ exit(1);
+ }
+
+ // Update image dimension
+ ImageDimension = 2;
+ if (m_dim[2] > 1)
+ ImageDimension = 3;
+
+ // Update parameters dimension
+ ParametersDimension = ImageDimension * ImageDimension + ImageDimension;
+
+ // Set offset vector
+ ImageOffset[0] = m_dim[0]; // XLen
+ ImageOffset[1] = m_dim[0] * m_dim[1]; // XLen*YLen
+
+ // Set center of rotation of the image
+ for (int i = 0; i < ImageDimension; i++)
+ vCenter[i] = 0; // m_dim[i]/2.0*m_pixdim[i] - m_origin[i];
+
+ // Delete possible old data vectors
+ if (vData)
+ delete[] vData;
+ if (vData_f)
+ delete[] vData_f;
+
+ // Check data type
+ switch (m_image_type)
+ {
+ case 15: // shorts (default)
+
+ // Initialize image vector
+ vData = new short[MaxLength];
+
+ // Read image data
+ myFile.read((char*)vData, sizeof(short) * MaxLength);
+
+ // Swap hi lo bytes
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData[i]);
+
+ // Get min and max gray value
+ GetMinMaxValue();
+
+ break;
+
+ case 64: // floats
+
+ // Initialize image vector
+ vData_f = new float[MaxLength];
+
+ // Read image data
+ myFile.read((char*)vData_f, sizeof(float) * MaxLength);
+
+ // Swap hi lo bytes
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData_f[i]);
+
+ break;
+
+ default:
+ break;
+ }
+
+ // int SourceIndex = ((int)(5 + 10*ImageOffset[0] + 20*ImageOffset[1]));
+ // float test = vData_f[SourceIndex];
+
+ // Initially no downscaling
+ vDownsample[0] = 1;
+ vDownsample[1] = 1;
+ vDownsample[2] = 1;
+
+ // Close file
+ myFile.close();
}
// -------------------------------------------------------------------------
/**
-* \brief Get minimum and maximum gray values in image.
-*/
+ * \brief Get minimum and maximum gray values in image.
+ */
// -------------------------------------------------------------------------
-void Image::GetMinMaxValue()
+void
+Image::GetMinMaxValue()
{
- iMax = -16959;
- iMin = +16959;
-
- // Check data type
- if(m_image_type == 15)
- {
- for( int i = 0; i < MaxLength; i++)
- {
- if (vData[i] > iMax)
- iMax = vData[i];
- if (vData[i] < iMin)
- iMin = vData[i];
- }
- }
- if(m_image_type == 64)
- {
- for( int i = 0; i < MaxLength; i++)
- {
- if (vData_f[i] > iMax)
- iMax = vData_f[i];
- if (vData_f[i] < iMin)
- iMin = vData_f[i];
- }
- }
+ iMax = -16959;
+ iMin = +16959;
+
+ // Check data type
+ if (m_image_type == 15)
+ {
+ for (int i = 0; i < MaxLength; i++)
+ {
+ if (vData[i] > iMax)
+ iMax = vData[i];
+ if (vData[i] < iMin)
+ iMin = vData[i];
+ }
+ }
+ if (m_image_type == 64)
+ {
+ for (int i = 0; i < MaxLength; i++)
+ {
+ if (vData_f[i] > iMax)
+ iMax = vData_f[i];
+ if (vData_f[i] < iMin)
+ iMin = vData_f[i];
+ }
+ }
}
// -------------------------------------------------------------------------
/**
-* \brief Read GIPL header.
-*
-* \param myFile Input file
-*/
+ * \brief Read GIPL header.
+ *
+ * \param myFile Input file
+ */
// -------------------------------------------------------------------------
-void Image::ReadGiplHeader(std::fstream* myFile)
+void
+Image::ReadGiplHeader(std::fstream* myFile)
{
- int i;
- for(i=0;i<4;i++)
- myFile->read(reinterpret_cast < char * > (&m_dim[i]), sizeof(m_dim[i]));
- myFile->read(reinterpret_cast < char * > (&m_image_type), sizeof(m_image_type));
- for(i=0;i<4;i++)
- myFile->read(reinterpret_cast < char * > (&m_pixdim[i]), sizeof(m_pixdim[i]));
- for(i=0;i<80;i++)
- myFile->read(reinterpret_cast < char * > (&m_patDesc[i]), sizeof(m_patDesc[i]));
- for(i=0;i<12;i++)
- myFile->read(reinterpret_cast < char * > (&m_matrixElements[i]), sizeof(m_matrixElements[i]));
- myFile->read(reinterpret_cast < char * > (&m_identifier), sizeof(m_identifier));
- for(i=0;i<28;i++)
- myFile->read(reinterpret_cast < char * > (&m_spare[i]), sizeof(m_spare[i]));
- myFile->read(reinterpret_cast < char * > (&m_orientationFlag), sizeof(m_orientationFlag));
- myFile->read(reinterpret_cast < char * > (&m_flag1), sizeof(m_flag1));
- myFile->read(reinterpret_cast < char * > (&m_min), sizeof(m_min));
- myFile->read(reinterpret_cast < char * > (&m_max), sizeof(m_max));
- for(i=0;i<4;i++)
- myFile->read(reinterpret_cast < char * > (&m_origin[i]), sizeof(m_origin[i]));
- myFile->read(reinterpret_cast < char * > (&m_pixval_offset), sizeof(m_pixval_offset));
- myFile->read(reinterpret_cast < char * > (&m_pixval_cal), sizeof(m_pixval_cal));
- myFile->read(reinterpret_cast < char * > (&m_user_def1), sizeof(m_user_def1));
- myFile->read(reinterpret_cast < char * > (&m_user_def2), sizeof(m_user_def2));
- myFile->read(reinterpret_cast < char * > (&m_magic_number), sizeof(m_magic_number));
-
+ int i;
+ for (i = 0; i < 4; i++)
+ myFile->read(reinterpret_cast(&m_dim[i]), sizeof(m_dim[i]));
+ myFile->read(reinterpret_cast(&m_image_type), sizeof(m_image_type));
+ for (i = 0; i < 4; i++)
+ myFile->read(reinterpret_cast(&m_pixdim[i]), sizeof(m_pixdim[i]));
+ for (i = 0; i < 80; i++)
+ myFile->read(reinterpret_cast(&m_patDesc[i]), sizeof(m_patDesc[i]));
+ for (i = 0; i < 12; i++)
+ myFile->read(reinterpret_cast(&m_matrixElements[i]), sizeof(m_matrixElements[i]));
+ myFile->read(reinterpret_cast(&m_identifier), sizeof(m_identifier));
+ for (i = 0; i < 28; i++)
+ myFile->read(reinterpret_cast(&m_spare[i]), sizeof(m_spare[i]));
+ myFile->read(reinterpret_cast(&m_orientationFlag), sizeof(m_orientationFlag));
+ myFile->read(reinterpret_cast(&m_flag1), sizeof(m_flag1));
+ myFile->read(reinterpret_cast(&m_min), sizeof(m_min));
+ myFile->read(reinterpret_cast(&m_max), sizeof(m_max));
+ for (i = 0; i < 4; i++)
+ myFile->read(reinterpret_cast(&m_origin[i]), sizeof(m_origin[i]));
+ myFile->read(reinterpret_cast(&m_pixval_offset), sizeof(m_pixval_offset));
+ myFile->read(reinterpret_cast(&m_pixval_cal), sizeof(m_pixval_cal));
+ myFile->read(reinterpret_cast(&m_user_def1), sizeof(m_user_def1));
+ myFile->read(reinterpret_cast(&m_user_def2), sizeof(m_user_def2));
+ myFile->read(reinterpret_cast(&m_magic_number), sizeof(m_magic_number));
}
// -------------------------------------------------------------------------
/**
-* \brief Write image data to GIPL output file.
-*
-* \param filename Output filename
-*/
+ * \brief Write image data to GIPL output file.
+ *
+ * \param filename Output filename
+ */
// -------------------------------------------------------------------------
-void Image::GiplWrite(const char* filename)
+void
+Image::GiplWrite(const char* filename)
{
- // Open file for writing
- std::fstream myFile(filename, std::ios_base::out | std::ios_base::binary);
-
- // Require big endian format
- ByteSwapHeader();
-
- // Write gipl header
- WriteGiplHeader(&myFile);
-
- // Go back to correct format
- ByteSwapHeader();
-
- // Check data type
- switch(m_image_type)
- {
- case 15: // shorts (default)
-
- // Swap bytes of data vector before writing to file
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData[i]);
-
- // Read image data
- myFile.write((char*)vData, sizeof(short)*MaxLength);
-
- // Swap hi lo bytes
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData[i]);
- break;
-
- // Swap bytes back
- for(int i = 0; i < MaxLength; i++)
- ByteSwap(&vData[i]);
-
- case 64: // floats
-
- // Swap bytes of data vector before writing to file
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData_f[i]);
-
- // Read image data
- myFile.write((char*)vData_f, sizeof(float)*MaxLength);
-
- // Swap hi lo bytes
- for( int i = 0; i < MaxLength; i++)
- ByteSwap(&vData_f[i]);
- break;
-
- // Swap bytes back
- for(int i = 0; i < MaxLength; i++)
- ByteSwap(&vData_f[i]);
-
- default:
- break;
- }
-
- // Close file
- myFile.close();
+ // Open file for writing
+ std::fstream myFile(filename, std::ios_base::out | std::ios_base::binary);
+
+ // Require big endian format
+ ByteSwapHeader();
+
+ // Write gipl header
+ WriteGiplHeader(&myFile);
+
+ // Go back to correct format
+ ByteSwapHeader();
+
+ // Check data type
+ switch (m_image_type)
+ {
+ case 15: // shorts (default)
+
+ // Swap bytes of data vector before writing to file
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData[i]);
+
+ // Read image data
+ myFile.write((char*)vData, sizeof(short) * MaxLength);
+
+ // Swap hi lo bytes
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData[i]);
+ break;
+
+ // Swap bytes back
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData[i]);
+
+ case 64: // floats
+
+ // Swap bytes of data vector before writing to file
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData_f[i]);
+
+ // Read image data
+ myFile.write((char*)vData_f, sizeof(float) * MaxLength);
+
+ // Swap hi lo bytes
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData_f[i]);
+ break;
+
+ // Swap bytes back
+ for (int i = 0; i < MaxLength; i++)
+ ByteSwap(&vData_f[i]);
+
+ default:
+ break;
+ }
+
+ // Close file
+ myFile.close();
}
// -------------------------------------------------------------------------
/**
-* \brief Write GIPL header to output file.
-*
-* \param myfile Output file
-*/
+ * \brief Write GIPL header to output file.
+ *
+ * \param myfile Output file
+ */
// -------------------------------------------------------------------------
-void Image::WriteGiplHeader(fstream* myFile)
+void
+Image::WriteGiplHeader(fstream* myFile)
{
- int i;
- for(i=0;i<4;i++)
- myFile->write(reinterpret_cast < char * > (&m_dim[i]), sizeof(m_dim[i]));
- myFile->write(reinterpret_cast < char * > (&m_image_type), sizeof(m_image_type));
- for(i=0;i<4;i++)
- myFile->write(reinterpret_cast < char * > (&m_pixdim[i]), sizeof(m_pixdim[i]));
- for(i=0;i<80;i++)
- myFile->write(reinterpret_cast < char * > (&m_patDesc[i]), sizeof(m_patDesc[i]));
- for(i=0;i<12;i++)
- myFile->write(reinterpret_cast < char * > (&m_matrixElements[i]), sizeof(m_matrixElements[i]));
- myFile->write(reinterpret_cast < char * > (&m_identifier), sizeof(m_identifier));
- for(i=0;i<28;i++)
- myFile->write(reinterpret_cast < char * > (&m_spare[i]), sizeof(m_spare[i]));
- myFile->write(reinterpret_cast < char * > (&m_orientationFlag), sizeof(m_orientationFlag));
- myFile->write(reinterpret_cast < char * > (&m_flag1), sizeof(m_flag1));
- myFile->write(reinterpret_cast < char * > (&m_min), sizeof(m_min));
- myFile->write(reinterpret_cast < char * > (&m_max), sizeof(m_max));
- for(i=0;i<4;i++)
- myFile->write(reinterpret_cast < char * > (&m_origin[i]), sizeof(m_origin[i]));
- myFile->write(reinterpret_cast < char * > (&m_pixval_offset), sizeof(m_pixval_offset));
- myFile->write(reinterpret_cast < char * > (&m_pixval_cal), sizeof(m_pixval_cal));
- myFile->write(reinterpret_cast < char * > (&m_user_def1), sizeof(m_user_def1));
- myFile->write(reinterpret_cast < char * > (&m_user_def2), sizeof(m_user_def2));
- myFile->write(reinterpret_cast < char * > (&m_magic_number), sizeof(m_magic_number));
+ int i;
+ for (i = 0; i < 4; i++)
+ myFile->write(reinterpret_cast(&m_dim[i]), sizeof(m_dim[i]));
+ myFile->write(reinterpret_cast(&m_image_type), sizeof(m_image_type));
+ for (i = 0; i < 4; i++)
+ myFile->write(reinterpret_cast(&m_pixdim[i]), sizeof(m_pixdim[i]));
+ for (i = 0; i < 80; i++)
+ myFile->write(reinterpret_cast(&m_patDesc[i]), sizeof(m_patDesc[i]));
+ for (i = 0; i < 12; i++)
+ myFile->write(reinterpret_cast(&m_matrixElements[i]), sizeof(m_matrixElements[i]));
+ myFile->write(reinterpret_cast(&m_identifier), sizeof(m_identifier));
+ for (i = 0; i < 28; i++)
+ myFile->write(reinterpret_cast(&m_spare[i]), sizeof(m_spare[i]));
+ myFile->write(reinterpret_cast(&m_orientationFlag), sizeof(m_orientationFlag));
+ myFile->write(reinterpret_cast(&m_flag1), sizeof(m_flag1));
+ myFile->write(reinterpret_cast(&m_min), sizeof(m_min));
+ myFile->write(reinterpret_cast(&m_max), sizeof(m_max));
+ for (i = 0; i < 4; i++)
+ myFile->write(reinterpret_cast(&m_origin[i]), sizeof(m_origin[i]));
+ myFile->write(reinterpret_cast(&m_pixval_offset), sizeof(m_pixval_offset));
+ myFile->write(reinterpret_cast(&m_pixval_cal), sizeof(m_pixval_cal));
+ myFile->write(reinterpret_cast(&m_user_def1), sizeof(m_user_def1));
+ myFile->write(reinterpret_cast(&m_user_def2), sizeof(m_user_def2));
+ myFile->write(reinterpret_cast(&m_magic_number), sizeof(m_magic_number));
}
// -------------------------------------------------------------------------
/**
-* \brief Swap bytes (little/big endian conversion).
-*/
+ * \brief Swap bytes (little/big endian conversion).
+ */
// -------------------------------------------------------------------------
-void Image::ByteSwapHeader() // for PC little endian
-{
- int i;
- for(i=0;i<4;i++) ByteSwap(&(m_dim[i]));
+void
+Image::ByteSwapHeader() // for PC little endian
+{
+ int i;
+ for (i = 0; i < 4; i++)
+ ByteSwap(&(m_dim[i]));
+
+ ByteSwap(&(m_image_type));
- ByteSwap(&(m_image_type));
+ for (i = 0; i < 4; i++)
+ ByteSwap(&(m_pixdim[i]));
- for(i=0;i<4;i++) ByteSwap(&(m_pixdim[i]));
+ for (i = 0; i < 12; i++)
+ ByteSwap(&(m_matrixElements[i]));
- for(i=0;i<12;i++) ByteSwap(&(m_matrixElements[i]));
+ ByteSwap(&(m_min));
+ ByteSwap(&(m_max));
- ByteSwap(&(m_min));
- ByteSwap(&(m_max));
+ for (i = 0; i < 4; i++)
+ ByteSwap(&(m_origin[i]));
- for(i=0;i<4;i++) ByteSwap(&(m_origin[i]));
-
- ByteSwap(&(m_pixval_offset));
- ByteSwap(&(m_pixval_cal));
- ByteSwap(&(m_user_def1));
- ByteSwap(&(m_user_def2));
- ByteSwap((int*)&(m_magic_number));
+ ByteSwap(&(m_pixval_offset));
+ ByteSwap(&(m_pixval_cal));
+ ByteSwap(&(m_user_def1));
+ ByteSwap(&(m_user_def2));
+ ByteSwap((int*)&(m_magic_number));
- return;
+ return;
}
// -------------------------------------------------------------------------
/**
-* \brief Initialize zero image
-*
-* \param Input Input image
-*/
+ * \brief Initialize zero image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Zeros(Image* Input, short iType)
+void
+Image::Zeros(Image* Input, short iType)
{
- // Initialize image with dimensions
- this->Initialize(Input, iType);
-
- if (iType == _SHORT)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData[i] = 0;
-
- // Update min max values
- //this->GetMinMaxValue();
- }
-
- if (iType == _FLOAT)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData_f[i] = 0;
- }
-
- iMin = 0;
- iMax = 0;
+ // Initialize image with dimensions
+ this->Initialize(Input, iType);
+
+ if (iType == _SHORT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData[i] = 0;
+
+ // Update min max values
+ // this->GetMinMaxValue();
+ }
+
+ if (iType == _FLOAT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData_f[i] = 0;
+ }
+
+ iMin = 0;
+ iMax = 0;
}
// -------------------------------------------------------------------------
/**
-* \brief Initialize zero image
-*
-* \param Input Input image
-*/
+ * \brief Initialize zero image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Zeros()
+void
+Image::Zeros()
{
- // Set image to 0
- if (this->m_image_type == 15)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData[i] = 0;
-
- // Update min max values
- //this->GetMinMaxValue();
- }
-
- if (this->m_image_type == 64)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData_f[i] = 0;
- }
-
- iMin = 0;
- iMax = 0;
+ // Set image to 0
+ if (this->m_image_type == 15)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData[i] = 0;
+
+ // Update min max values
+ // this->GetMinMaxValue();
+ }
+
+ if (this->m_image_type == 64)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData_f[i] = 0;
+ }
+
+ iMin = 0;
+ iMax = 0;
}
// -------------------------------------------------------------------------
/**
-* \brief Initialize one image
-*
-* \param Input Input image
-*/
+ * \brief Initialize one image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Ones(Image* Input, short iType)
+void
+Image::Ones(Image* Input, short iType)
{
- // Initialize image with dimensions
- this->Initialize(Input, iType);
-
- if (iType == _SHORT)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData[i] = 1;
- }
-
- if (iType == _FLOAT)
- {
- // Initialize with zeros
- for (int i = 0; i < this->MaxLength; i++)
- this->vData_f[i] = 1;
- }
-
- iMax = 1;
- iMin = 1;
+ // Initialize image with dimensions
+ this->Initialize(Input, iType);
+
+ if (iType == _SHORT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData[i] = 1;
+ }
+
+ if (iType == _FLOAT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < this->MaxLength; i++)
+ this->vData_f[i] = 1;
+ }
+
+ iMax = 1;
+ iMin = 1;
}
// -------------------------------------------------------------------------
/**
-* \brief Initialize image data by with the dimension of the input image
-*
-* \param Input Input image
-*/
+ * \brief Initialize image data by with the dimension of the input image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Initialize(Image* Input, short iType)
+void
+Image::Initialize(Image* Input, short iType)
{
- // Set new data properties
- for (int i = 0; i < 4; i++)
- {
- this->m_dim[i] = Input->m_dim[i];
- this->m_pixdim[i] = Input->m_pixdim[i];
- this->vCenter[i] = Input->vCenter[i];
- this->m_origin[i] = Input->m_origin[i];
- }
-
- this->m_orientationFlag = Input->m_orientationFlag;
- this->m_flag1 = Input->m_flag1;
- this->m_min = Input->m_min;
- this->m_max = Input->m_max;
- this->m_pixval_offset = Input->m_pixval_offset;
- this->m_pixval_cal = Input->m_pixval_cal;
- this->m_user_def1 = Input->m_user_def1;
- this->m_user_def2 = Input->m_user_def2;
- this->m_magic_number = Input->m_magic_number;
- for (int i = 0; i < 2; i++)
- this->ImageOffset[i] = Input->ImageOffset[i]; // XLen and XLen*YLen
- this->MaxLength = Input->MaxLength;
-
- // Initialize data vactor
- if (iType == _SHORT)
- {
- // Delete old vector and initialize data vector
- if (this->vData)
- delete[] this->vData;
-
- this->m_image_type = 15;
- this->vData = new short[Input->MaxLength];
- }
-
- if (iType == _FLOAT)
- {
- // Delete old vector and initialize data vector
- if (this->vData_f)
- delete[] this->vData_f;
-
- this->m_image_type = 64;
- this->vData_f = new float[Input->MaxLength];
- }
+ // Set new data properties
+ for (int i = 0; i < 4; i++)
+ {
+ this->m_dim[i] = Input->m_dim[i];
+ this->m_pixdim[i] = Input->m_pixdim[i];
+ this->vCenter[i] = Input->vCenter[i];
+ this->m_origin[i] = Input->m_origin[i];
+ }
+
+ this->m_orientationFlag = Input->m_orientationFlag;
+ this->m_flag1 = Input->m_flag1;
+ this->m_min = Input->m_min;
+ this->m_max = Input->m_max;
+ this->m_pixval_offset = Input->m_pixval_offset;
+ this->m_pixval_cal = Input->m_pixval_cal;
+ this->m_user_def1 = Input->m_user_def1;
+ this->m_user_def2 = Input->m_user_def2;
+ this->m_magic_number = Input->m_magic_number;
+ for (int i = 0; i < 2; i++)
+ this->ImageOffset[i] = Input->ImageOffset[i]; // XLen and XLen*YLen
+ this->MaxLength = Input->MaxLength;
+
+ // Initialize data vactor
+ if (iType == _SHORT)
+ {
+ // Delete old vector and initialize data vector
+ if (this->vData)
+ delete[] this->vData;
+
+ this->m_image_type = 15;
+ this->vData = new short[Input->MaxLength];
+ }
+
+ if (iType == _FLOAT)
+ {
+ // Delete old vector and initialize data vector
+ if (this->vData_f)
+ delete[] this->vData_f;
+
+ this->m_image_type = 64;
+ this->vData_f = new float[Input->MaxLength];
+ }
}
// -------------------------------------------------------------------------
/**
-* \brief Initialize image data by with the dimension of the input image
-*
-* \param Input Input image
-*/
+ * \brief Initialize image data by with the dimension of the input image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Initialize(Image* Input)
+void
+Image::Initialize(Image* Input)
{
- // Set new data properties
- for (int i = 0; i < 4; i++)
- {
- this->m_dim[i] = Input->m_dim[i];
- this->m_pixdim[i] = Input->m_pixdim[i];
- this->vCenter[i] = Input->vCenter[i];
- this->m_origin[i] = Input->m_origin[i];
- }
-
- this->m_orientationFlag = Input->m_orientationFlag;
- this->m_flag1 = Input->m_flag1;
- this->m_min = Input->m_min;
- this->m_max = Input->m_max;
- this->m_pixval_offset = Input->m_pixval_offset;
- this->m_pixval_cal = Input->m_pixval_cal;
- this->m_user_def1 = Input->m_user_def1;
- this->m_user_def2 = Input->m_user_def2;
- this->m_magic_number = Input->m_magic_number;
- for (int i = 0; i < 2; i++)
- this->ImageOffset[i] = Input->ImageOffset[i]; // XLen and XLen*YLen
- this->MaxLength = Input->MaxLength;
-
- // Initialize data vactor
- if (this->vData)
- delete[] this->vData;
-
- if (this->vData_f)
- delete[] this->vData_f;
+ // Set new data properties
+ for (int i = 0; i < 4; i++)
+ {
+ this->m_dim[i] = Input->m_dim[i];
+ this->m_pixdim[i] = Input->m_pixdim[i];
+ this->vCenter[i] = Input->vCenter[i];
+ this->m_origin[i] = Input->m_origin[i];
+ }
+
+ this->m_orientationFlag = Input->m_orientationFlag;
+ this->m_flag1 = Input->m_flag1;
+ this->m_min = Input->m_min;
+ this->m_max = Input->m_max;
+ this->m_pixval_offset = Input->m_pixval_offset;
+ this->m_pixval_cal = Input->m_pixval_cal;
+ this->m_user_def1 = Input->m_user_def1;
+ this->m_user_def2 = Input->m_user_def2;
+ this->m_magic_number = Input->m_magic_number;
+ for (int i = 0; i < 2; i++)
+ this->ImageOffset[i] = Input->ImageOffset[i]; // XLen and XLen*YLen
+ this->MaxLength = Input->MaxLength;
+
+ // Initialize data vactor
+ if (this->vData)
+ delete[] this->vData;
+
+ if (this->vData_f)
+ delete[] this->vData_f;
}
// -------------------------------------------------------------------------
/**
-* \brief Copy image
-*
-* \param Input Input image
-*/
+ * \brief Copy image
+ *
+ * \param Input Input image
+ */
// -------------------------------------------------------------------------
-void Image::Copy(Image* Input, short iType)
+void
+Image::Copy(Image* Input, short iType)
{
- // Initialize image with dimensions
- this->Initialize(Input, iType);
-
- if (iType == _SHORT)
- {
- // Initialize with zeros
- for (int i = 0; i < Input->MaxLength; i++)
- this->vData[i] = Input->vData[i];
-
- // Update min max values
- this->GetMinMaxValue();
- }
-
- if (iType == _FLOAT)
- {
- // Initialize with zeros
- for (int i = 0; i < Input->MaxLength; i++)
- this->vData_f[i] = Input->vData_f[i];
- }
+ // Initialize image with dimensions
+ this->Initialize(Input, iType);
+
+ if (iType == _SHORT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < Input->MaxLength; i++)
+ this->vData[i] = Input->vData[i];
+
+ // Update min max values
+ this->GetMinMaxValue();
+ }
+
+ if (iType == _FLOAT)
+ {
+ // Initialize with zeros
+ for (int i = 0; i < Input->MaxLength; i++)
+ this->vData_f[i] = Input->vData_f[i];
+ }
}
// -------------------------------------------------------------------------
// Byte swap functions
// -------------------------------------------------------------------------
-void Image::ByteSwap(int *i) // for PC little endian
-{
- typedef struct {
+void
+Image::ByteSwap(int* i) // for PC little endian
+{
+ typedef struct
+ {
unsigned char byte1;
unsigned char byte2;
unsigned char byte3;
unsigned char byte4;
} Bytes;
- typedef union {
+ typedef union
+ {
int integer;
Bytes bytes;
} intUnion;
@@ -681,19 +698,22 @@ void Image::ByteSwap(int *i) // for PC little endian
intUS.bytes.byte2 = intU.bytes.byte3;
intUS.bytes.byte3 = intU.bytes.byte2;
intUS.bytes.byte4 = intU.bytes.byte1;
-
+
*i = intUS.integer;
return;
}
-void Image::ByteSwap(short *s) // for PC little endian
-{
- typedef struct {
+void
+Image::ByteSwap(short* s) // for PC little endian
+{
+ typedef struct
+ {
unsigned char byte1;
unsigned char byte2;
} Bytes;
- typedef union {
+ typedef union
+ {
short shortint;
Bytes bytes;
} shortUnion;
@@ -708,16 +728,19 @@ void Image::ByteSwap(short *s) // for PC little endian
return;
}
-void Image::ByteSwap(float *f) // for PC little endian
-{
- typedef struct {
+void
+Image::ByteSwap(float* f) // for PC little endian
+{
+ typedef struct
+ {
unsigned char byte1;
unsigned char byte2;
unsigned char byte3;
unsigned char byte4;
} Bytes;
- typedef union {
+ typedef union
+ {
float floatnum;
Bytes bytes;
} floatUnion;
@@ -734,9 +757,11 @@ void Image::ByteSwap(float *f) // for PC little endian
return;
}
-void Image::ByteSwap(double *d) // for PC little endian
-{
- typedef struct {
+void
+Image::ByteSwap(double* d) // for PC little endian
+{
+ typedef struct
+ {
unsigned char byte1;
unsigned char byte2;
unsigned char byte3;
@@ -747,7 +772,8 @@ void Image::ByteSwap(double *d) // for PC little endian
unsigned char byte8;
} Bytes;
- typedef union {
+ typedef union
+ {
double doublenum;
Bytes bytes;
} doubleUnion;
diff --git a/src/IO/IO_registries.cxx b/src/IO/IO_registries.cxx
index 63f006168..e969d3ff0 100644
--- a/src/IO/IO_registries.cxx
+++ b/src/IO/IO_registries.cxx
@@ -3,7 +3,7 @@
Copyright (C) 2012, Kris Thielemans
Copyright (C) 2013, Institute for Bioengineering of Catalonia
Copyright (C) 2013, University College London
-
+
This file is part of STIR.
SPDX-License-Identifier: Apache-2.0
@@ -33,48 +33,47 @@
#include "stir/IO/MultiParametricDiscretisedDensityInputFileFormat.h"
#include "stir/IO/MultiParametricDiscretisedDensityOutputFileFormat.h"
#ifdef HAVE_LLN_MATRIX
-#include "stir/IO/ECAT6OutputFileFormat.h"
-#include "stir/IO/ECAT7OutputFileFormat.h"
-#include "stir/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.h"
-#include "stir/IO/ECAT7ParametricDensityOutputFileFormat.h"
-#include "stir/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.h"
+# include "stir/IO/ECAT6OutputFileFormat.h"
+# include "stir/IO/ECAT7OutputFileFormat.h"
+# include "stir/IO/ECAT7DynamicDiscretisedDensityOutputFileFormat.h"
+# include "stir/IO/ECAT7ParametricDensityOutputFileFormat.h"
+# include "stir/IO/ECAT7DynamicDiscretisedDensityInputFileFormat.h"
#endif
-
#if 1
-#include "stir/IO/InputFileFormatRegistry.h"
-#include "stir/IO/InterfileImageInputFileFormat.h"
-#ifdef HAVE_LLN_MATRIX
-#include "stir/IO/ECAT6ImageInputFileFormat.h"
-#include "stir/IO/ECAT7ImageInputFileFormat.h"
-#include "stir/IO/ECAT966ListmodeInputFileFormat.h"
-#include "stir/IO/ECAT962ListmodeInputFileFormat.h"
-#endif
+# include "stir/IO/InputFileFormatRegistry.h"
+# include "stir/IO/InterfileImageInputFileFormat.h"
+# ifdef HAVE_LLN_MATRIX
+# include "stir/IO/ECAT6ImageInputFileFormat.h"
+# include "stir/IO/ECAT7ImageInputFileFormat.h"
+# include "stir/IO/ECAT966ListmodeInputFileFormat.h"
+# include "stir/IO/ECAT962ListmodeInputFileFormat.h"
+# endif
#endif
#include "stir/IO/ECAT8_32bitListmodeInputFileFormat.h"
#ifdef HAVE_HDF5
-#include "stir/IO/GEHDF5ListmodeInputFileFormat.h"
+# include "stir/IO/GEHDF5ListmodeInputFileFormat.h"
#endif
//! Addition for SAFIR listmode input file format
#include "stir/IO/SAFIRCListmodeInputFileFormat.h"
//! Addition for ROOT support - Nikos Efthimiou
#ifdef HAVE_CERN_ROOT
-#include "stir/IO/ROOTListmodeInputFileFormat.h"
-#include "stir/IO/InputStreamFromROOTFileForCylindricalPET.h"
-#include "stir/IO/InputStreamFromROOTFileForECATPET.h"
+# include "stir/IO/ROOTListmodeInputFileFormat.h"
+# include "stir/IO/InputStreamFromROOTFileForCylindricalPET.h"
+# include "stir/IO/InputStreamFromROOTFileForECATPET.h"
#endif
#ifdef HAVE_UPENN
-#include "stir/IO/PENNListmodeInputFileFormat.h"
-#include "stir/IO/InputStreamWithRecordsFromUPENNbin.h"
-#include "stir/IO/InputStreamWithRecordsFromUPENNtxt.h"
+# include "stir/IO/PENNListmodeInputFileFormat.h"
+# include "stir/IO/InputStreamWithRecordsFromUPENNbin.h"
+# include "stir/IO/InputStreamWithRecordsFromUPENNtxt.h"
#endif
#ifdef HAVE_ITK
-#include "stir/IO/ITKOutputFileFormat.h"
-#include "stir/IO/ITKImageInputFileFormat.h"
+# include "stir/IO/ITKOutputFileFormat.h"
+# include "stir/IO/ITKImageInputFileFormat.h"
#endif
START_NAMESPACE_STIR
@@ -102,8 +101,6 @@ static InputStreamFromROOTFileForCylindricalPET::RegisterIt dummy60606;
static InputStreamFromROOTFileForECATPET::RegisterIt dummy606062;
#endif
-
-
#ifdef HAVE_LLN_MATRIX
START_NAMESPACE_ECAT
START_NAMESPACE_ECAT6
@@ -117,12 +114,11 @@ END_NAMESPACE_ECAT7
END_NAMESPACE_ECAT
#endif
-
static RegisterInputFileFormat idummy0(0);
#ifdef HAVE_LLN_MATRIX
static RegisterInputFileFormat idummy2(4);
-// ECAT6 very low priority it doesn't have a signature
+// ECAT6 very low priority it doesn't have a signature
static RegisterInputFileFormat idummy4(100000);
static RegisterInputFileFormat dynidummy(0);
@@ -130,15 +126,14 @@ static RegisterInputFileFormat > > idummy6(10000);
-static RegisterInputFileFormat > > > idummy7(10000);
+static RegisterInputFileFormat>> idummy6(10000);
+static RegisterInputFileFormat>>> idummy7(10000);
#endif
static RegisterInputFileFormat dyndummy_intf(1);
static RegisterInputFileFormat paradummy_intf(1);
static RegisterInputFileFormat dynim_dummy_multi(1);
static RegisterInputFileFormat parim_dummy_multi(1);
-
/*************************** listmode data **********************/
#ifdef HAVE_LLN_MATRIX
static RegisterInputFileFormat LMdummyECAT966(4);
@@ -153,8 +148,8 @@ static RegisterInputFileFormat LMdu
static RegisterInputFileFormat LMdummyPENN(8);
static InputStreamWithRecordsFromUPENNbin::RegisterIt dummy68606;
static InputStreamWithRecordsFromUPENNtxt::RegisterIt dummy686062;
-//static RegisterInputFileFormat LMdummyPENNbin(9);
-//static RegisterInputFileFormat idummy1(2);
+// static RegisterInputFileFormat LMdummyPENNbin(9);
+// static RegisterInputFileFormat idummy1(2);
#endif
END_NAMESPACE_STIR
diff --git a/src/IO/ITKImageInputFileFormat.cxx b/src/IO/ITKImageInputFileFormat.cxx
index 8835b22b7..c429e3609 100644
--- a/src/IO/ITKImageInputFileFormat.cxx
+++ b/src/IO/ITKImageInputFileFormat.cxx
@@ -48,74 +48,64 @@ START_NAMESPACE_STIR
*/
-typedef itk::Image ITKImageSingle;
-typedef itk::VectorImage ITKImageMulti;
-typedef DiscretisedDensity<3, float> STIRImageSingle;
-typedef VoxelsOnCartesianGrid STIRImageSingleConcrete;
-typedef DiscretisedDensity<3, CartesianCoordinate3D > STIRImageMulti;
-typedef VoxelsOnCartesianGrid > STIRImageMultiConcrete;
-typedef itk::MetaDataObject< std::string > MetaDataStringType;
+typedef itk::Image ITKImageSingle;
+typedef itk::VectorImage ITKImageMulti;
+typedef DiscretisedDensity<3, float> STIRImageSingle;
+typedef VoxelsOnCartesianGrid STIRImageSingleConcrete;
+typedef DiscretisedDensity<3, CartesianCoordinate3D> STIRImageMulti;
+typedef VoxelsOnCartesianGrid> STIRImageMultiConcrete;
+typedef itk::MetaDataObject MetaDataStringType;
// internal function to do the conversion. Note that it can throw an exception.
-template
-static
-STIRImageType *
-read_file_itk(const std::string &filename);
+template
+static STIRImageType* read_file_itk(const std::string& filename);
template
-bool
-ITKImageInputFileFormat::actual_can_read(const FileSignature& signature,
- std::istream& input) const
+bool
+ITKImageInputFileFormat::actual_can_read(const FileSignature& signature, std::istream& input) const
{
return false;
}
template
bool
-ITKImageInputFileFormat::can_read(const FileSignature& signature,
- std::istream& input) const
+ITKImageInputFileFormat::can_read(const FileSignature& signature, std::istream& input) const
{
return this->actual_can_read(signature, input);
}
template
-bool
-ITKImageInputFileFormat::can_read(const FileSignature& /*signature*/,
- const std::string& filename) const
+bool
+ITKImageInputFileFormat::can_read(const FileSignature& /*signature*/, const std::string& filename) const
{
typedef itk::ImageFileReader ReaderType;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(filename);
- try
- {
- reader->Update();
+ try
+ {
+ reader->Update();
return true;
- }
- catch( itk::ExceptionObject & /*err*/ )
- {
-
+ }
+ catch (itk::ExceptionObject& /*err*/)
+ {
+
return false;
- }
+ }
}
template
unique_ptr
ITKImageInputFileFormat::read_from_file(std::istream& input) const
{
- error("read_from_file for ITK with istream not implemented %s:%d. Sorry",
- __FILE__, __LINE__);
- return
- unique_ptr