create_empty_interfile_binary improvements using InterfileImageHeader…

… and InterfilePDFSHeader

Split functionality of create_image_and_header_from(). This was to grant access to create_image_from_header, where the argument is a InterfileImageHeader object with parsed data

src/IO/interfile.cxx

@@ -20,6 +20,7 @@
   \author PARAPET project
   \author Richard Brown
   \author Parisa Khateri
+  \author Robert Twyman
 */
 //   Pretty horrible implementations at the moment...
 
@@ -87,53 +88,58 @@ is_interfile_signature(const char * const signature)
 
 
 // help function
-static
-VoxelsOnCartesianGrid<float> *
-create_image_and_header_from(InterfileImageHeader& hdr,
-                             char * full_data_file_name, // preallocated
-                             istream& input, 
-                             const string&  directory_for_data)
+
+
+VoxelsOnCartesianGrid<float>
+create_image_from_header(stir::InterfileImageHeader& hdr)
 {
-  if (!hdr.parse(input))
-  {
-      return 0; //KT 10/12/2001 do not call ask_parameters anymore 
-    }
-
-  // prepend directory_for_data to the data_file_name from the header
-  strcpy(full_data_file_name, hdr.data_file_name.c_str());
-  prepend_directory_name(full_data_file_name, directory_for_data.c_str());  
-
-
-  CartesianCoordinate3D<float> voxel_size(static_cast<float>(hdr.pixel_sizes[2]), 
-					  static_cast<float>(hdr.pixel_sizes[1]), 
-					  static_cast<float>(hdr.pixel_sizes[0]));
-
+  CartesianCoordinate3D<float> voxel_size(static_cast<float>(hdr.pixel_sizes[2]),
+                                          static_cast<float>(hdr.pixel_sizes[1]),
+                                          static_cast<float>(hdr.pixel_sizes[0]));
+
   const int z_size =  hdr.matrix_size[2][0];
   const int y_size =  hdr.matrix_size[1][0];
   const int x_size =  hdr.matrix_size[0][0];
-  const BasicCoordinate<3,int> min_indices = 
-    make_coordinate(0, -y_size/2, -x_size/2);
-  const BasicCoordinate<3,int> max_indices = 
-    min_indices + make_coordinate(z_size, y_size, x_size) - 1;
+  const BasicCoordinate<3,int> min_indices = make_coordinate(0, -y_size/2, -x_size/2);
+  const BasicCoordinate<3,int> max_indices = min_indices + make_coordinate(z_size, y_size, x_size) - 1;
 
   CartesianCoordinate3D<float> origin(0,0,0);
   if (hdr.first_pixel_offsets[2] != InterfileHeader::double_value_not_set)
+  {
+    // make sure that origin is such that
+    // first_pixel_offsets =  min_indices*voxel_size + origin
+    origin =
+        make_coordinate(float(hdr.first_pixel_offsets[2]),
+                        float(hdr.first_pixel_offsets[1]),
+                        float(hdr.first_pixel_offsets[0]))
+        - voxel_size * BasicCoordinate<3,float>(min_indices);
+  }
+
+  return VoxelsOnCartesianGrid<float>
+          (hdr.get_exam_info_sptr(),
+           IndexRange<3>(min_indices, max_indices),
+           origin,
+           voxel_size);
+}
+
+// help function
+static
+VoxelsOnCartesianGrid<float> *
+create_image_and_header_from(InterfileImageHeader& hdr,
+                             char * full_data_file_name, // preallocated
+                             istream& input,
+                             const string&  directory_for_data)
+{
+  if (!hdr.parse(input))
     {
-      // make sure that origin is such that 
-      // first_pixel_offsets =  min_indices*voxel_size + origin
-      origin =
-	make_coordinate(float(hdr.first_pixel_offsets[2]),
-			float(hdr.first_pixel_offsets[1]),
-			float(hdr.first_pixel_offsets[0]))
-	- voxel_size * BasicCoordinate<3,float>(min_indices);
+      return 0; // KT 10/12/2001 do not call ask_parameters anymore
     }
 
-  return
-    new VoxelsOnCartesianGrid<float>
-    (hdr.get_exam_info_sptr(),
-     IndexRange<3>(min_indices, max_indices),
-     origin,
-     voxel_size);
+  // prepend directory_for_data to the data_file_name from the header
+  strcpy(full_data_file_name, hdr.data_file_name.c_str());
+  prepend_directory_name(full_data_file_name, directory_for_data.c_str());
+//  return create_image_from_header(hdr);
+  return new VoxelsOnCartesianGrid<float>(create_image_from_header(hdr));
 }
 
 VoxelsOnCartesianGrid<float> *

src/include/stir/IO/interfile.h

@@ -6,7 +6,7 @@
 /*
     Copyright (C) 2000 PARAPET partners
     Copyright (C) 2000- 2011, Hammersmith Imanet Ltd
-    Copyright (C) 2018, University College London
+    Copyright (C) 2018, 2024 University College London
     This file is part of STIR.
     SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
 
@@ -22,6 +22,7 @@
   \author Sanida Mustafovic
   \author PARAPET project
   \author Richard Brown
+  \author Robert Twyman
 
 */
 
@@ -31,6 +32,7 @@
 // has to include Succeeded.h (even if it doesn't use the return value).
 #include "stir/Succeeded.h"
 #include "stir/ByteOrder.h"
