Skip to content

Commit

Permalink
Merge pull request #123 from aglowacki/master
Browse files Browse the repository at this point in the history 
Updated netcdf to be able to load 8 element detectors
  • Loading branch information
Arthur Glowacki authored Oct 18, 2021
2 parents 36226af + bf24d2e commit ae3f5cf
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 48 deletions.
2 changes: 1 addition & 1 deletion src/core/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,7 +401,7 @@ int main(int argc, char *argv[])
}
else
{
for (size_t det = 0; det < 4; det++)
for (size_t det = 0; det < 8; det++)
{
analysis_job.detector_num_arr.push_back(det);
}
Expand Down
23 changes: 20 additions & 3 deletions src/io/file/hl_file_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1004,7 +1004,11 @@ bool load_spectra_volume(std::string dataset_directory,
full_filename = dataset_directory + "flyXRF"+ DIR_END_CHAR + tmp_dataset_file + file_middle + std::to_string(i) + ".nc";
//todo: add verbose option
//logI<<"Loading file "<<full_filename<<"\n";
io::file::NetCDF_IO::inst()->load_spectra_line(full_filename, detector_num, &(*spectra_volume)[i]);
size_t spec_size = io::file::NetCDF_IO::inst()->load_spectra_line(full_filename, detector_num, &(*spectra_volume)[i]);
if (detector_num > 3 && spec_size == -1) // this netcdf file only has 4 element detectors
{
return false;
}
}
}
else
Expand Down Expand Up @@ -1033,6 +1037,11 @@ bool load_spectra_volume(std::string dataset_directory,
full_filename = dataset_directory + "flyXRF"+ DIR_END_CHAR + bnp_netcdf_base_name + row_idx_str_full + ".nc";
size_t prev_size = 0;
size_t spec_size = io::file::NetCDF_IO::inst()->load_spectra_line(full_filename, detector_num, &(*spectra_volume)[i]);
//
if (detector_num > 3 && spec_size == -1) // this netcdf file only has 4 element detectors
{
return false;
}
//if we failed to load and it isn't the first row, copy the previous one
if (i > 0)
{
Expand Down Expand Up @@ -1330,7 +1339,11 @@ bool load_and_integrate_spectra_volume(std::string dataset_directory,
{
full_filename = dataset_directory + "flyXRF"+ DIR_END_CHAR + tmp_dataset_file + file_middle + std::to_string(i) + ".nc";
//logI<<"Loading file "<<full_filename<<"\n";
io::file::NetCDF_IO::inst()->load_spectra_line_integrated(full_filename, detector_num, dims[1], integrated_spectra);
size_t spec_size = io::file::NetCDF_IO::inst()->load_spectra_line_integrated(full_filename, detector_num, dims[1], integrated_spectra);
if (detector_num > 3 && spec_size == -1) // this netcdf file only has 4 element detectors
{
return false;
}
}
}
else
Expand All @@ -1357,7 +1370,11 @@ bool load_and_integrate_spectra_volume(std::string dataset_directory,
}
row_idx_str_full += row_idx_str;
full_filename = dataset_directory + "flyXRF"+ DIR_END_CHAR + bnp_netcdf_base_name + row_idx_str_full + ".nc";
io::file::NetCDF_IO::inst()->load_spectra_line_integrated(full_filename, detector_num, dims[1], integrated_spectra);
size_t spec_size = io::file::NetCDF_IO::inst()->load_spectra_line_integrated(full_filename, detector_num, dims[1], integrated_spectra);
if (detector_num > 3 && spec_size == -1) // this netcdf file only has 4 element detectors
{
return false;
}
}
}
else
Expand Down
139 changes: 95 additions & 44 deletions src/io/file/netcdf_io.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ NetCDF_IO* NetCDF_IO::_this_inst(nullptr);
#define INPUT_COUNTS_OFFSET 36
#define OUTPUT_COUNTS_OFFSET 38

#define MAX_NUM_SUPPORTED_DETECOTRS 4
#define MAX_NUM_SUPPORTED_DETECOTRS_PER_COL 4

//-----------------------------------------------------------------------------

