Skip to content

Commit

Permalink
Merge pull request #64 from ornlneutronimaging/standalone_spectra_out
Browse files Browse the repository at this point in the history 
Standalone spectra out
  • Loading branch information
@KedoKudo
KedoKudo authored Sep 10, 2024
2 parents 02818a4 + 6c73b1e commit 2eccd1c
Show file tree
Hide file tree
Showing 4 changed files with 124 additions and 32 deletions.
5 changes: 5 additions & 0 deletions sophiread/SophireadCLI/include/sophiread_core.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,4 +46,9 @@ void updateTOFImages(
std::vector<std::vector<std::vector<unsigned int>>>& tof_images,
const TPX3& batch, double super_resolution,
const std::vector<double>& tof_bin_edges, const std::string& mode);
std::vector<uint64_t> calculateSpectralCounts(
const std::vector<std::vector<std::vector<unsigned int>>>& tof_images);
void writeSpectralFile(const std::string& filename,
const std::vector<uint64_t>& spectral_counts,
const std::vector<double>& tof_bin_edges);
} // namespace sophiread
39 changes: 29 additions & 10 deletions sophiread/SophireadCLI/src/sophiread.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ struct ProgramOptions {
std::string output_tof_imaging;
std::string tof_filename_base = "tof_image";
std::string tof_mode = "neutron";
std::string spectra_filen = "Spectra";
size_t chunk_size = 5ULL * 1024 * 1024 * 1024; // Default 5GB
bool debug_logging = false;
bool verbose = false;
Expand Down Expand Up @@ -68,6 +69,7 @@ void print_usage(const char* program_name) {
spdlog::info(
" -m <tof_mode> TOF mode: 'hit' or 'neutron' (default: "
"neutron)");
spdlog::info(" -s <spectra_filename> Output filename for spectra");
spdlog::info(" -c <chunk_size> Chunk size in MB (default: 5120)");
spdlog::info(" -d Enable debug logging");
spdlog::info(" -v Enable verbose logging");
Expand All @@ -83,7 +85,7 @@ ProgramOptions parse_arguments(int argc, char* argv[]) {
ProgramOptions options;
int opt;

while ((opt = getopt(argc, argv, "i:H:E:u:T:f:m:c:dv")) != -1) {
while ((opt = getopt(argc, argv, "i:H:E:u:T:f:m:s:c:dv")) != -1) {
switch (opt) {
case 'i':
options.input_tpx3 = optarg;
Expand All @@ -106,6 +108,9 @@ ProgramOptions parse_arguments(int argc, char* argv[]) {
case 'm':
options.tof_mode = optarg;
break;
case 's':
options.spectra_filen = optarg;
break;
case 'c':
options.chunk_size = static_cast<size_t>(std::stoull(optarg)) * 1024 *
1024; // Convert MB to bytes
Expand Down Expand Up @@ -210,7 +215,9 @@ int main(int argc, char* argv[]) {

// Initialize TOF images if needed
std::vector<std::vector<std::vector<unsigned int>>> tof_images;
if (!options.output_tof_imaging.empty()) {
const bool needs_tof_images =
!options.output_tof_imaging.empty() || !options.spectra_filen.empty();
if (needs_tof_images) {
tof_images = sophiread::initializeTOFImages(config->getSuperResolution(),
config->getTOFBinEdges());
}
Expand Down Expand Up @@ -258,7 +265,7 @@ int main(int argc, char* argv[]) {
}

// Update TOF images
if (!options.output_tof_imaging.empty()) {
if (needs_tof_images) {
sophiread::updateTOFImages(
tof_images, batch, config->getSuperResolution(),
config->getTOFBinEdges(), options.tof_mode);
Expand Down Expand Up @@ -293,13 +300,6 @@ int main(int argc, char* argv[]) {
eventsFile.close();
}

// Save TOF images if needed
if (!options.output_tof_imaging.empty()) {
sophiread::timedSaveTOFImagingToTIFF(options.output_tof_imaging,
tof_images, config->getTOFBinEdges(),
options.tof_filename_base);
}

auto end = std::chrono::high_resolution_clock::now();
auto elapsed =
std::chrono::duration_cast<std::chrono::microseconds>(end - start)
Expand All @@ -309,6 +309,25 @@ int main(int argc, char* argv[]) {
spdlog::info("Total hits: {}", totalHits);
spdlog::info("Total neutrons: {}", totalNeutrons);

// Save TOF images if needed
if (!options.output_tof_imaging.empty()) {
spdlog::info("Saving TOF imaging to TIFF: {}",
options.output_tof_imaging);
sophiread::timedSaveTOFImagingToTIFF(options.output_tof_imaging,
tof_images, config->getTOFBinEdges(),
options.tof_filename_base);
}

// Save spectra if needed
if (!options.spectra_filen.empty()) {
spdlog::info("Saving spectra to file: {}", options.spectra_filen);
std::string spectral_filename = options.spectra_filen + ".txt";
std::vector<uint64_t> spectral_counts =
sophiread::calculateSpectralCounts(tof_images);
sophiread::writeSpectralFile(spectral_filename, spectral_counts,
config->getTOFBinEdges());
}

} catch (const std::exception& e) {
spdlog::error("Error: {}", e.what());
return 1;
Expand Down
55 changes: 33 additions & 22 deletions sophiread/SophireadCLI/src/sophiread_core.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -365,10 +365,7 @@ void timedSaveTOFImagingToTIFF(
spdlog::info("Created output directory: {}", out_tof_imaging);
}

// 2. Initialize vector for spectral data
std::vector<uint64_t> spectral_counts(tof_images.size(), 0);

// 3. Iterate through each TOF bin and save TIFF files
// 2. Iterate through each TOF bin and save TIFF files
tbb::parallel_for(
tbb::blocked_range<size_t>(0, tof_images.size()),
[&](const tbb::blocked_range<size_t> &range) {
Expand Down Expand Up @@ -439,39 +436,53 @@ void timedSaveTOFImagingToTIFF(
} else {
spdlog::error("Failed to open TIFF file for writing: {}", filename);
}

// Accumulate counts for spectral file
spectral_counts[bin] = std::accumulate(
accumulated_image.cbegin(), accumulated_image.cend(),
static_cast<uint64_t>(0), [](const auto sum, const auto &row) {
return sum + std::accumulate(row.cbegin(), row.cend(), 0ULL);
});
}
});

// 3. Calculate spectral counts
std::vector<uint64_t> spectral_counts = calculateSpectralCounts(tof_images);
// 4. Write spectral file
std::string spectral_filename =
fmt::format("{}/{}_Spectra.txt", out_tof_imaging, tof_filename_base);
// Write spectral data to file
std::ofstream spectral_file(spectral_filename);
writeSpectralFile(spectral_filename, spectral_counts, tof_bin_edges);

auto end = std::chrono::high_resolution_clock::now();
auto duration =
std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
spdlog::info("TIFF and spectra file writing completed in {} ms",
duration.count());
}

std::vector<uint64_t> calculateSpectralCounts(
const std::vector<std::vector<std::vector<unsigned int>>> &tof_images) {
std::vector<uint64_t> spectral_counts(tof_images.size(), 0);

for (size_t bin = 0; bin < tof_images.size(); ++bin) {
spectral_counts[bin] = std::accumulate(
tof_images[bin].cbegin(), tof_images[bin].cend(),
static_cast<uint64_t>(0), [](const auto sum, const auto &row) {
return sum + std::accumulate(row.cbegin(), row.cend(), 0ULL);
});
}

return spectral_counts;
}

void writeSpectralFile(const std::string &filename,
const std::vector<uint64_t> &spectral_counts,
const std::vector<double> &tof_bin_edges) {
std::ofstream spectral_file(filename);
if (spectral_file.is_open()) {
spectral_file << "shutter_time,counts\n";
for (size_t bin = 0; bin < tof_bin_edges.size() - 1; ++bin) {
spectral_file << tof_bin_edges[bin + 1] << "," << spectral_counts[bin]
<< "\n";
}
spectral_file.close();
spdlog::info("Wrote spectral file: {}", spectral_filename);
spdlog::info("Wrote spectral file: {}", filename);
} else {
spdlog::error("Failed to open spectra file for writing: {}",
spectral_filename);
spdlog::error("Failed to open spectra file for writing: {}", filename);
}

auto end = std::chrono::high_resolution_clock::now();
auto duration =
std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
spdlog::info("TIFF and spectra file writing completed in {} ms",
duration.count());
}

} // namespace sophiread
57 changes: 57 additions & 0 deletions sophiread/SophireadCLI/tests/test_sophiread_core.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,23 @@ TEST_F(SophireadCoreTest, TimedSaveTOFImagingToTIFF) {
EXPECT_TRUE(std::filesystem::exists("test_tof/test_bin_0002.tiff"));
EXPECT_TRUE(std::filesystem::exists("test_tof/test_bin_0003.tiff"));
EXPECT_TRUE(std::filesystem::exists("test_tof/test_Spectra.txt"));

// Check spectra file contents
std::ifstream spectra_file("test_tof/test_Spectra.txt");
std::string line;
std::vector<std::string> lines;
while (std::getline(spectra_file, line)) {
lines.push_back(line);
}
// close the file
spectra_file.close();

EXPECT_EQ(lines.size(), 4); // Header + 3 data lines
EXPECT_EQ(lines[0], "shutter_time,counts");
EXPECT_EQ(lines[1], "0.1,100"); // 10 * 10 * 1 = 100 counts for each bin
EXPECT_EQ(lines[2], "0.2,100");
EXPECT_EQ(lines[3], "0.3,100");

std::filesystem::remove_all("test_tof");
}

Expand Down Expand Up @@ -180,6 +197,46 @@ TEST_F(SophireadCoreTest, UpdateTOFImages) {
EXPECT_TRUE(updated);
}

TEST_F(SophireadCoreTest, CalculateSpectralCounts) {
std::vector<std::vector<std::vector<unsigned int>>> tof_images(
3, std::vector<std::vector<unsigned int>>(
5, std::vector<unsigned int>(
5, 2))); // 3 bins, 5x5 images, all values 2

auto spectral_counts = sophiread::calculateSpectralCounts(tof_images);

EXPECT_EQ(spectral_counts.size(), 3); // 3 bins
for (const auto& count : spectral_counts) {
EXPECT_EQ(count, 50); // 5 * 5 * 2 = 50 for each bin
}
}

TEST_F(SophireadCoreTest, WriteSpectralFile) {
std::vector<uint64_t> spectral_counts = {10, 20, 30};
std::vector<double> tof_bin_edges = {0.0, 0.1, 0.2, 0.3};
std::string filename = "test_spectra.txt";

sophiread::writeSpectralFile(filename, spectral_counts, tof_bin_edges);

EXPECT_TRUE(std::filesystem::exists(filename));

// Read and check file contents
std::ifstream file(filename);
std::string line;
std::vector<std::string> lines;
while (std::getline(file, line)) {
lines.push_back(line);
}

EXPECT_EQ(lines.size(), 4); // Header + 3 data lines
EXPECT_EQ(lines[0], "shutter_time,counts");
EXPECT_EQ(lines[1], "0.1,10");
EXPECT_EQ(lines[2], "0.2,20");
EXPECT_EQ(lines[3], "0.3,30");

std::filesystem::remove(filename);
}

int main(int argc, char** argv) {
::testing::InitGoogleTest(&argc, argv);
return RUN_ALL_TESTS();
Expand Down

0 comments on commit 2eccd1c

Please sign in to comment.