+#include "stir/IO/InterfileHeader.h"
 #include <iostream>
 #include <string>
 
@@ -82,6 +84,8 @@ bool is_interfile_signature(const char * const signature);
 VoxelsOnCartesianGrid<float>* read_interfile_image(std::istream& input, 
 						   const std::string& directory_for_data = "");
 
+VoxelsOnCartesianGrid<float> create_image_from_header(InterfileImageHeader& hdr);
+
 //! This reads the first 3d image in an Interfile header file, given as a filename.
 /*!
   \ingroup InterfileIO

src/utilities/CMakeLists.txt

@@ -74,7 +74,7 @@ set(${dir_EXE_SOURCES}
         find_sum_projection_of_viewgram_and_sinogram.cxx
         separate_true_from_random_scatter_for_necr.cxx
         stir_timings.cxx
-		create_empty_interfile_binary.cxx
+        create_empty_interfile_binary.cxx
 )
 
 if (HAVE_ITK)

src/utilities/create_empty_interfile_binary.cxx

@@ -15,6 +15,10 @@
   be used for collaborative development. If one user experiences issues with STIR,
   the only the header file needs to be transferred to the other user for debugging issues.
 
+  This utility will not overwrite existing binary files.
+
+  Currently, requires header file to have a valid extension (.hs or .hv) to determine the type of data.
+
   \author Robert Twyman
 */
 /*
@@ -28,10 +32,85 @@
 
 #include <iostream>
 #include <fstream>
+#include "stir/ProjData.h"
+#include "stir/ProjDataInterfile.h"
+#include "stir/IO/InterfileHeader.h"
+#include "stir/Succeeded.h"
+#include "stir/IO/interfile.h"
+#include "stir/utilities.h"
+#include "stir/VoxelsOnCartesianGrid.h"
+#include "stir/IO/write_to_file.h"
+
+using namespace stir;
+
+Succeeded
+create_empty_projdata_interfile_binary(char* header_file, float fill_value)
+{
+  InterfilePDFSHeader hdr;
+
+  if (!hdr.parse(header_file))
+    return Succeeded::no;
+
+  if (std::ifstream(hdr.data_file_name).good())
+    {
+      warning("create_empty_projdata_interfile_binary: data file already exists");
+      return Succeeded::no;
+    }
+  shared_ptr<ProjData> proj_data_sptr(new ProjDataInterfile(hdr.get_exam_info_sptr(), hdr.data_info_sptr, hdr.data_file_name));
+  proj_data_sptr->fill(fill_value);
+  return Succeeded::yes;
+}
+
+Succeeded
+create_empty_discretised_density_interfile_binary(char* header_file, float fill_value)
+{
+  InterfileImageHeader hdr;
+  if (!hdr.parse(header_file))
+    return Succeeded::no;
+
+  if (std::ifstream(hdr.data_file_name).good())
+    {
+      warning("create_empty_discretised_density_interfile_binary: data file already exists");
+      return Succeeded::no;
+    }
+
+  VoxelsOnCartesianGrid<float> density = create_image_from_header(hdr);
+  density.fill(fill_value);
+  write_to_file(hdr.data_file_name, density);
+  return Succeeded::yes;
+}
 
-#ifndef STIR_NO_NAMESPACES
-using std::cerr;
-#endif
+Succeeded
+process_header_file(char* header_file, float fill_value)
+{
+  // Verify that the header file is valid and can be parsed
+  if (!std::ifstream(std::string(header_file).c_str()).is_open())
+    error("read_interfile_image: couldn't open file %s\n", header_file);
+
+  if (strlen(header_file) < 3)
+    {
+      std::cerr << "Warning: Header file name is too short.\n";
+      return Succeeded::no;
+    }
+
+  std::string header_file_ext = std::string(header_file).substr(strlen(header_file) - 3);
+  if (header_file_ext == ".hs")
+    {
+      std::cerr << "File extension indicates a ProjData. Creating empty ProjData with uniform value: " << fill_value << '\n';
+      return create_empty_projdata_interfile_binary(header_file, fill_value);
+    }
+  else if (header_file_ext == ".hv")
+    {
+      std::cerr << "File extension indicates a DiscretisedDensity. Creating empty DiscretisedDensity with uniform value: "
+                << fill_value << '\n';
+      return create_empty_discretised_density_interfile_binary(header_file, fill_value);
+    }
+  else
+    {
+      std::cerr << "Error: Invalid header file extension: " << header_file_ext << '\n';
+      return Succeeded::no;
+    }
+}
 
 int
 main(int argc, char* argv[])
@@ -63,8 +142,15 @@ main(int argc, char* argv[])
         }
     }
 
-  std::cerr << "Valid configuration: " << argv[1] << " " << fill_value << '\n';
-
+  std::cerr << "Valid configuration:\n\tHeader File:\t" << argv[1] << "\n\tFill Value:\t\t" << fill_value << '\n';
+
+  if (process_header_file(argv[1], fill_value) == Succeeded::no)
+    {
+      std::cerr << "Error: Failed to create empty interfile binary file.\n";
+      return EXIT_FAILURE;
+    }
+  std::cerr << "Successfully created empty interfile binary file:\n\tHeader File:\t" << argv[1] << "\n\tFill Value:\t\t"
+            << fill_value << '\n';
 
   return EXIT_SUCCESS;
 }