Expand Down Expand Up @@ -105,7 +105,6 @@ size_t NetCDF_IO::load_spectra_line(std::string path, size_t detector, data_stru
ptrdiff_t stride[] = {1, 1, 1};
real_t data_in[1][1][10000];
size_t spectra_size;

nc_type rh_type;
int rh_ndims;
int rh_dimids[NC_MAX_VAR_DIMS] = {0};
Expand Down Expand Up @@ -148,6 +147,17 @@ size_t NetCDF_IO::load_spectra_line(std::string path, size_t detector, data_stru
}
}

if (detector > 3)
{
if (dim2size[1] != 2)
{
logE << "NetCDF dims: [" << dim2size[0] <<"]["<< dim2size[1] <<"]["<< dim2size[2] <<"] needs to be [x][2][x] for detector "<<detector<<" " << path << "\n";
nc_close(ncid);
return -1;
}
start[1] = 1;
}

if( (retval = nc_get_vars_real(ncid, varid, start, count, stride, &data_in[0][0][0]) ) != 0)
{
logE<< path << " :: " << nc_strerror(retval)<<"\n";
Expand All @@ -157,7 +167,7 @@ size_t NetCDF_IO::load_spectra_line(std::string path, size_t detector, data_stru

if (data_in[0][0][0] != 21930 || data_in[0][0][1] != -21931)
{
logE<<"NetCDF header not found! Stopping load : "<<path<<"\n";
logE<<"NetCDF header [0][0][0] not found! Stopping load : "<<path<<"\n";
nc_close(ncid);
return 0;
}
Expand All @@ -176,6 +186,11 @@ size_t NetCDF_IO::load_spectra_line(std::string path, size_t detector, data_stru
count[2] = header_size;
size_t j=0;

if (detector > 3)
{
detector -= MAX_NUM_SUPPORTED_DETECOTRS_PER_COL; // 4,5,6,7 = 0,1,2,3
}

for(; j<spec_line->size(); j++)
{
(*spec_line)[j].resize(spectra_size); // should be renames to resize
Expand All @@ -202,6 +217,7 @@ size_t NetCDF_IO::load_spectra_line(std::string path, size_t detector, data_stru
return j;
}


unsigned short i1 = data_in[0][0][ELAPSED_LIVETIME_OFFSET+(detector*8)];
unsigned short i2 = data_in[0][0][ELAPSED_LIVETIME_OFFSET+(detector*8)+1];
unsigned int ii = i1 | i2<<16;
Expand Down Expand Up @@ -296,30 +312,25 @@ bool NetCDF_IO::load_spectra_line_with_callback(std::string path,
size_t header_size = 256;
int ncid, varid, retval;
size_t start[] = {0, 0, 0};
size_t count[] = {1, 1, header_size};
size_t count[] = {1, 2, header_size};
ptrdiff_t stride[] = {1, 1, 1};
real_t data_in[1][1][10000];
//real_t data_in[1][2][18000];
real_t *data_in;
size_t spectra_size;
size_t num_detectors = detector_num_arr.size();
int dataidx = 0;
real_t elapsed_livetime = 0.;
real_t elapsed_realtime = 0.;
real_t input_counts = 0.;
real_t output_counts = 0.;

nc_type rh_type;
int rh_ndims;
int rh_dimids[NC_MAX_VAR_DIMS] = {0};
int rh_natts;

real_t elapsed_livetime[MAX_NUM_SUPPORTED_DETECOTRS];
real_t elapsed_realtime[MAX_NUM_SUPPORTED_DETECOTRS];
real_t input_counts[MAX_NUM_SUPPORTED_DETECOTRS];
real_t output_counts[MAX_NUM_SUPPORTED_DETECOTRS];

size_t dim2size[NC_MAX_VAR_DIMS] = {0};

if( num_detectors > 4)
{
logE<<"Max detectors supported is 4, requesting :"<< num_detectors<<"\n";
return false;
}

if( (retval = nc_open(path.c_str(), NC_NOWRITE, &ncid)) != 0 )
{
logE<< path << " :: " << nc_strerror(retval)<<"\n";
Expand Down Expand Up @@ -350,43 +361,52 @@ bool NetCDF_IO::load_spectra_line_with_callback(std::string path,
}
}

data_in = new real_t[dim2size[0] * dim2size[1] * dim2size[2]];
count[2] = dim2size[2];
//read in last col sector to get total number of cols
//start[0] = dim2size[0] - 1;
if( (retval = nc_get_vars_real(ncid, varid, start, count, stride, &data_in[0][0][0]) ) != 0)
if( (retval = nc_get_vars_real(ncid, varid, start, count, stride, &data_in[0]) ) != 0)
{
delete[] data_in;
logE<< path << " :: " << nc_strerror(retval)<<"\n";
nc_close(ncid);
return false;
}

if (data_in[0][0][0] != 21930 || data_in[0][0][1] != -21931)
if (data_in[0] != 21930 || data_in[1] != -21931)
{
delete[] data_in;
logE<<"NetCDF header not found! Stopping load : "<<path<<"\n";
nc_close(ncid);
return false;
}

//can't read from file because it can change inbetween rows ...
//max_cols = (124 * (dim2size[0] - 1) ) + data_in[0][0][8];
header_size = data_in[0][0][2];
header_size = data_in[2];
//num_cols = data_in[][0][8]; //sum all across the first dim looking at value 8
spectra_size = data_in[0][0][20];
spectra_size = data_in[20];

start[2] += header_size;
count[2] = header_size + (spectra_size * MAX_NUM_SUPPORTED_DETECOTRS_PER_COL); //only 4 element detector supported per col

start[2] += count[2];
count[2] = header_size + (spectra_size * MAX_NUM_SUPPORTED_DETECOTRS); //only 4 element detector supported
//loop through col sectors
for(size_t j = 0; j < max_cols; j++)
{
/*
//read header
if( (retval = nc_get_vars_real(ncid, varid, start, count, stride, &data_in[0][0][0]) ) != 0)
{
logE<< path << " :: " << nc_strerror(retval)<<"\n";
nc_close(ncid);
return false;
}
*/
int midx = start[2];

if (data_in[0][0][0] != 13260 || data_in[0][0][1] != -13261)
if (data_in[midx] != 13260 || data_in[midx + 1] != -13261)
{
delete[] data_in;
if(j < max_cols -2)
{
logE<<"NetCDF sub header not found! Stopping load at Col: "<<j<<" path :"<<path<<"\n";
Expand All @@ -401,55 +421,68 @@ bool NetCDF_IO::load_spectra_line_with_callback(std::string path,

for(size_t detector_num : detector_num_arr)
{
unsigned short i1 = data_in[0][0][ELAPSED_LIVETIME_OFFSET+(detector_num*8)];
unsigned short i2 = data_in[0][0][ELAPSED_LIVETIME_OFFSET+(detector_num*8)+1];

if (detector_num > 3)
{
dataidx = 1;
detector_num -= MAX_NUM_SUPPORTED_DETECOTRS_PER_COL;
}
else
{
dataidx = 0;
}

midx = (dataidx * dim2size[2]) + start[2];

unsigned short i1 = data_in[midx+(ELAPSED_LIVETIME_OFFSET+(detector_num*8))];
unsigned short i2 = data_in[midx + (ELAPSED_LIVETIME_OFFSET+(detector_num*8)+1)];
unsigned int ii = i1 | i2<<16;
elapsed_livetime[detector_num] = ((float)ii) * 320e-9f; // need to multiply by this value becuase of the way it is saved
if(elapsed_livetime[detector_num] == 0)
elapsed_livetime = ((float)ii) * 320e-9f; // need to multiply by this value becuase of the way it is saved
if(elapsed_livetime == 0)
{
if(j < max_cols-2) // usually the last two are missing which spams the log ouput.
{
logW<<"Reading in elapsed lifetime for Col:"<<j<<" is 0. Setting it to 1.0 " << path <<"\n";
elapsed_livetime[detector_num] = 1.0;
elapsed_livetime = 1.0;
}
}

i1 = data_in[0][0][ELAPSED_REALTIME_OFFSET+(detector_num*8)];
i2 = data_in[0][0][ELAPSED_REALTIME_OFFSET+(detector_num*8)+1];
i1 = data_in[midx + (ELAPSED_REALTIME_OFFSET+(detector_num*8))];
i2 = data_in[midx + (ELAPSED_REALTIME_OFFSET+(detector_num*8)+1)];
ii = i1 | i2<<16;
elapsed_realtime[detector_num] = ((float)ii) * 320e-9f; // need to multiply by this value becuase of the way it is saved
if(elapsed_realtime[detector_num] == 0)
elapsed_realtime = ((float)ii) * 320e-9f; // need to multiply by this value becuase of the way it is saved
if(elapsed_realtime == 0)
{
if(j < max_cols-2) // usually the last two are missing which spams the log ouput.
{
logW<<"Reading in elapsed realtime for Col:"<<j<<" is 0. Setting it to 1.0 "<<path<<"\n";
elapsed_realtime[detector_num] = 1.0;
elapsed_realtime = 1.0;
}
}

i1 = data_in[0][0][INPUT_COUNTS_OFFSET+(detector_num*8)];
i2 = data_in[0][0][INPUT_COUNTS_OFFSET+(detector_num*8)+1];
i1 = data_in[midx + (INPUT_COUNTS_OFFSET+(detector_num*8))];
i2 = data_in[midx + (INPUT_COUNTS_OFFSET+(detector_num*8)+1)];
ii = i1 | i2<<16;
input_counts[detector_num] = ((float)ii) / elapsed_livetime[detector_num];
input_counts = ((float)ii) / elapsed_livetime;

i1 = data_in[0][0][OUTPUT_COUNTS_OFFSET+(detector_num*8)];
i2 = data_in[0][0][OUTPUT_COUNTS_OFFSET+(detector_num*8)+1];
i1 = data_in[midx + (OUTPUT_COUNTS_OFFSET+(detector_num*8))];
i2 = data_in[midx + (OUTPUT_COUNTS_OFFSET+(detector_num*8)+1)];
ii = i1 | i2<<16;
output_counts[detector_num] = ((float)ii) / elapsed_realtime[detector_num];
output_counts = ((float)ii) / elapsed_realtime;

data_struct::Spectra * spectra = new data_struct::Spectra(spectra_size);

spectra->elapsed_livetime(elapsed_livetime[detector_num]);
spectra->elapsed_realtime(elapsed_realtime[detector_num]);
spectra->input_counts(input_counts[detector_num]);
spectra->output_counts(output_counts[detector_num]);
spectra->elapsed_livetime(elapsed_livetime);
spectra->elapsed_realtime(elapsed_realtime);
spectra->input_counts(input_counts);
spectra->output_counts(output_counts);
spectra->recalc_elapsed_livetime();

int idx = header_size + (detector_num*spectra_size);

for(size_t k=0; k<spectra_size; k++)
{
(*spectra)[k] = data_in[0][0][idx+k];
(*spectra)[k] = data_in[midx + (idx+k)];
}

callback_fun(row, j, max_rows, max_cols, detector_num, spectra, user_data);
Expand All @@ -465,6 +498,8 @@ bool NetCDF_IO::load_spectra_line_with_callback(std::string path,
}
}

delete[] data_in;

if((retval = nc_close(ncid)) != 0)
{
logE<< path << " :: " << nc_strerror(retval)<<"\n";
Expand Down Expand Up @@ -531,6 +566,17 @@ size_t NetCDF_IO::load_spectra_line_integrated(std::string path, size_t detector
}
}

if (detector > 3)
{
if (dim2size[1] != 2)
{
logE << "NetCDF dims: [" << dim2size[0] << "][" << dim2size[1] << "][" << dim2size[2] << "] needs to be [x][2][x] for detector " << detector << " " << path << "\n";
nc_close(ncid);
return 0;
}
start[1] = 1;
}

if ((retval = nc_get_vars_real(ncid, varid, start, count, stride, &data_in[0][0][0])) != 0)
{
logE << path << " :: " << nc_strerror(retval) << "\n";
Expand Down Expand Up @@ -559,6 +605,11 @@ size_t NetCDF_IO::load_spectra_line_integrated(std::string path, size_t detector
count[2] = header_size;
size_t j = 0;

if (detector > 3)
{
detector -= MAX_NUM_SUPPORTED_DETECOTRS_PER_COL; // 4,5,6,7 = 0,1,2,3
}

for (; j < line_size; j++)
{
//(*spec_line)[j].resize(spectra_size); // should be renames to resize
Expand Down

0 comments on commit ae3f5cf

Please sign in to comment.