diff --git a/cmake/CMakeDebugFlags.cmake b/cmake/CMakeDebugFlags.cmake index d5e8fa220..42022a620 100644 --- a/cmake/CMakeDebugFlags.cmake +++ b/cmake/CMakeDebugFlags.cmake @@ -7,5 +7,5 @@ endif() #-DDATA_ANALYTICS=1 to enable if(${DATA_ANALYTICS}) - add_definitions(-DDATA_ANALYTICS) + add_compile_definitions(DATA_ANALYTICS=${DATA_ANALYTICS}) endif() \ No newline at end of file diff --git a/pytests/test_demo_scenes.py b/pytests/test_demo_scenes.py index be9d7542c..43a090074 100644 --- a/pytests/test_demo_scenes.py +++ b/pytests/test_demo_scenes.py @@ -120,9 +120,9 @@ def test_arbaro_tls_pyh(): def eval_arbaro_tls(dirname): assert (dirname / 'leg000_points.las').exists() - assert abs((dirname / 'leg000_points.las').stat().st_size - 21_656_853) < MAX_DIFFERENCE_BYTES + assert abs((dirname / 'leg000_points.las').stat().st_size - 21_755_791) < MAX_DIFFERENCE_BYTES assert (dirname / 'leg001_points.las').exists() - assert abs((dirname / 'leg001_points.las').stat().st_size - 13_908_757) < MAX_DIFFERENCE_BYTES + assert abs((dirname / 'leg001_points.las').stat().st_size - 13_936_877) < MAX_DIFFERENCE_BYTES with open(dirname / 'leg000_trajectory.txt', 'r') as f: line = f.readline() assert line.startswith('1.0000 25.5000 0.0000') @@ -231,7 +231,7 @@ def test_interpolated_traj_pyh(): def eval_interpolated_traj(dirname): assert (dirname / 'leg000_points.laz').exists() assert (dirname / 'leg000_trajectory.txt').exists() - assert abs((dirname / 'leg000_points.laz').stat().st_size - 873_554) < MAX_DIFFERENCE_BYTES + assert abs((dirname / 'leg000_points.laz').stat().st_size - 875_054) < MAX_DIFFERENCE_BYTES with open(dirname / 'leg000_trajectory.txt', 'r') as f: for _ in range(3): next(f) diff --git a/pytests/test_gpsStartTimeFlag.py b/pytests/test_gpsStartTimeFlag.py index 684491810..4b4e1d114 100644 --- a/pytests/test_gpsStartTimeFlag.py +++ b/pytests/test_gpsStartTimeFlag.py @@ -66,7 +66,7 @@ def test_gpsStartTimeFlag_exe(): r3_sum = sha256sum(r3 / 'leg000_points.xyz') assert r2_sum == r3_sum assert r2_sum == 'e1daca137b066eca4bd29b62f3cb5ecde4bcc41650f2cdfc44110c87b97f23dc' or \ - r2_sum == '59b9d09d6e2a6b8a8c25112f5b06fc08430c4058216efeea6ec7324dff0c7797' # linux checksum + r2_sum == 'e5ff320c318bd8cb87865ff8a3fd5d8835ca7a66736cc7469e59ad2c365a3a31' # linux checksum assert r1_sum != r2_sum if DELETE_FILES_AFTER: @@ -95,7 +95,7 @@ def test_gpsStartTimeFlag_pyh(): r3_sum = sha256sum(r3 / 'leg000_points.xyz') assert r2_sum == r3_sum assert r2_sum == 'e1daca137b066eca4bd29b62f3cb5ecde4bcc41650f2cdfc44110c87b97f23dc' or \ - r2_sum == '59b9d09d6e2a6b8a8c25112f5b06fc08430c4058216efeea6ec7324dff0c7797' # linux checksum + r2_sum == 'e5ff320c318bd8cb87865ff8a3fd5d8835ca7a66736cc7469e59ad2c365a3a31' # linux checksum assert r1_sum != r2_sum if DELETE_FILES_AFTER: @@ -104,4 +104,4 @@ def test_gpsStartTimeFlag_pyh(): shutil.rmtree(r2) shutil.rmtree(r3) except Exception as e: - print(f"Error cleaning up: {e}") \ No newline at end of file + print(f"Error cleaning up: {e}") diff --git a/pytests/test_pyhelios.py b/pytests/test_pyhelios.py index d7b39b0ed..19e660e0b 100644 --- a/pytests/test_pyhelios.py +++ b/pytests/test_pyhelios.py @@ -278,7 +278,7 @@ def test_create_survey(): output = simB.join() meas, traj = pyhelios.outputToNumpy(output) # check length of output - assert meas.shape == (10397, 17) + assert meas.shape == (10470, 17) assert traj.shape == (6670, 7) # compare individual points np.testing.assert_allclose(meas[100, :3], np.array([83.32, -66.44204, 0.03114649])) @@ -373,7 +373,7 @@ def test_output(export_to_file): measurements_array, trajectory_array = pyhelios.outputToNumpy(output) np.testing.assert_allclose(measurements_array[0, :3], np.array([474500.3, 5473530.0, 106.0988]), rtol=0.000001) - assert measurements_array.shape == (2433, 17) + assert measurements_array.shape == (2435, 17) assert trajectory_array.shape == (9, 7) if export_to_file: assert Path(output.outpath).parent.parent == Path(WORKING_DIR) / "output" / "als_hd_demo" diff --git a/scripts/debug/fullwave_plotter.py b/scripts/debug/fullwave_plotter.py old mode 100644 new mode 100755 diff --git a/scripts/debug/hda_diff_report.py b/scripts/debug/hda_diff_report.py old mode 100644 new mode 100755 diff --git a/scripts/debug/hda_pulse_calc_intensity_csv_to_laz.sh b/scripts/debug/hda_pulse_calc_intensity_csv_to_laz.sh new file mode 100755 index 000000000..46fca7b98 --- /dev/null +++ b/scripts/debug/hda_pulse_calc_intensity_csv_to_laz.sh @@ -0,0 +1,52 @@ +#!/bin/bash + +# ######################################################################## +# # +# AUTHOR : Alberto M. Esmorís Pena # +# # +# Script to transform an input CSV file, which rows are pulse records # +# written by the HDA_PulseRecorder component of HELIOS++, to a LAS/LAZ # +# file. # +# # +# ######################################################################## + +# Path or alias to TXT2LAS util +TXT2LAS='txt2las' + +# Function to show help +function show_help { + cat << EOF + +hda_pulse_calc_intensity_csv_to_laz.sh + + Argument 1: Path to the input CSV file + + Argument 2: Path to the output LAZ file + + +EOF +} + +# Handle input arguments +if [[ $# -ne 2 ]]; then + echo -e 'ERROR: Exactly two arguments must be provided\n' + show_help + exit 1 +fi + + +echo -e "Transforming \"$1\" to \"$2\" ..." +mkdir -p $(dirname "$2") +${TXT2LAS} -i "$1" \ + -set_version 1.4 -parse xyz01234567 -rescale 1e-5 1e-5 1e-5 \ + -add_attribute 10 'incidence_angle_rad' 'incidence_angle_rad' \ + -add_attribute 10 'target_range_m' 'target_range_m' \ + -add_attribute 10 'target_area_m2' 'target_area_m2' \ + -add_attribute 10 'radius_m' 'radius_m' \ + -add_attribute 10 'bdrf' 'bdrf' \ + -add_attribute 10 'cross_section_m2' 'cross_section_m2' \ + -add_attribute 10 'received_power' 'received_power' \ + -add_attribute 2 'captured' 'captured' \ + -o "$2" + + diff --git a/scripts/debug/hda_pulse_records_plotter.py b/scripts/debug/hda_pulse_records_plotter.py new file mode 100755 index 000000000..e96252a40 --- /dev/null +++ b/scripts/debug/hda_pulse_records_plotter.py @@ -0,0 +1,1025 @@ +import sys +import os +import time +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl + + +# --- FUNCTIONS --- # +# ------------------- # +def print_help(): + print( +''' +Input arguments: + 1 -> Path to the first directory containing simulation records + 2 -> Path to the second directory containing simulation records + 3 -> Path to the directory where plots will be stored +''' + ) + + +def parse_args(helpf=print_help): + """Parse input arguments. Raise an exception if not correct arguments were + given""" + if len(sys.argv) == 1: + helpf() + exit(0) + elif len(sys.argv) < 4: + raise Exception( + '{m} arguments were given but 3 are required' + .format(m=len(sys.argv)-1) + ) + dira_path = sys.argv[1] + if not validate_directory(dira_path): + raise Exception( + 'The directory "{d}"\n' + 'was given as the first directory of records, but it is not valid' + ) + dirb_path = sys.argv[2] + if not validate_directory(dirb_path): + raise Exception( + 'The directory "{d}"\n' + 'was given as the second directory of records, but it is not valid' + ) + dirout_path = sys.argv[3] + if not validate_directory(dirout_path): + raise Exception( + 'The directory "{d}"\n' + 'was given as the third directory for plots, but it is not valid' + ) + return { + 'dira_path': dira_path, + 'dirb_path': dirb_path, + 'dirout_path': dirout_path + } + + +def validate_directory(path): + """Check path points to a valid existent directory. + If it does not exist and it cannot be created, then it is not valid""" + if os.path.exists(path): + if os.path.isdir(path): + return True + else: + return False + else: + os.mkdir(path, mode=0o775) + if os.path.exists(path): + return True + return False + + +def read_records(path, sep=','): + """Read all record files contained in the directory pointed by given + path""" + # Read vectorial records + intensity_calc = read_record(os.path.join( + path, 'intensity_calc.csv' + ), sep) + subray_sim = read_record(os.path.join( + path, 'subray_sim.csv' + ), sep) + # Return key-word records + return { + # Intensity calculation records + 'incidence_angle_rad': intensity_calc[:, 3], + 'target_range_m': intensity_calc[:, 4], + 'target_area_m2': intensity_calc[:, 5], + 'radius_m': intensity_calc[:, 6], + 'bdrf': intensity_calc[:, 7], + 'cross_section': intensity_calc[:, 8], + 'received_power': intensity_calc[:, 9], + # Subray simulation records + 'subray_hit': subray_sim[:, 0].astype(bool), + 'radius_step': subray_sim[:, 1], + 'circle_steps': subray_sim[:, 2], + 'circle_step': subray_sim[:, 3], + 'divergence_angle_rad': subray_sim[:, 4], + 'ray_dir_norm': subray_sim[:, 5], + 'subray_dir_norm': subray_sim[:, 6], + 'ray_subray_angle_rad': subray_sim[:, 7], + 'ray_subray_sign_check': subray_sim[:, 8], + 'subray_tmin': subray_sim[:, 9], + 'subray_tmax': subray_sim[:, 10], + 'subray_dir_x': subray_sim[:, 11], + 'subray_dir_y': subray_sim[:, 12], + 'subray_dir_z': subray_sim[:, 13], + 'ray_dir_x': subray_sim[:, 14], + 'ray_dir_y': subray_sim[:, 15], + 'ray_dir_z': subray_sim[:, 16] + } + + +def read_record(path, sep): + """Read given record file + :param path: The path to the record file to be read + :param sep: The separator used in the record file + :return: None if the record file could not be read, the data as a numpy + array otherwise + """ + if os.path.exists(path) and os.path.isfile(path): + return np.loadtxt(path, delimiter=sep) + return None + + +def plot_records(arec, brec, outdir): + """Do plots for each case and export them in given output directory + :param arec: The records of the first case + :param brec: The records of the second case + :param outdir: The output directory where generated plots must be + written + """ + do_incidence_angle_plots(arec, brec, outdir) + do_by_incidence_angle_plots(arec, brec, outdir) + do_subray_hit_plots(arec, brec, outdir) + do_ray_subray_plots(arec, brec, outdir) + do_ray_subray_dir_plots(arec, brec, outdir) + + +def validate_record(key, rec, recid): + """Check that the record with given key is available + :return: True if the record is valid (available data), False otherwise + """ + if key not in rec or rec.get(key, None) is None: + print( + 'Record "{key}" is not available for {recid} records' + .format( + key=key, + recid=recid + ) + ) + return False + return True + + +def init_figure(figsize=(20, 12), constrained_layout=False): + """Initialize a matplotlib's figure context""" + fig = plt.figure( + figsize=figsize, + constrained_layout=constrained_layout + ) + return fig + + +def do_incidence_angle_subplot( + fig, ax, phi, label=None, title=None, xlabel=None, ylabel=None, bins=32, + log=False, relative=False, label_fontsize=16 +): + if title is not None: + ax.set_title(title, fontsize=20) + weights = 100*np.ones_like(phi)/len(phi) if relative else None + hist = ax.hist(phi, bins=bins, label=label, log=log, weights=weights) + ax.axvline(x=np.mean(phi), color='tab:orange', lw=3, label='$\\mu$') + if xlabel is not None: + ax.set_xlabel(xlabel, fontsize=label_fontsize) + if ylabel is not None: + ax.set_ylabel(ylabel, fontsize=label_fontsize) + ax.tick_params(axis='both', which='both', labelsize=14) + ax.legend(loc='upper right', fontsize=14) + ax.grid('both') + ax.set_axisbelow(True) + + +def do_incidence_angle_plots(arec, brec, outdir): + # Validate incidence angle data + if not validate_record('incidence_angle_rad', arec, 'a') or \ + not validate_record('incidence_angle_rad', brec, 'b'): + print('Cannot do incidence angle plots') + return + + # Do the incidence angle plots (rads) + fig = init_figure() # Initialize figure + ax = fig.add_subplot(2, 2, 1) # Initialize phi(a) subplot + do_incidence_angle_subplot( + fig, ax, arec['incidence_angle_rad'], + label='$\\varphi(a)$', + title='A-Incidence angle ($\\varphi$) in rad', + xlabel='$\\varphi(a)$', + ylabel='cases' + ) + ax = fig.add_subplot(2, 2, 3) # Initialize phi(a) log subplot + do_incidence_angle_subplot( + fig, ax, arec['incidence_angle_rad'], + label='$\\varphi(a)$', + title='A-Incidence angle ($\\varphi$) in rad (logarithmic)', + xlabel='$\\varphi(a)$', + ylabel='cases', + log=True + ) + ax = fig.add_subplot(2, 2, 2) # Initialize phi(b) subplot + do_incidence_angle_subplot( + fig, ax, brec['incidence_angle_rad'], + label='$\\varphi(b)$', + title='B-Incidence angle ($\\varphi$) in rad', + xlabel='$\\varphi(b)$', + ylabel='cases' + ) + ax = fig.add_subplot(2, 2, 4) # Initialize phi(a) log subplot + do_incidence_angle_subplot( + fig, ax, brec['incidence_angle_rad'], + label='$\\varphi(b)$', + title='B-Incidence angle ($\\varphi$) in rad (logarithmic)', + xlabel='$\\varphi(b)$', + ylabel='cases', + log=True + ) + fig.tight_layout() + # Save figure to file and remove it from memory + fig.savefig( + os.path.join(outdir, 'incidence_angle_distribution_rad.png') + ) + fig.clear() + plt.close(fig) + # Do the incidence angle plots (degrees) + fig = init_figure() # Initialize figure + ax = fig.add_subplot(2, 2, 1) # Initialize phi(a) subplot + do_incidence_angle_subplot( + fig, ax, arec['incidence_angle_rad']*180/np.pi, + label='$\\varphi(a)$', + title='A-Incidence angle ($\\varphi$) in rad', + xlabel='$\\varphi(a)$', + ylabel='cases' + ) + ax = fig.add_subplot(2, 2, 3) # Initialize phi(a) log subplot + do_incidence_angle_subplot( + fig, ax, arec['incidence_angle_rad']*180/np.pi, + label='$\\varphi(a)$', + title='A-Incidence angle ($\\varphi$) in rad (logarithmic)', + xlabel='$\\varphi(a)$', + ylabel='cases', + log=True + ) + ax = fig.add_subplot(2, 2, 2) # Initialize phi(b) subplot + do_incidence_angle_subplot( + fig, ax, brec['incidence_angle_rad']*180/np.pi, + label='$\\varphi(b)$', + title='B-Incidence angle ($\\varphi$) in rad', + xlabel='$\\varphi(b)$', + ylabel='cases' + ) + ax = fig.add_subplot(2, 2, 4) # Initialize phi(a) log subplot + do_incidence_angle_subplot( + fig, ax, brec['incidence_angle_rad']*180/np.pi, + label='$\\varphi(b)$', + title='B-Incidence angle ($\\varphi$) in rad (logarithmic)', + xlabel='$\\varphi(b)$', + ylabel='cases', + log=True + ) + fig.tight_layout() + # Save figure to file and remove it from memory + fig.savefig( + os.path.join(outdir, 'incidence_angle_distribution_degrees.png') + ) + fig.clear() + plt.close(fig) + + +def do_y_by_x_subplot( + fig, ax, x, y, title=None, xlabel=None, ylabel=None, + color='black' +): + if title is not None: + ax.set_title(title, fontsize=14) + ax.scatter(x, y, c=color, s=8) + if xlabel is not None: + ax.set_xlabel(xlabel, fontsize=12) + if ylabel is not None: + ax.set_ylabel(ylabel, fontsize=12) + ax.tick_params(axis='both', which='both', labelsize=12) + ax.grid('both') + ax.set_axisbelow(True) + + +def do_by_incidence_angle_plots(arec, brec, outdir): + # Validate classification calculation data + if not validate_record('incidence_angle_rad', arec, 'a') or \ + not validate_record('target_range_m', arec, 'a') or \ + not validate_record('target_area_m2', arec, 'a') or \ + not validate_record('radius_m', arec, 'a') or \ + not validate_record('bdrf', arec, 'a') or \ + not validate_record('cross_section', arec, 'a') or \ + not validate_record('received_power', arec, 'a') or \ + not validate_record('incidence_angle_rad', brec, 'b') or \ + not validate_record('target_range_m', brec, 'b') or \ + not validate_record('target_area_m2', brec, 'b') or \ + not validate_record('radius_m', brec, 'b') or \ + not validate_record('bdrf', brec, 'b') or \ + not validate_record('cross_section', brec, 'b') or \ + not validate_record('received_power', brec, 'b'): + print('Cannot do by incidence angle plots') + return + # Plot by radians + _do_by_incidence_angle_plots( + arec['incidence_angle_rad'], + brec['incidence_angle_rad'], + arec, + brec, + outdir, + fname='plot_by_incidence_angle_rad.png' + ) + # Plot by degrees + _do_by_incidence_angle_plots( + arec['incidence_angle_rad']*180/np.pi, + brec['incidence_angle_rad']*180/np.pi, + arec, + brec, + outdir, + fname='plot_by_incidence_angle_degrees.png' + ) + + +def _do_by_incidence_angle_plots( + incidence_angle_a, incidence_angle_b, arec, brec, outdir, fname +): + # Do the "by incidence angle" plots + fig = init_figure() # Initialize figure + ax = fig.add_subplot(3, 4, 1) # Initialize target range A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['target_range_m'], + title='A-Target range (m)', + xlabel='Incidence angle (rad)', + ylabel='Target range (m)', + color='black' + ) + ax = fig.add_subplot(3, 4, 2) # Initialize target area A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['target_area_m2'], + title='A-Target area ($m^2$)', + xlabel='Incidence angle (rad)', + ylabel='Target area ($m^2$)', + color='tab:blue' + ) + ax = fig.add_subplot(3, 4, 5) # Initialize radius A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['radius_m'], + title='A-Radius (m)', + xlabel='Incidence angle (rad)', + ylabel='Radius (m)', + color='tab:red' + ) + ax = fig.add_subplot(3, 4, 6) # Initialize BDRF A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['bdrf'], + title='A-BDRF', + xlabel='Incidence angle (rad)', + ylabel='BDRF', + color='tab:green' + ) + ax = fig.add_subplot(3, 4, 9) # Initialize Cross-section A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['cross_section'], + title='A-Cross-section ($m^2$)', + xlabel='Incidence angle (rad)', + ylabel='Cross-section ($m^2$)', + color='tab:orange' + ) + ax = fig.add_subplot(3, 4, 10) # Initialize received power A subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_a, arec['received_power'], + title='A-Received power', + xlabel='Incidence angle (rad)', + ylabel='Received power', + color='tab:purple' + ) + ax = fig.add_subplot(3, 4, 3) # Initialize target range B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['target_range_m'], + title='B-Target range (m)', + xlabel='Incidence angle (rad)', + ylabel='Target range (m)', + color='black' + ) + ax = fig.add_subplot(3, 4, 4) # Initialize target area B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['target_area_m2'], + title='B-Target area ($m^2$)', + xlabel='Incidence angle (rad)', + ylabel='Target area ($m^2$)', + color='tab:blue' + ) + ax = fig.add_subplot(3, 4, 7) # Initialize radius B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['radius_m'], + title='B-Radius (m)', + xlabel='Incidence angle (rad)', + ylabel='Radius (m)', + color='tab:red' + ) + ax = fig.add_subplot(3, 4, 8) # Initialize BDRF B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['bdrf'], + title='B-BDRF', + xlabel='Incidence angle (rad)', + ylabel='BDRF', + color='tab:green' + ) + ax = fig.add_subplot(3, 4, 11) # Initialize Cross-section B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['cross_section'], + title='B-Cross-section ($m^2$)', + xlabel='Incidence angle (rad)', + ylabel='Cross-section ($m^2$)', + color='tab:orange' + ) + ax = fig.add_subplot(3, 4, 12) # Initialize received power B subplot + do_y_by_x_subplot( + fig, ax, incidence_angle_b, brec['received_power'], + title='B-Received power', + xlabel='Incidence angle (rad)', + ylabel='Received power', + color='tab:purple' + ) + fig.tight_layout() + # Save figure to file and remove it from memory + fig.savefig( + os.path.join(outdir, fname) + ) + fig.clear() + plt.close(fig) + + +def do_subray_hit_subplot_hist2d( + fig, ax, x, y, title=None, xlabel=None, ylabel=None, bins="auto" +): + if title is not None: + ax.set_title(title, fontsize=15) + if bins == "auto": + bins = [len(np.unique(x)), len(np.unique(y))] + hist2d = ax.hist2d( + x, y, bins=bins, cmap='viridis', + weights=100*np.ones_like(x)/len(x), + edgecolors='black' + ) + fig.colorbar(hist2d[3]) + if xlabel is not None: + ax.set_xlabel(xlabel, fontsize=14) + if ylabel is not None: + ax.set_ylabel(ylabel, fontsize=14) + ax.tick_params(axis='both', which='both', labelsize=12) + ax.grid('both') + ax.set_axisbelow(True) + + +def do_subray_hit_subplot_hist( + fig, ax, hit, x, title=None, xlabel=None, ylabel=None, bins=7, + relative=False +): + if title is not None: + ax.set_title(title, fontsize=15) + x_hit = x[hit] + x_nohit = x[~hit] + weights = [ + 100*np.ones_like(x_hit)/len(x_hit), + 100*np.ones_like(x_nohit)/len(x_nohit) + ] if relative else None + hist = ax.hist( + [x_hit, x_nohit], bins=bins, label=['hit', 'miss'], weights=weights + ) + if xlabel is not None: + ax.set_xlabel(xlabel, fontsize=14) + if ylabel is not None: + ax.set_ylabel(ylabel, fontsize=14) + ax.tick_params(axis='both', which='both', labelsize=14) + ax.legend(loc='upper right', fontsize=12) + ax.grid('both') + ax.set_axisbelow(True) + + +def do_subray_hit_plots(arec, brec, outdir): + # Validate subray hit data + if(not validate_record('subray_hit', arec, 'a') or + not validate_record('radius_step', arec, 'a') or + not validate_record('circle_steps', arec, 'a') or + not validate_record('circle_step', arec, 'a') or + not validate_record('divergence_angle_rad', arec, 'a') or + not validate_record('subray_hit', brec, 'b') or + not validate_record('radius_step', brec, 'b') or + not validate_record('circle_steps', brec, 'b') or + not validate_record('circle_step', brec, 'b') or + not validate_record('divergence_angle_rad', brec, 'b') + ): + print('Cannot do subray hit plots') + return + + # Do the subray hit plots + fig = init_figure() # Initialize figure + # CASE A + ax = fig.add_subplot(4, 5, 1) # Initialize hit2Dhist on (radstep,circstep) + do_subray_hit_subplot_hist2d( + fig, ax, + arec['circle_step'][arec['subray_hit']], + arec['radius_step'][arec['subray_hit']], + title='Hit distribution (100%) (A)', + ) + ax = fig.add_subplot(4, 5, 2) # Initialize a hist on radius step by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['radius_step'], + ylabel='Absolute' + ) + ax = fig.add_subplot(4, 5, 3) # Initialize a hist on circle steps by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['circle_steps'], + ) + ax = fig.add_subplot(4, 5, 4) # Initialize a hist on circle step by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['circle_step'], + ) + ax = fig.add_subplot(4, 5, 5) # Initialize a hist on div. angle by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], + 1e03*arec['divergence_angle_rad']*180/np.pi, + ) + ax = fig.add_subplot(4, 5, 6) # Initialize non-hit 2D hist on (rs, cs) + do_subray_hit_subplot_hist2d( + fig, ax, + arec['circle_step'][~arec['subray_hit']], + arec['radius_step'][~arec['subray_hit']], + title='No-hit distribution (100%)', + xlabel='Circle step', + ylabel='Radius step' + ) + ax = fig.add_subplot(4, 5, 7) # Initialize a hist on radius step by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['radius_step'], + ylabel='Relative ($100\\%$)', + relative=True, + xlabel='Radius step' + ) + ax = fig.add_subplot(4, 5, 8) # Initialize a hist on circle steps by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['circle_steps'], + relative=True, + xlabel='Circle steps' + ) + ax = fig.add_subplot(4, 5, 9) # Initialize a hist on circle step by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['circle_step'], + relative=True, + xlabel='Circle step' + ) + ax = fig.add_subplot(4, 5, 10) # Initialize a hist on div. angle by hit + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], + 1e03*arec['divergence_angle_rad']*180/np.pi, + relative=True, + xlabel='Divergence angle (deg $\\times 10^{-3}$)' + ) + # CASE B + ax = fig.add_subplot(4, 5, 11) # Initialize hit2Dhist on (radstep,circstep) + do_subray_hit_subplot_hist2d( + fig, ax, + brec['circle_step'][brec['subray_hit']], + brec['radius_step'][brec['subray_hit']], + title='Hit distribution (100%) (B)', + ) + ax = fig.add_subplot(4, 5, 12) # Initialize a hist on radius step by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['radius_step'], + ylabel='Absolute' + ) + ax = fig.add_subplot(4, 5, 13) # Initialize a hist on circle steps by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['circle_steps'], + ) + ax = fig.add_subplot(4, 5, 14) # Initialize a hist on circle step by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['circle_step'], + ) + ax = fig.add_subplot(4, 5, 15) # Initialize a hist on div. angle by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], + 1e03*brec['divergence_angle_rad']*180/np.pi, + ) + ax = fig.add_subplot(4, 5, 16) # Initialize non-hit 2D hist on (rs, cs) + do_subray_hit_subplot_hist2d( + fig, ax, + brec['circle_step'][~brec['subray_hit']], + brec['radius_step'][~brec['subray_hit']], + title='No-hit distribution (100%)', + xlabel='Circle step', + ylabel='Radius step' + ) + ax = fig.add_subplot(4, 5, 17) # Initialize a hist on radius step by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['radius_step'], + ylabel='Relative ($100\\%$)', + relative=True, + xlabel='Radius step' + ) + ax = fig.add_subplot(4, 5, 18) # Initialize a hist on circle steps by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['circle_steps'], + relative=True, + xlabel='Circle steps' + ) + ax = fig.add_subplot(4, 5, 19) # Initialize a hist on circle step by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['circle_step'], + relative=True, + xlabel='Circle step' + ) + ax = fig.add_subplot(4, 5, 20) # Initialize a hist on div. angle by hit + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], + 1e03*brec['divergence_angle_rad']*180/np.pi, + relative=True, + xlabel='Divergence angle (deg $\\times 10^{-3}$)' + ) + fig.tight_layout() + # Save figure to file and remove it from memory + fig.savefig(os.path.join(outdir, 'subray_hit.png')) + fig.clear() + plt.close(fig) + + +def do_ray_subray_plots(arec, brec, outdir): + # Validate ray subray data + if( + not validate_record('ray_dir_norm', arec, 'a') or + not validate_record('subray_dir_norm', arec, 'a') or + not validate_record('ray_subray_angle_rad', arec, 'a') or + not validate_record('ray_subray_sign_check', arec, 'a') or + not validate_record('subray_tmin', arec, 'a') or + not validate_record('subray_tmax', arec, 'a') or + not validate_record('ray_dir_norm', brec, 'b') or + not validate_record('subray_dir_norm', brec, 'b') or + not validate_record('ray_subray_angle_rad', brec, 'b') or + not validate_record('ray_subray_sign_check', brec, 'b') or + not validate_record('subray_tmin', brec, 'b') or + not validate_record('subray_tmax', brec, 'b') + ): + print('Cannot do ray-subray plots') + return + # Do the ray subray plots + fig = init_figure() # Initialize figure + # CASE A + ax = fig.add_subplot(3, 4, 1) # Initialize ray norm subplot + do_incidence_angle_subplot( + fig, ax, arec['ray_dir_norm'], + xlabel='Ray direction norm' + ) + ax = fig.add_subplot(3, 4, 2) # Initialize subray norm subplot + do_incidence_angle_subplot( + fig, ax, arec['subray_dir_norm'], + xlabel='Subay direction norm' + ) + ax = fig.add_subplot(3, 4, 5) # Initialize ray-subray angle subplot + do_incidence_angle_subplot( + fig, ax, arec['ray_subray_angle_rad']*180/np.pi, + xlabel='Ray-subray angle (deg)' + ) + ax = fig.add_subplot(3, 4, 6) # Initialize ray-subray sign check subplot + do_incidence_angle_subplot( + fig, ax, arec['ray_subray_sign_check'], + xlabel='Sign equality check' + ) + ax = fig.add_subplot(3, 4, 9) # Initialize subray tmin subplot + do_incidence_angle_subplot( + fig, ax, arec['subray_tmin'], + xlabel='Subray $t_{\\mathrm{min}}$' + ) + ax = fig.add_subplot(3, 4, 10) # Initialize subray tmin subplot + do_incidence_angle_subplot( + fig, ax, arec['subray_tmax'], + xlabel='Subray $t_{\\mathrm{max}}$' + ) + # CASE B + ax = fig.add_subplot(3, 4, 3) # Initialize ray norm subplot + do_incidence_angle_subplot( + fig, ax, brec['ray_dir_norm'], + xlabel='Ray direction norm' + ) + ax = fig.add_subplot(3, 4, 4) # Initialize subray norm subplot + do_incidence_angle_subplot( + fig, ax, brec['subray_dir_norm'], + xlabel='Subay direction norm' + ) + ax = fig.add_subplot(3, 4, 7) # Initialize ray-subray angle subplot + do_incidence_angle_subplot( + fig, ax, brec['ray_subray_angle_rad']*180/np.pi, + xlabel='Ray-subray angle (deg)' + ) + ax = fig.add_subplot(3, 4, 8) # Initialize ray-subray sign check subplot + do_incidence_angle_subplot( + fig, ax, brec['ray_subray_sign_check'], + xlabel='Sign equality check' + ) + ax = fig.add_subplot(3, 4, 11) # Initialize subray tmin subplot + do_incidence_angle_subplot( + fig, ax, brec['subray_tmin'], + xlabel='Subray $t_{\\mathrm{min}}$' + ) + ax = fig.add_subplot(3, 4, 12) # Initialize subray tmin subplot + do_incidence_angle_subplot( + fig, ax, brec['subray_tmax'], + xlabel='Subray $t_{\\mathrm{max}}$' + ) + fig.tight_layout() + # Save figure to file and remove it from memory + fig.savefig( + os.path.join(outdir, 'ray_subray.png') + ) + fig.clear() + plt.close(fig) + + +def do_dir_2d_subplot( + fig, ax, x, y, hit=None, title=None, xlabel=None, ylabel=None, legend=False +): + if title is not None: + ax.set_title(title, fontsize=15) + # Plot unitary circumference + theta = np.linspace(-np.pi, np.pi) + ax.plot(np.cos(theta), np.sin(theta), color='black', lw=3, zorder=7) + if hit is not None: + ax.scatter(x[hit], y[hit], s=64, c='tab:red', zorder=6, label='hit') + ax.scatter(x[~hit], y[~hit], s=64, c='tab:blue', zorder=5, label='hit') + if legend: + ax.legend(loc='upper right').set_zorder(11) + else: + ax.scatter(x, y, s=64, c='tab:green', zorder=6) + if xlabel is not None: + ax.set_xlabel(xlabel, fontsize=14) + if ylabel is not None: + ax.set_ylabel(ylabel, fontsize=14) + ax.axis('equal') + + +def do_ray_subray_dir_plots(arec, brec, outdir): + # Validate ray and subray direction data + if( + not validate_record('subray_dir_x', arec, 'a') or + not validate_record('subray_dir_y', arec, 'a') or + not validate_record('subray_dir_z', arec, 'a') or + not validate_record('ray_dir_x', arec, 'a') or + not validate_record('ray_dir_y', arec, 'a') or + not validate_record('ray_dir_z', arec, 'a') or + not validate_record('subray_dir_x', brec, 'b') or + not validate_record('subray_dir_y', brec, 'b') or + not validate_record('subray_dir_z', brec, 'b') or + not validate_record('ray_dir_x', brec, 'b') or + not validate_record('ray_dir_y', brec, 'b') or + not validate_record('ray_dir_z', brec, 'b') + ): + print('Cannot do ray and subray direction plots') + return + # Do the ray subray direction plots + fig = init_figure() # Initialize figure + gs = plt.GridSpec(4, 1, figure=fig) + gs0 = mpl.gridspec.GridSpecFromSubplotSpec( + 1, 6, subplot_spec=gs[0], wspace=0.3 + ) + gs1 = mpl.gridspec.GridSpecFromSubplotSpec( + 1, 4, subplot_spec=gs[1] + ) + gs2 = mpl.gridspec.GridSpecFromSubplotSpec( + 1, 6, subplot_spec=gs[2], wspace=0.3 + ) + gs3 = mpl.gridspec.GridSpecFromSubplotSpec( + 1, 4, subplot_spec=gs[3] + ) + # CASE A + ax = fig.add_subplot(gs0[0, 0]) # Initialize subray dir xy subplot + subray_xynorm = np.linalg.norm(np.array([ + arec['subray_dir_x'], arec['subray_dir_y'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['subray_dir_x']/subray_xynorm, + arec['subray_dir_y']/subray_xynorm, + hit=arec['subray_hit'], + title='A subray (x, y)', + xlabel='$x$', ylabel='$y$', + legend=True + ) + ax = fig.add_subplot(gs0[0, 1]) # Initialize subray dir xz subplot + xznorm = np.linalg.norm(np.array([ + arec['subray_dir_x'], arec['subray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['subray_dir_x']/xznorm, + arec['subray_dir_z']/xznorm, + hit=arec['subray_hit'], + title='subray (x, z)', + xlabel='$x$', ylabel='$z$' + ) + ax = fig.add_subplot(gs0[0, 2]) # Initialize subray dir yz subplot + yznorm = np.linalg.norm(np.array([ + arec['subray_dir_y'], arec['subray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['subray_dir_y']/yznorm, + arec['subray_dir_z']/yznorm, + hit=arec['subray_hit'], + title='subray (y, z)', + xlabel='$y$', ylabel='$z$' + ) + ax = fig.add_subplot(gs0[0, 3]) # Initialize ray dir xy subplot + ray_xynorm = np.linalg.norm(np.array([ + arec['ray_dir_x'], arec['ray_dir_y'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['ray_dir_x']/ray_xynorm, + arec['ray_dir_y']/ray_xynorm, + title='ray (x, y)', + xlabel='$x$', ylabel='$y$', + legend=True + ) + ax = fig.add_subplot(gs0[0, 4]) # Initialize ray dir xz subplot + xznorm = np.linalg.norm(np.array([ + arec['ray_dir_x'], arec['ray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['ray_dir_x']/xznorm, + arec['ray_dir_z']/xznorm, + title='ray (x, z)', + xlabel='$x$', ylabel='$z$' + ) + ax = fig.add_subplot(gs0[0, 5]) # Initialize ray dir yz subplot + yznorm = np.linalg.norm(np.array([ + arec['ray_dir_y'], arec['ray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + arec['ray_dir_y']/yznorm, + arec['ray_dir_z']/yznorm, + title='ray (y, z)', + xlabel='$y$', ylabel='$z$' + ) + + ax = fig.add_subplot(gs1[0, 0]) # Initialize subray dir xy histogram + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], subray_xynorm, + xlabel='Subray dir norm on (x, y)', + ylabel='A Relative (100%)', + relative=True + ) + ax = fig.add_subplot(gs1[0, 1]) # Initialize subray dir xy histogram + do_subray_hit_subplot_hist( + fig, ax, arec['subray_hit'], arec['subray_dir_z'], + xlabel='Subray dir z', + relative=True + ) + ax = fig.add_subplot(gs1[0, 2]) # Initialize subray dir xy histogram + do_incidence_angle_subplot( + fig, ax, ray_xynorm, + xlabel='Ray dir norm on (x, y)', + relative=True, + label_fontsize=14 + ) + ax = fig.add_subplot(gs1[0, 3]) # Initialize subray dir xy histogram + do_incidence_angle_subplot( + fig, ax, arec['ray_dir_z'], + xlabel='Ray dir z', + relative=True, + label_fontsize=14 + ) + # CASE B + ax = fig.add_subplot(gs2[0, 0]) # Initialize subray dir xy subplot + subray_xynorm = np.linalg.norm(np.array([ + brec['subray_dir_x'], brec['subray_dir_y'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['subray_dir_x']/subray_xynorm, + brec['subray_dir_y']/subray_xynorm, + hit=brec['subray_hit'], + title='B subray (x, y)', + xlabel='$x$', ylabel='$y$', + legend=True + ) + ax = fig.add_subplot(gs2[0, 1]) # Initialize subray dir xz subplot + xznorm = np.linalg.norm(np.array([ + brec['subray_dir_x'], brec['subray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['subray_dir_x']/xznorm, + brec['subray_dir_z']/xznorm, + hit=brec['subray_hit'], + title='subray (x, z)', + xlabel='$x$', ylabel='$z$' + ) + ax = fig.add_subplot(gs2[0, 2]) # Initialize subray dir yz subplot + yznorm = np.linalg.norm(np.array([ + brec['subray_dir_y'], brec['subray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['subray_dir_y']/yznorm, + brec['subray_dir_z']/yznorm, + hit=brec['subray_hit'], + title='subray (y, z)', + xlabel='$y$', ylabel='$z$' + ) + ax = fig.add_subplot(gs2[0, 3]) # Initialize ray dir xy subplot + ray_xynorm = np.linalg.norm(np.array([ + brec['ray_dir_x'], brec['ray_dir_y'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['ray_dir_x']/ray_xynorm, + brec['ray_dir_y']/ray_xynorm, + title='ray (x, y)', + xlabel='$x$', ylabel='$y$', + legend=True + ) + ax = fig.add_subplot(gs2[0, 4]) # Initialize ray dir xz subplot + xznorm = np.linalg.norm(np.array([ + brec['ray_dir_x'], brec['ray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['ray_dir_x']/xznorm, + brec['ray_dir_z']/xznorm, + title='ray (x, z)', + xlabel='$x$', ylabel='$z$' + ) + ax = fig.add_subplot(gs2[0, 5]) # Initialize ray dir yz subplot + yznorm = np.linalg.norm(np.array([ + brec['ray_dir_y'], brec['ray_dir_z'] + ]), axis=0) + do_dir_2d_subplot( + fig, ax, + brec['ray_dir_y']/yznorm, + brec['ray_dir_z']/yznorm, + title='ray (y, z)', + xlabel='$y$', ylabel='$z$' + ) + ax = fig.add_subplot(gs3[0, 0]) # Initialize subray dir xy histogram + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], subray_xynorm, + xlabel='Subray dir norm on (x, y)', + ylabel='B Relative (100%)', + relative=True + ) + ax = fig.add_subplot(gs3[0, 1]) # Initialize subray dir xy histogram + do_subray_hit_subplot_hist( + fig, ax, brec['subray_hit'], brec['subray_dir_z'], + xlabel='Subray dir z', + relative=True + ) + ax = fig.add_subplot(gs3[0, 2]) # Initialize subray dir xy histogram + do_incidence_angle_subplot( + fig, ax, ray_xynorm, + xlabel='Ray dir norm on (x, y)', + relative=True, + label_fontsize=14 + ) + ax = fig.add_subplot(gs3[0, 3]) # Initialize subray dir xy histogram + do_incidence_angle_subplot( + fig, ax, brec['ray_dir_z'], + xlabel='Ray dir z', + relative=True, + label_fontsize=14 + ) + # Save figure to file and remove it from memory + fig.tight_layout() + fig.savefig( + os.path.join(outdir, 'ray_subray_dir.png') + ) + fig.clear() + plt.close(fig) + + +# --- M A I N --- # +# ------------------- # +if __name__ == '__main__': + # Prepare plotter + args = parse_args() + sep = ',' + # Read A records + print( + 'Reading A-records from "{path}" ...' + .format(path=args['dira_path']) + ) + start = time.perf_counter() + arec = read_records(args['dira_path'], sep=sep) + end = time.perf_counter() + print('Read A-records in {t} seconds'.format(t=end-start)) + # Read B records + print( + 'Reading B-records from "{path}" ...' + .format(path=args['dirb_path']) + ) + start = time.perf_counter() + brec = read_records(args['dirb_path'], sep=sep) + end = time.perf_counter() + print('Read B-records in {t} seconds'.format(t=end-start)) + # Plot records + print( + 'Generating plots at "{path}" ...' + .format( + path=args['dirout_path'] + ) + ) + start = time.perf_counter() + plot_records(arec, brec, args['dirout_path']) + end = time.perf_counter() + print('Generated plots in {t} seconds'.format(t=end-start)) diff --git a/scripts/debug/hda_simstep_plotter.py b/scripts/debug/hda_simstep_plotter.py old mode 100644 new mode 100755 index 92d8e0e43..9ce472fd7 --- a/scripts/debug/hda_simstep_plotter.py +++ b/scripts/debug/hda_simstep_plotter.py @@ -45,7 +45,7 @@ def parse_args(helpf=print_help): if not validate_directory(dirout_path): raise Exception( 'The directory "{d}"\n' - 'was given as the third directory of records, but it is not valid' + 'was given as the third directory for plots, but it is not valid' ) return { 'dira_path': dira_path, diff --git a/scripts/debug/hda_xyzcloud_report.py b/scripts/debug/hda_xyzcloud_report.py old mode 100644 new mode 100755 diff --git a/scripts/debug/plot_budding_metrics.py b/scripts/debug/plot_budding_metrics.py old mode 100644 new mode 100755 diff --git a/scripts/debug/plot_log_data.py b/scripts/debug/plot_log_data.py old mode 100644 new mode 100755 diff --git a/src/alg/raycast/KDGroveRaycaster.cpp b/src/alg/raycast/KDGroveRaycaster.cpp index 002fcd6fb..20342dde0 100644 --- a/src/alg/raycast/KDGroveRaycaster.cpp +++ b/src/alg/raycast/KDGroveRaycaster.cpp @@ -29,7 +29,7 @@ RaySceneIntersection * KDGroveRaycaster::search( ){ std::map out; size_t const m = grove->getNumTrees(); - RaySceneIntersection *bestRSI = nullptr, *rsi = nullptr; + RaySceneIntersection *bestRSI = nullptr, *rsi; for(size_t i = 0 ; i < m ; ++i){ rsi = grove->getTreeShared(i)->search( rayOrigin, rayDir, tmin, tmax, groundOnly diff --git a/src/alg/raycast/KDTreeRaycaster.cpp b/src/alg/raycast/KDTreeRaycaster.cpp index 70bdc048a..1c3e9367c 100644 --- a/src/alg/raycast/KDTreeRaycaster.cpp +++ b/src/alg/raycast/KDTreeRaycaster.cpp @@ -176,7 +176,6 @@ Primitive* KDTreeRaycaster::search_recursive( // ######### BEGIN If node is a leaf, perform ray-primitive intersection on all primitives in the leaf's bucket ########### if (node->splitAxis == -1) { - //logging::DEBUG("leaf node:"); for (auto prim : *node->primitives) { double const newDistance = prim->getRayIntersectionDistance( search.rayOrigin, search.rayDir @@ -198,10 +197,9 @@ Primitive* KDTreeRaycaster::search_recursive( // (in other leaves, with other primitives) that are *closer* to // the ray originWaypoint. If this was the case, the returned intersection // would be wrong. - - if( - (newDistance > 0 && newDistance < search.closestHitDistance) && - (newDistance >= tmin && newDistance <= tmax) + if( + newDistance > 0 && newDistance < search.closestHitDistance && + newDistance >= tmin && newDistance <= tmax ){ if( !search.groundOnly || @@ -257,22 +255,30 @@ Primitive* KDTreeRaycaster::search_recursive( // thit >= tmax means that the ray crosses the split plane *after it has already left the volume*. // In this case, it never enters the second half. if (thit >= tmax) { - hitPrim = search_recursive(first, tmin, tmax, search); + hitPrim = search_recursive( + first, tmin, tmax, search + ); } // thit <= tmin means that the ray crosses the split plane *before it enters the volume*. // In this case, it never enters the first half: else if (thit <= tmin) { - hitPrim = search_recursive(second, tmin, tmax, search); + hitPrim = search_recursive( + second, tmin, tmax, search + ); } // Otherwise, the ray crosses the split plane within the volume. // This means that it passes through both sides: else { - hitPrim = search_recursive(first, tmin, thit+epsilon, search); + hitPrim = search_recursive( + first, tmin, thit+epsilon, search + ); if (hitPrim == nullptr) { - hitPrim = search_recursive(second, thit-epsilon, tmax, search); + hitPrim = search_recursive( + second, thit-epsilon, tmax, search + ); } } diff --git a/src/dataanalytics/HDA_GlobalVars.cpp b/src/dataanalytics/HDA_GlobalVars.cpp new file mode 100644 index 000000000..32534f72e --- /dev/null +++ b/src/dataanalytics/HDA_GlobalVars.cpp @@ -0,0 +1,136 @@ +#if DATA_ANALYTICS >= 2 +#include + +namespace helios { namespace analytics{ + + +// *** GLOBAL OBJECT *** // +// *********************** // +HDA_GlobalVars HDA_GV; + +// *** WRITE METHODS *** // +// *********************** // +HDA_GlobalVars & HDA_GlobalVars::incrementGeneratedSubraysCount(){ + std::unique_lock lock(generatedSubraysCount_mutex); + ++generatedSubraysCount; + return *this; +} + +HDA_GlobalVars & + HDA_GlobalVars::incrementGeneratedRaysBeforeEarlyAbortCount(){ + std::unique_lock lock( + generatedRaysBeforeEarlyAbortCount_mutex + ); + ++generatedRaysBeforeEarlyAbortCount; + return *this; +} + +HDA_GlobalVars & + HDA_GlobalVars::incrementGeneratedRaysAfterEarlyAbortCount(){ + std::unique_lock lock( + generatedRaysAfterEarlyAbortCount_mutex + ); + ++generatedRaysAfterEarlyAbortCount; + return *this; +} + +HDA_GlobalVars & + HDA_GlobalVars::incrementIntensityComputationsCount() { + std::unique_lock lock(intensityComputationsCount_mutex); + ++intensityComputationsCount; + return *this; +} + +HDA_GlobalVars & HDA_GlobalVars::incrementIntersectiveSubraysCount(){ + std::unique_lock lock(intersectiveSubraysCount_mutex); + ++intersectiveSubraysCount; + return *this; +} + +HDA_GlobalVars & HDA_GlobalVars::incrementNonIntersectiveSubraysCount(){ + std::unique_lock lock(nonIntersectiveSubraysCount_mutex); + ++nonIntersectiveSubraysCount; + return *this; +} + +HDA_GlobalVars & +HDA_GlobalVars::incrementNonIntersectiveSubraysDueToNullTimeCount(){ + std::unique_lock lock( + nonIntersectiveSubraysDueToNullTimeCount_mutex + ); + ++nonIntersectiveSubraysDueToNullTimeCount; + return *this; + +} + +HDA_GlobalVars & HDA_GlobalVars::incrementSubrayIntersectionCount() { + std::unique_lock lock(subrayIntersectionCount_mutex); + ++subrayIntersectionCount; + return *this; +} + +HDA_GlobalVars & HDA_GlobalVars::incrementSubrayNonIntersectionCount() { + std::unique_lock lock(subrayNonIntersectionCount_mutex); + ++subrayNonIntersectionCount; + return *this; +} + + +// *** READ METHODS *** // +// ********************** // +std::size_t HDA_GlobalVars::getGeneratedSubraysCount(){ + std::unique_lock lock(generatedSubraysCount_mutex); + return generatedSubraysCount; +} + +std::size_t HDA_GlobalVars::getGeneratedRaysBeforeEarlyAbortCount(){ + std::unique_lock lock( + generatedRaysBeforeEarlyAbortCount_mutex + ); + return generatedRaysBeforeEarlyAbortCount; +} + +std::size_t HDA_GlobalVars::getGeneratedRaysAfterEarlyAbortCount(){ + std::unique_lock lock( + generatedRaysAfterEarlyAbortCount_mutex + ); + return generatedRaysAfterEarlyAbortCount; +} + +std::size_t HDA_GlobalVars::getIntersectiveSubraysCount() { + std::unique_lock lock(intersectiveSubraysCount_mutex); + return intersectiveSubraysCount; +} + +std::size_t HDA_GlobalVars::getNonIntersectiveSubraysCount() { + std::unique_lock lock(nonIntersectiveSubraysCount_mutex); + return nonIntersectiveSubraysCount; +} + +std::size_t HDA_GlobalVars::getNonIntersectiveSubraysDueToNullTimeCount(){ + std::unique_lock lock( + nonIntersectiveSubraysDueToNullTimeCount_mutex + ); + return nonIntersectiveSubraysDueToNullTimeCount; +} + +std::size_t HDA_GlobalVars::getSubrayIntersectionCount() { + std::unique_lock lock(subrayIntersectionCount_mutex); + return subrayIntersectionCount; +} + +std::size_t HDA_GlobalVars::getSubrayNonIntersectionCount() { + std::unique_lock lock(subrayNonIntersectionCount_mutex); + return subrayNonIntersectionCount; +} + +std::size_t HDA_GlobalVars::getIntensityComputationsCount(){ + std::unique_lock lock(intensityComputationsCount_mutex); + return intensityComputationsCount; +} + + + +}} + +#endif diff --git a/src/dataanalytics/HDA_GlobalVars.h b/src/dataanalytics/HDA_GlobalVars.h new file mode 100644 index 000000000..20912982e --- /dev/null +++ b/src/dataanalytics/HDA_GlobalVars.h @@ -0,0 +1,95 @@ +#if DATA_ANALYTICS >=2 +#pragma once +#include +#include + +namespace helios { namespace analytics{ + +// *** GLOBAL OBJECT *** // +// *********************** // +extern class HDA_GlobalVars HDA_GV; + +/** + * @author Alberto M. Esmoris Pena + * @version 1.0 + * + * @brief Class to handle global variables together with methods for safe + * readings and updates. + * + * The HDA_GlobalVars class provides a series of variables and corresponding + * methods for safe operations. For example, variables that must be + * accessed through parallel threads are handled with a mutex such that, as + * long as the proper methods are used, read and write operations are + * thread safe. + */ +class HDA_GlobalVars{ +public: + // *** ATTRIBUTES *** // + // ******************** // + std::size_t generatedRaysBeforeEarlyAbortCount; + std::size_t generatedRaysAfterEarlyAbortCount; + std::size_t generatedSubraysCount; + std::size_t intersectiveSubraysCount; + std::size_t nonIntersectiveSubraysCount; + std::size_t nonIntersectiveSubraysDueToNullTimeCount; + std::size_t subrayIntersectionCount; + std::size_t subrayNonIntersectionCount; + std::size_t intensityComputationsCount; + +protected: + // *** CONCURRENCY HANDLING ATTRIBUTES *** // + // ***************************************** // + std::mutex generatedRaysBeforeEarlyAbortCount_mutex; + std::mutex generatedRaysAfterEarlyAbortCount_mutex; + std::mutex generatedSubraysCount_mutex; + std::mutex intersectiveSubraysCount_mutex; + std::mutex nonIntersectiveSubraysCount_mutex; + std::mutex nonIntersectiveSubraysDueToNullTimeCount_mutex; + std::mutex subrayIntersectionCount_mutex; + std::mutex subrayNonIntersectionCount_mutex; + std::mutex intensityComputationsCount_mutex; + +public: + // *** CONSTRUCTION / DESTRUCTION *** // + // ************************************ // + HDA_GlobalVars() : + generatedRaysBeforeEarlyAbortCount(0), + generatedRaysAfterEarlyAbortCount(0), + generatedSubraysCount(0), + intersectiveSubraysCount(0), + nonIntersectiveSubraysCount(0), + nonIntersectiveSubraysDueToNullTimeCount(0), + subrayIntersectionCount(0), + subrayNonIntersectionCount(0), + intensityComputationsCount(0) + {} + virtual ~HDA_GlobalVars() = default; + + // *** WRITE METHODS *** // + // *********************** // + HDA_GlobalVars & incrementGeneratedRaysBeforeEarlyAbortCount(); + HDA_GlobalVars & incrementGeneratedRaysAfterEarlyAbortCount(); + HDA_GlobalVars & incrementGeneratedSubraysCount(); + HDA_GlobalVars & incrementIntersectiveSubraysCount(); + HDA_GlobalVars & incrementNonIntersectiveSubraysCount(); + HDA_GlobalVars & incrementNonIntersectiveSubraysDueToNullTimeCount(); + HDA_GlobalVars & incrementSubrayIntersectionCount(); + HDA_GlobalVars & incrementSubrayNonIntersectionCount(); + HDA_GlobalVars & incrementIntensityComputationsCount(); + + // *** READ METHODS *** // + // ********************** // + std::size_t getGeneratedSubraysCount(); + std::size_t getGeneratedRaysBeforeEarlyAbortCount(); + std::size_t getGeneratedRaysAfterEarlyAbortCount(); + std::size_t getIntersectiveSubraysCount(); + std::size_t getNonIntersectiveSubraysCount(); + std::size_t getNonIntersectiveSubraysDueToNullTimeCount(); + std::size_t getSubrayIntersectionCount(); + std::size_t getSubrayNonIntersectionCount(); + std::size_t getIntensityComputationsCount(); +}; + +}} + +#endif diff --git a/src/dataanalytics/HDA_GlobalVarsReporter.h b/src/dataanalytics/HDA_GlobalVarsReporter.h new file mode 100644 index 000000000..d25a72230 --- /dev/null +++ b/src/dataanalytics/HDA_GlobalVarsReporter.h @@ -0,0 +1,64 @@ +#if DATA_ANALYTICS >= 2 +#pragma once +#include +#include + +#include +#include + +namespace helios { namespace analytics{ + +/** + * @author Alberto M. Esmoris Pena + * @version 1.0 + * + * @brief Class to handle reports related to global variables. + * + * @see helios::analytics::HDA_GlobalVars + */ +class HDA_GlobalVarsReporter{ +protected: + // *** ATTRIBUTES *** // + // ******************** // + HDA_GlobalVars & gv; + +public: + // *** CONSTRUCTION / DESTRUCTION *** // + // ************************************ // + explicit HDA_GlobalVarsReporter(HDA_GlobalVars &gv) : gv(gv) {} + virtual ~HDA_GlobalVarsReporter() = default; + + // *** PRINT *** // + // ***************** // + void print(){ + // Initialize string stream to build the print + std::stringstream ss; + + // Print global variables + ss << "HDA GLOBAL VARS REPORT:\n\n"; + ss << "Generated rays before early abort: " + << gv.getGeneratedRaysBeforeEarlyAbortCount() << "\n"; + ss << "Generated rays after early abort: " + << gv.getGeneratedRaysAfterEarlyAbortCount() << "\n"; + ss << "Generated subrays: " << gv.getGeneratedSubraysCount() << "\n"; + ss << "Intersective subrays: " + << gv.getIntersectiveSubraysCount() << "\n"; + ss << "Non-intersective subrays: " + << gv.getNonIntersectiveSubraysCount() << "\n"; + ss << "\tNon-intersective subrays due to null time: " + << gv.getNonIntersectiveSubraysDueToNullTimeCount() << "\n"; + ss << "Subray intersections: " + << gv.getSubrayIntersectionCount() << "\n"; + ss << "Subray non-intersections: " + << gv.getSubrayNonIntersectionCount() << "\n"; + ss << "Number of computed intensities: " + << gv.getIntensityComputationsCount() << "\n"; + // Print through info logging level system + std::string text = ss.str(); + logging::INFO(ss.str()); + } +}; + +}} + +#endif diff --git a/src/dataanalytics/HDA_OfstreamWrapper.h b/src/dataanalytics/HDA_OfstreamWrapper.h new file mode 100644 index 000000000..88df13ac2 --- /dev/null +++ b/src/dataanalytics/HDA_OfstreamWrapper.h @@ -0,0 +1,96 @@ +#ifdef DATA_ANALYTICS + +#include +#include +#include + + +namespace helios { namespace analytics{ + +class HDA_OfstreamWrapper{ +protected: + // *** ATTRIBUTES *** // + // ******************** // + /** + * @brief The output file stream used to write to a file. + */ + std::ofstream ofs; + /** + * @brief The separator between recorded elements. + */ + std::string sep; + /** + * @brief Boolean flag to control whether the HDA_OfstreamWrapper has + * written something or not. Once something has been written, the flag + * will be set to true. By default, it is initialized to false. + */ + bool hasWrittenSomething; + +public: + // *** CONSTRUCTION / DESTRUCTION *** // + // ************************************ // + /** + * @brief Build a HDA_OfstreamWrapper to handle writing records to files. + * @param outpath The path to the output file. + * @see HDA_OfstreamWrapper::sep + */ + HDA_OfstreamWrapper( + std::string const &outpath, std::string const &sep="," + ) : + ofs(outpath, std::ios_base::out), + sep(sep), + hasWrittenSomething(false) + {} + virtual ~HDA_OfstreamWrapper() = default; + + // *** OUTPUT FILE STREAM METHODS *** // + // ************************************ // + /** + * @brief Wraps the is_open method of the output file stream class. + */ + inline bool is_open() {return ofs.is_open();} + /** + * @brief Wraps the close method of the output file stream class. + */ + inline void close() {return ofs.close();} + + // *** OUTPUT FILE STREAM OPERATORS *** // + // ************************************** // + /** + * @brief Default stream input operator, i.e., <<. + * @tparam T The type of element to be written. + * @param elem The element to be written. + * @return A reference to the HDA_OfstreamWrapper itself for fluent + * programming purposes. + */ + template + HDA_OfstreamWrapper & operator <<(T const &elem){ + if(hasWrittenSomething){ + ofs << sep << elem; + } + else{ + ofs << elem; + hasWrittenSomething = true; + } + return *this; + } + + /** + * @brief Stream input operator to handle vectors of doubles as element. + * @param elem The vector of doubles as element. + * @return A reference to the HDA_OfstreamWrapper itself for fluent + * programming purposes. + */ + HDA_OfstreamWrapper & operator <<(std::vector const &elem){ + ofs << elem[0]; + for(size_t i = 1 ; i < elem.size(); ++i){ + ofs << sep << elem[i]; + } + ofs << "\n"; + return *this; + } + +}; + +}} +#endif \ No newline at end of file diff --git a/src/dataanalytics/HDA_PulseRecorder.cpp b/src/dataanalytics/HDA_PulseRecorder.cpp new file mode 100644 index 000000000..0b27f94e5 --- /dev/null +++ b/src/dataanalytics/HDA_PulseRecorder.cpp @@ -0,0 +1,74 @@ +#ifdef DATA_ANALYTICS +#include + +using namespace helios::analytics; + + +// *** RECORDER METHODS *** // +// ************************** // +bool HDA_PulseRecorder::isAnyBufferOpen(){ + bool anyOpen = false; + anyOpen |= intensityCalc->isOpen(); + anyOpen |= subraySim->isOpen(); + return anyOpen; +} + +void HDA_PulseRecorder::openBuffers(){ + // Open subray related buffers + size_t const maxSize = 256; + std::string const sep = ","; + intensityCalc = std::make_shared>>( + craftOutputPath("intensity_calc.csv"), + maxSize, + sep, + true // vectorial flag + ); + subraySim = std::make_shared>>( + craftOutputPath("subray_sim.csv"), + maxSize, + sep, + true // vectorial flag + ); +} + +void HDA_PulseRecorder::closeBuffers(){ + // Close subray buffers + std::unique_lock lock(intensityCalcMutex); + intensityCalc->close(); + subraySim->close(); +} + + +// *** RECORD METHODS *** // +// ************************ // +void HDA_PulseRecorder::recordIntensityCalculation( + std::vector const &record +){ + std::unique_lock lock(intensityCalcMutex); + intensityCalc->push(record); +} +void HDA_PulseRecorder::recordIntensityCalculation( + std::vector> const &records +){ + std::unique_lock lock(intensityCalcMutex); + for(std::vector const & record : records) { + intensityCalc->push(record); + } +} + +void HDA_PulseRecorder::recordSubraySimulation( + std::vector const &record +){ + std::unique_lock lock(subraySimMutex); + subraySim->push(record); +} + +void HDA_PulseRecorder::recordSubraySimulation( + std::vector> const &records +){ + std::unique_lock lock(subraySimMutex); + for(std::vector const & record : records){ + subraySim->push(record); + } +} +#endif diff --git a/src/dataanalytics/HDA_PulseRecorder.h b/src/dataanalytics/HDA_PulseRecorder.h new file mode 100644 index 000000000..dce8fb55e --- /dev/null +++ b/src/dataanalytics/HDA_PulseRecorder.h @@ -0,0 +1,170 @@ +#ifdef DATA_ANALYTICS +#pragma once + +#include +#include + +#include +#include +#include + +namespace helios { namespace analytics{ + +/** + * @author Alberto M. Esmoris Pena + * @version 1.0 + * @brief The HeliosDataAnalytics recorder for pulses (more concretely, pulse + * tasks / runnables). It is a class which records relevant data from the + * many pulse computations. + * @see FullWaveformPulseRunnable + */ +class HDA_PulseRecorder : public HDA_Recorder{ +protected: + // *** ATTRIBUTES *** // + // ******************** // + /** + * @brief The vectors which components are variables involved on a + * particular intensity calculation for a given subray. + * + * [0, 1, 2] -> \f$(x, y, z)\f$ + * + * [3] -> Incidence angle in radians. + * + * [4] -> The target range in meters, i.e., the distance between the beam's + * origin and the intersection point. + * + * [5] -> The target area in squared meters. + * + * [6] -> The radius in meters, i.e., the distance between the beam's + * center line and the intersection point. + * + * [7] -> The bidirectional reflectance function (BDRF). + * + * [8] -> The cross-section in squared meters. + * + * [9] -> The calculated received power, i.e., intensity. + * + * [10] -> 1 if the point was captured, 0 otherwise. + */ + std::shared_ptr>> intensityCalc; + /** + * @brief The vectors which components are variables involved on the + * subray simulation. + * + * [0] -> Subray hit (0 does not hit, 1 hit) + * + * [1] -> Radius step + * + * [2] -> Circle steps + * + * [3] -> Circle step + * + * [4] -> Divergence angle (in rad) + * + * [5] -> Ray direction norm + * + * [6] -> Subray direction norm + * + * [7] -> Angle between ray and subray (in rad) + * + * [8] -> Ray-subray sign check (1 if sign match, 0 otherwise) + * + * [9] -> Min time for subray intersection + * + * [10] -> Max time for subray intersection + * + * [11] -> Subray direction (x component) + * + * [12] -> Subray direction (y component) + * + * [13] -> Subray direction (z component) + * + * [14] -> Ray direction (x component) + * + * [15] -> Ray direction (y component) + * + * [16] -> Ray direction (z component) + */ + std::shared_ptr>> subraySim; + + /** + * @brief The mutex to handle concurrent writes to the buffers related to + * intensity calculation. + */ + std::mutex intensityCalcMutex; + /** + * @brief The mutex to handle concurrent writes to the buffers related to + * subray simulation. + */ + std::mutex subraySimMutex; + +public: + // *** CONSTRUCTION / DESTRUCTION *** // + // ************************************ // + /** + * @brief Build a HDA_PulseRecorder so pulse computations are written to + * data files in the directory specified through given path. + * @param path Path of the directory where simulation records will be + * written. + */ + HDA_PulseRecorder(std::string const &path) : HDA_Recorder(path){ + openBuffers(); + } + + virtual ~HDA_PulseRecorder() { + if(isAnyBufferOpen()) closeBuffers(); + } + + // *** RECORDER METHODS *** // + // ************************** // + /** + * @brief Check whether there are opened buffers or not + * @return True if there is at least one opened buffer, false if there is + * not even a single opened buffer + */ + bool isAnyBufferOpen(); + /** + * @brief Open all the record buffers so the HDA_SimStepRecord can record + */ + void openBuffers(); + /** + * @brief Close all the record buffers. Once it is done, the + * HDA_SimStepRecorder will not be able to properly handle any new + * record. + */ + void closeBuffers(); + + + // *** RECORD METHODS *** // + // ************************ // + /** + * @brief Handle all the records for the current simulation step. + */ + virtual void recordIntensityCalculation(std::vector const &record); + /** + * @brief Like + * HDA_PulseRecorder::recordIntensityCalculation(std::vector) + * but receiving many records at once. + * @see HDA_PulseRecorder::recordIntensityCalculation(std::vector) + */ + virtual void recordIntensityCalculation( + std::vector> const &records + ); + /** + * @brief Handle all the records for the current subray simulation. + */ + virtual void recordSubraySimulation(std::vector const &record); + /** + * @brief Like + * HDA_PulseRecorder::recordSubraySimulation(std::vector) + * but receiving many records at once. + * @see HDA_PulseRecorder::recordSubraySimulation(std::vector) + */ + virtual void recordSubraySimulation( + std::vector> const &records + ); + +}; + +}} +#endif \ No newline at end of file diff --git a/src/dataanalytics/HDA_RecordBuffer.h b/src/dataanalytics/HDA_RecordBuffer.h index 10283743d..6ee1de8ff 100644 --- a/src/dataanalytics/HDA_RecordBuffer.h +++ b/src/dataanalytics/HDA_RecordBuffer.h @@ -1,6 +1,8 @@ #ifdef DATA_ANALYTICS #pragma once +#include + #include #include #include @@ -34,18 +36,9 @@ class HDA_RecordBuffer { */ std::string outpath; /** - * @brief The output stream to write the contents of the buffer - */ - std::ofstream ofs; - /** - * @brief The separator between recorded elements - */ - std::string sep; - /** - * @brief The function to write the content of the buffer through the - * output stream + * @brief The wrapped output stream to write the contents of the buffer */ - std::function write; + HDA_OfstreamWrapper ofsw; public: // *** CONSTRUCTION / DESTRUCTION *** // @@ -59,15 +52,13 @@ class HDA_RecordBuffer { HDA_RecordBuffer( std::string const &outpath, size_t const maxSize=256, - std::string const &sep="," + std::string const &sep=",", + bool const vectorial=false ) : maxSize(maxSize), outpath(outpath), - ofs(outpath, std::ios_base::out), - sep(sep) - { - write = [&](void) -> void {this->firstWrite();}; - } + ofsw(outpath, sep) + {} virtual ~HDA_RecordBuffer(){ if(isOpen()) close(); @@ -80,32 +71,17 @@ class HDA_RecordBuffer { * not * @return True if the output stream is opened, false otherwise */ - inline bool isOpen() {return ofs.is_open();} - /** - * @brief Write the contents of the buffer through the output stream for - * the first time - */ - inline void firstWrite() { - size_t numElems = buff.size(); - ofs << buff[0]; - for(size_t i = 1 ; i < numElems ; ++i){ - ofs << sep << buff[i]; - } - this->write = [&] (void) -> void {this->nextWrite();}; - } + inline bool isOpen() {return ofsw.is_open();} /** - * @brief Write the contents of the buffer through the output stream after - * the first time + * @brief Method to write the elements of the buffer to a file. */ - inline void nextWrite() { - for(T const & elem : buff) ofs << sep << elem; - } + inline void write() {for(T const & elem : buff) ofsw << elem;} /** * @brief Write all the contents of the buffer and then make it empty */ inline void flush(){ if(buff.size() > 0){ - this->write(); + write(); buff.clear(); } } @@ -115,7 +91,7 @@ class HDA_RecordBuffer { */ inline void close(){ flush(); - if(isOpen()) ofs.close(); + if(isOpen()) ofsw.close(); } /** * @brief Insert given element in the buffer. If the buffer is full, it diff --git a/src/dataanalytics/HDA_Recorder.cpp b/src/dataanalytics/HDA_Recorder.cpp new file mode 100644 index 000000000..3618c668c --- /dev/null +++ b/src/dataanalytics/HDA_Recorder.cpp @@ -0,0 +1,36 @@ +#ifdef DATA_ANALYTICS + +#include +#include + +#include + +#include + + +using namespace helios::analytics; + +// *** RECORDER METHODS *** // +// ************************** // +void HDA_Recorder::validateOutDir(){ + // Check directory exists + if(!boost::filesystem::exists(outdir)){ + if(!boost::filesystem::create_directory(outdir)){ + std::stringstream ss; + ss << "HDA_SimStepRecorder::validateOutDir thrown an exception " + << "because the output directory does not exist and cannot be" + << "created.\noutdir: \"" << outdir << "\""; + throw HeliosException(ss.str()); + } + } +} + +std::string HDA_Recorder::craftOutputPath(std::string const &fname){ + std::stringstream ss; + ss << outdir + << boost::filesystem::path::preferred_separator + << fname; + return ss.str(); +} + +#endif diff --git a/src/dataanalytics/HDA_Recorder.h b/src/dataanalytics/HDA_Recorder.h new file mode 100644 index 000000000..dc51e69fa --- /dev/null +++ b/src/dataanalytics/HDA_Recorder.h @@ -0,0 +1,61 @@ +#ifdef DATA_ANALYTICS +#pragma once + +#include + +namespace helios { namespace analytics{ + +/** + * @author Alberto M. Esmoris Pena + * @version 1.0 + * @brief The HeliosDataAnalytics abstract recorder. It is, an abstract class + * that handle only the common baseline logic for any recorder. + */ +class HDA_Recorder{ +protected: + // *** ATTRIBUTES *** // + // ******************** // + /** + * @brief Path to the output directory where the different outputs will be + * stored + */ + std::string outdir; + + // *** CONSTRUCTION / DESTRUCTION *** // + // ************************************ // + /** + * @brief Build a HDA_Recorder to write data files inside given output + * directory. + * @param path Path to the output directory. + */ + HDA_Recorder(std::string const &path) : outdir(path) { + validateOutDir(); + } + + virtual ~HDA_Recorder(){} + + // *** RECORDER METHODS *** // + // ************************** // + /** + * @brief Check whether the output directory is a valid one or not. If not, + * an exception will be thrown. + * + * This method can be overriden by derived classes when necessary, but it + * is not a must. + */ + virtual void validateOutDir(); + + /** + * @brief Craft the full output path considering the output directory and + * the given file name + * @param fname The given filename + * @return Full output path considering the output directory and the given + * file name + */ + virtual std::string craftOutputPath(std::string const &fname); + +}; + +}} + +#endif \ No newline at end of file diff --git a/src/dataanalytics/HDA_SimStepRecorder.cpp b/src/dataanalytics/HDA_SimStepRecorder.cpp index 9f164ea07..cd8a174e4 100644 --- a/src/dataanalytics/HDA_SimStepRecorder.cpp +++ b/src/dataanalytics/HDA_SimStepRecorder.cpp @@ -34,21 +34,9 @@ void HDA_SimStepRecorder::delayedRecord(){ recordStochastic(); } + // *** RECORDER METHODS *** // // ************************** // -void HDA_SimStepRecorder::validateOutDir(){ - // Check directory exists - if(!boost::filesystem::exists(outdir)){ - if(!boost::filesystem::create_directory(outdir)){ - std::stringstream ss; - ss << "HDA_SimStepRecorder::validateOutDir thrown an exception " - << "because the output directory does not exist and cannot be" - << "created.\noutdir: \"" << outdir << "\""; - throw HeliosException(ss.str()); - } - } -} - bool HDA_SimStepRecorder::isAnyBufferOpen(){ bool anyOpen = false; anyOpen |= platformPositionX->isOpen(); @@ -326,9 +314,16 @@ void HDA_SimStepRecorder::recordScanner(){ scannerPositionZ->push(pos.z); // Record scanner angles double roll, pitch, yaw; - s.getHeadRelativeEmitterAttitude().getAngles( - &RotationOrder::XYZ, roll, pitch, yaw - ); + try { + s.getHeadRelativeEmitterAttitude().getAngles( + &RotationOrder::XYZ, roll, pitch, yaw + ); + } + catch(HeliosException &hex){ + roll = std::numeric_limits::quiet_NaN(); + pitch = std::numeric_limits::quiet_NaN(); + yaw = std::numeric_limits::quiet_NaN(); + } scannerRoll->push(roll); scannerPitch->push(pitch); scannerYaw->push(yaw); @@ -398,15 +393,18 @@ void HDA_SimStepRecorder::recordStochastic(){ // Obtain parallel randomness generator RandomnessGenerator *rg1Par = nullptr; + size_t nthreads = 0; WarehouseScanningPulseProcess * wspp = dynamic_cast(spp); if(wspp != nullptr){ rg1Par = wspp->pool.randGens; + nthreads = wspp->pool.getPoolSize(); } else{ BuddingScanningPulseProcess *bspp = dynamic_cast(spp); rg1Par = bspp->pool.randGens; + nthreads = wspp->pool.getPoolSize(); } // Prepare fake pulse @@ -434,23 +432,16 @@ void HDA_SimStepRecorder::recordStochastic(){ } // Record parallel measurement error - for(size_t i = 0 ; i < N_SAMPLES ; ++i) { - err = ERR_BASE; - fwpr.applyMeasurementError( - *rg1Par, err, fakePulse.getOriginRef(), fakePulse.getOriginRef() - ); - measErrPar->push(err-ERR_BASE); + if(nthreads > 0) { // There is no parallelism for empty thread pools + for (size_t i = 0; i < N_SAMPLES; ++i) { + err = ERR_BASE; + fwpr.applyMeasurementError( + *rg1Par, err, fakePulse.getOriginRef(), fakePulse.getOriginRef() + ); + measErrPar->push(err - ERR_BASE); + } } - - } -std::string HDA_SimStepRecorder::craftOutputPath(std::string const &fname){ - std::stringstream ss; - ss << outdir - << boost::filesystem::path::preferred_separator - << fname; - return ss.str(); -} #endif \ No newline at end of file diff --git a/src/dataanalytics/HDA_SimStepRecorder.h b/src/dataanalytics/HDA_SimStepRecorder.h index d2e474d21..fcb15a0a3 100644 --- a/src/dataanalytics/HDA_SimStepRecorder.h +++ b/src/dataanalytics/HDA_SimStepRecorder.h @@ -1,6 +1,7 @@ #ifdef DATA_ANALYTICS #pragma once +#include #include #include @@ -16,7 +17,7 @@ namespace helios { namespace analytics{ * which records the data representing the state of simulation components at * a certain simulation step. */ -class HDA_SimStepRecorder{ +class HDA_SimStepRecorder : public HDA_Recorder { protected: // *** ATTRIBUTES *** // // ******************** // @@ -24,11 +25,6 @@ class HDA_SimStepRecorder{ * @brief The SurveyPlayback which simulation steps must be recorded */ SurveyPlayback *sp; - /** - * @brief Path to the output directory where the different outputs will be - * stored - */ - std::string outdir; /** * @brief The record buffer for the x coordinate of the platform position */ @@ -198,23 +194,20 @@ class HDA_SimStepRecorder{ SurveyPlayback *sp, std::string const &path ) : - sp(sp), - outdir(path) + HDA_Recorder(path), + sp(sp) { - validateOutDir(); openBuffers(); } virtual ~HDA_SimStepRecorder() { - if(isAnyBufferOpen()){ - closeBuffers(); - } + if(isAnyBufferOpen()) closeBuffers(); } // *** MAIN RECORD METHOD *** // // **************************** // /** - * @brief Handle all the records for the current simulation step + * @brief Handle all the records for the current simulation step. */ virtual void record(); /** @@ -226,27 +219,22 @@ class HDA_SimStepRecorder{ // *** RECORDER METHODS *** // // ************************** // - /** - * @brief Check whether the output directory is a valid one or not. If not, - * an exception will be thrown - */ - virtual void validateOutDir(); /** * @brief Check whether there are opened buffers or not * @return True if there is at least one opened buffer, false if there is * not even a single opened buffer */ - virtual bool isAnyBufferOpen(); + bool isAnyBufferOpen(); /** * @brief Open all the record buffers so the HDA_SimStepRecord can record */ - virtual void openBuffers(); + void openBuffers(); /** * @brief Close all the record buffers. Once it is done, the * HDA_SimStepRecorder will not be able to properly handle any new - * record + * record. */ - virtual void closeBuffers(); + void closeBuffers(); protected: // *** CONCRETE RECORD METHODS *** // @@ -286,14 +274,6 @@ class HDA_SimStepRecorder{ * current simulation step */ virtual void recordStochastic(); - /** - * @brief Craft the full output path considering the output directory and - * the given file name - * @param fname The given filename - * @return Full output path considering the output directory and the given - * file name - */ - virtual std::string craftOutputPath(std::string const &fname); }; diff --git a/src/main/Main.cpp b/src/main/Main.cpp index 7644a18f8..9e73d4289 100644 --- a/src/main/Main.cpp +++ b/src/main/Main.cpp @@ -9,6 +9,11 @@ #include #endif +#ifdef DATA_ANALYTICS +#include +#include +#endif + #include #include @@ -161,6 +166,13 @@ int main(int argc, char** argv) { ); } +#if DATA_ANALYTICS >= 2 + helios::analytics::HDA_GlobalVarsReporter reporter( + helios::analytics::HDA_GV + ); + reporter.print(); +#endif + // Return successful exit status code (0) return EXIT_SUCCESS; } diff --git a/src/main/helios_version.cpp b/src/main/helios_version.cpp index 1d7f655bb..bfc9db9ca 100644 --- a/src/main/helios_version.cpp +++ b/src/main/helios_version.cpp @@ -4,7 +4,7 @@ const char * HELIOS_VERSION = "1.2.0"; -const char * HELIOS_GIT_HASH = "03b19a3a"; +const char * HELIOS_GIT_HASH = "5cba90e3"; const char * getHeliosVersion(){ return HELIOS_VERSION; diff --git a/src/scanner/BuddingScanningPulseProcess.cpp b/src/scanner/BuddingScanningPulseProcess.cpp index 410f0b756..4df3fc589 100644 --- a/src/scanner/BuddingScanningPulseProcess.cpp +++ b/src/scanner/BuddingScanningPulseProcess.cpp @@ -10,6 +10,9 @@ BuddingScanningPulseProcess::BuddingScanningPulseProcess( RandomnessGenerator &randGen1, RandomnessGenerator &randGen2, UniformNoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) : ScanningPulseProcess(scanner), dropper(dropper), @@ -17,6 +20,9 @@ BuddingScanningPulseProcess::BuddingScanningPulseProcess( randGen1(randGen1), randGen2(randGen2), intersectionHandlingNoiseSource(intersectionHandlingNoiseSource) +#if DATA_ANALYTICS >= 2 + ,pulseRecorder(pulseRecorder) +#endif { if(pool.getPoolSize() > 0){ if(pool.isDynamic()){ // Dynamic chunk schedule @@ -51,6 +57,9 @@ void BuddingScanningPulseProcess::onLegComplete(){ randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); #ifdef BUDDING_METRICS ofsBudding.flush(); @@ -79,6 +88,9 @@ void BuddingScanningPulseProcess::handlePulseComputationSequential( randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } void BuddingScanningPulseProcess::handlePulseComputationParallelDynamic( @@ -132,6 +144,9 @@ void BuddingScanningPulseProcess::handlePulseComputationParallelDynamic( randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } } @@ -177,6 +192,9 @@ void BuddingScanningPulseProcess::handlePulseComputationParallelStatic( randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } } diff --git a/src/scanner/BuddingScanningPulseProcess.h b/src/scanner/BuddingScanningPulseProcess.h index f0e5ea279..f4f20520b 100644 --- a/src/scanner/BuddingScanningPulseProcess.h +++ b/src/scanner/BuddingScanningPulseProcess.h @@ -18,7 +18,15 @@ class BuddingScanningPulseProcess : public ScanningPulseProcess { #ifdef DATA_ANALYTICS public: #else -protected: + protected: +#endif +#if DATA_ANALYTICS >= 2 + /** + * @brief The helios::analytics::PulseRecorder to be used to handle the + * records representing the computed pulse tasks. + * @see helios::analytics::PulseRecorder + */ + std::shared_ptr pulseRecorder; #endif // *** ATTRIBUTES *** // // ******************** // @@ -109,6 +117,9 @@ class BuddingScanningPulseProcess : public ScanningPulseProcess { RandomnessGenerator &randGen1, RandomnessGenerator &randGen2, UniformNoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ); virtual ~BuddingScanningPulseProcess() = default; diff --git a/src/scanner/MultiScanner.cpp b/src/scanner/MultiScanner.cpp index 85c455774..b1cb5ed4f 100644 --- a/src/scanner/MultiScanner.cpp +++ b/src/scanner/MultiScanner.cpp @@ -158,7 +158,6 @@ Rotation MultiScanner::calcAbsoluteBeamAttitude(size_t const idx){ } void MultiScanner::computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -166,19 +165,27 @@ void MultiScanner::computeSubrays( NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ){ scanDevs[idx].computeSubrays( handleSubray, - tMinMax, intersectionHandlingNoiseSource, reflections, intersects +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } @@ -212,6 +219,9 @@ double MultiScanner::calcIntensity( double const targetArea, double const radius, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const{ return scanDevs[idx].calcIntensity( incidenceAngle, @@ -219,6 +229,9 @@ double MultiScanner::calcIntensity( mat, targetArea, radius +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords +#endif ); } diff --git a/src/scanner/MultiScanner.h b/src/scanner/MultiScanner.h index 2b9cf43b8..b96b731db 100644 --- a/src/scanner/MultiScanner.h +++ b/src/scanner/MultiScanner.h @@ -157,7 +157,6 @@ class MultiScanner : public Scanner{ */ void computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -165,12 +164,18 @@ class MultiScanner : public Scanner{ NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) override; /** * @see Scanner::initializeFullWaveform @@ -196,6 +201,9 @@ class MultiScanner : public Scanner{ double const targetArea, double const radius, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const override; /** * @see Scanner::calcIntensity diff --git a/src/scanner/Scanner.cpp b/src/scanner/Scanner.cpp index 1a6a2410f..aba242ff4 100644 --- a/src/scanner/Scanner.cpp +++ b/src/scanner/Scanner.cpp @@ -403,6 +403,9 @@ void Scanner::buildScanningPulseProcess( *randGen1, *randGen2, *intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pool->getPulseRecorder() +#endif ) ); } @@ -415,6 +418,9 @@ void Scanner::buildScanningPulseProcess( *randGen1, *randGen2, *intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pool->getPulseRecorder() +#endif ) ); } diff --git a/src/scanner/Scanner.h b/src/scanner/Scanner.h index c27a091c1..865fd588b 100644 --- a/src/scanner/Scanner.h +++ b/src/scanner/Scanner.h @@ -429,8 +429,6 @@ class Scanner : public Asset { * @brief Perform ray casting to find intersections * @param[in] handleSubray The function where computed subrays must be * delegated to - * @param[in] tMinMax Minimum and maximum time to intersection with respect - * to the axis aligned bounding box that bounds the scene * @param[out] reflections Where reflections must be stored when a hit is * registered * @param[out] intersects Where intersections must be stored when a hit is @@ -442,7 +440,6 @@ class Scanner : public Asset { */ virtual void computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -450,12 +447,18 @@ class Scanner : public Asset { NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) = 0; /** @@ -499,6 +502,9 @@ class Scanner : public Asset { double const targetArea, double const radius, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const = 0; /** * @brief Handle to which scanning device request the intensity computation diff --git a/src/scanner/ScanningDevice.cpp b/src/scanner/ScanningDevice.cpp index b39356ec3..6653c37b8 100644 --- a/src/scanner/ScanningDevice.cpp +++ b/src/scanner/ScanningDevice.cpp @@ -3,6 +3,10 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +using namespace helios::analytics; +#endif // *** CONSTRUCTION / DESTRUCTION *** // // ************************************ // @@ -214,7 +218,6 @@ Rotation ScanningDevice::calcExactAbsoluteBeamAttitude( void ScanningDevice::computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -222,12 +225,21 @@ void ScanningDevice::computeSubrays( NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >=2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - std::vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, std::vector &intersects +#if DATA_ANALYTICS >=2 + ,std::shared_ptr pulseRecorder +#endif ){ +#if DATA_ANALYTICS >=2 + bool subrayHit; +#endif int const beamSampleQuality = FWF_settings.beamSampleQuality; double const radiusStep_rad = beamDivergence_rad/beamSampleQuality; @@ -251,8 +263,12 @@ void ScanningDevice::computeSubrays( // # Loop over sub-rays along the circle for (int circleStep = 0; circleStep < circleSteps; circleStep++){ +#if DATA_ANALYTICS >=2 + std::vector subraySimRecord( + 17, std::numeric_limits::quiet_NaN() + ); +#endif handleSubray( - tMinMax, circleStep, circleStep_rad, r1, @@ -260,7 +276,20 @@ void ScanningDevice::computeSubrays( intersectionHandlingNoiseSource, reflections, intersects +#if DATA_ANALYTICS >=2 + ,subrayHit, + subraySimRecord +#endif ); +#if DATA_ANALYTICS >=2 + HDA_GV.incrementGeneratedSubraysCount(); + subraySimRecord[0] = (double) subrayHit; + subraySimRecord[1] = (double) radiusStep; + subraySimRecord[2] = (double) circleSteps; + subraySimRecord[3] = (double) circleStep; + subraySimRecord[4] = subrayDivergenceAngle_rad; + pulseRecorder->recordSubraySimulation(subraySimRecord); +#endif } } } @@ -322,20 +351,33 @@ double ScanningDevice::calcIntensity( Material const &mat, double const targetArea, double const radius +#if DATA_ANALYTICS >=2 + ,std::vector> &calcIntensityRecords +#endif ) const { - double bdrf = 0; + double bdrf = 0, sigma = 0; if(mat.isPhong()) { bdrf = mat.reflectance * EnergyMaths::phongBDRF( incidenceAngle, mat.specularity, mat.specularExponent ); + sigma = EnergyMaths::calcCrossSection( + bdrf, targetArea, incidenceAngle + ); } else if(mat.isLambert()){ bdrf = mat.reflectance * std::cos(incidenceAngle); + sigma = EnergyMaths::calcCrossSection( + bdrf, targetArea, incidenceAngle + ); } else if(mat.isUniform()){ bdrf = mat.reflectance; + sigma = EnergyMaths::calcCrossSection( + bdrf, targetArea, incidenceAngle + ); + if(sigma < 0) sigma = -sigma; } else{ std::stringstream ss; @@ -343,10 +385,7 @@ double ScanningDevice::calcIntensity( << mat.name << "\""; logging::ERR(ss.str()); } - double const sigma = EnergyMaths::calcCrossSection( - bdrf, targetArea, incidenceAngle - ); - return EnergyMaths::calcReceivedPower( + double const receivedPower = EnergyMaths::calcReceivedPower( averagePower_w, wavelength_m, targetRange, @@ -359,6 +398,21 @@ double ScanningDevice::calcIntensity( atmosphericExtinction, sigma ) * 1000000000.0; +#if DATA_ANALYTICS >= 2 + std::vector calcIntensityRecord( + 11, std::numeric_limits::quiet_NaN() + ); + calcIntensityRecord[3] = incidenceAngle; + calcIntensityRecord[4] = targetRange; + calcIntensityRecord[5] = targetArea; + calcIntensityRecord[6] = radius; + calcIntensityRecord[7] = bdrf; + calcIntensityRecord[8] = sigma; + calcIntensityRecord[9] = receivedPower; + calcIntensityRecord[10] = 0; // By default, assume the point isn't captured + calcIntensityRecords.push_back(calcIntensityRecord); +#endif + return receivedPower; } double ScanningDevice::calcIntensity( double const targetRange, diff --git a/src/scanner/ScanningDevice.h b/src/scanner/ScanningDevice.h index c30be2321..754727175 100644 --- a/src/scanner/ScanningDevice.h +++ b/src/scanner/ScanningDevice.h @@ -11,6 +11,10 @@ class AbstractDetector; #include #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_PulseRecorder; +#endif #include @@ -296,7 +300,6 @@ class ScanningDevice : public Asset { */ void computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -304,11 +307,17 @@ class ScanningDevice : public Asset { NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - std::vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, std::vector &intersects +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ); /** * @see Scanner::initializeFullWaveform @@ -377,6 +386,9 @@ class ScanningDevice : public Asset { Material const &mat, double const targetArea, double const radius +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const; /** diff --git a/src/scanner/SingleScanner.cpp b/src/scanner/SingleScanner.cpp index a7dd9dd7b..8b7245d7f 100644 --- a/src/scanner/SingleScanner.cpp +++ b/src/scanner/SingleScanner.cpp @@ -194,7 +194,6 @@ Rotation SingleScanner::calcAbsoluteBeamAttitude(size_t const idx) { } void SingleScanner::computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -202,19 +201,27 @@ void SingleScanner::computeSubrays( NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ){ scanDev.computeSubrays( handleSubray, - tMinMax, intersectionHandlingNoiseSource, reflections, intersects +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } @@ -248,6 +255,9 @@ double SingleScanner::calcIntensity( double const targetArea, double const radius, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const { return scanDev.calcIntensity( incidenceAngle, @@ -255,6 +265,9 @@ double SingleScanner::calcIntensity( mat, targetArea, radius +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords +#endif ); } double SingleScanner::calcIntensity( diff --git a/src/scanner/SingleScanner.h b/src/scanner/SingleScanner.h index a21b6043d..becee3d3a 100644 --- a/src/scanner/SingleScanner.h +++ b/src/scanner/SingleScanner.h @@ -130,7 +130,6 @@ class SingleScanner : public Scanner{ */ void computeSubrays( std::function const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -138,12 +137,18 @@ class SingleScanner : public Scanner{ NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif )> handleSubray, - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) override; /** * @see Scanner::initializeFullWaveform @@ -169,6 +174,9 @@ class SingleScanner : public Scanner{ double const targetArea, double const radius, size_t const idx +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords +#endif ) const override; /** * @see Scanner::calcIntensity diff --git a/src/scanner/WarehouseScanningPulseProcess.cpp b/src/scanner/WarehouseScanningPulseProcess.cpp index 561d20757..a27e0b949 100644 --- a/src/scanner/WarehouseScanningPulseProcess.cpp +++ b/src/scanner/WarehouseScanningPulseProcess.cpp @@ -11,6 +11,9 @@ WarehouseScanningPulseProcess::WarehouseScanningPulseProcess( RandomnessGenerator &randGen1, RandomnessGenerator &randGen2, UniformNoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) : ScanningPulseProcess(scanner), dropper(dropper), @@ -18,8 +21,11 @@ WarehouseScanningPulseProcess::WarehouseScanningPulseProcess( randGen1(randGen1), randGen2(randGen2), intersectionHandlingNoiseSource(intersectionHandlingNoiseSource) -{ - if(pool.getPoolSize() > 0){ // Parallel computation +#if DATA_ANALYTICS >= 2 + ,pulseRecorder(pulseRecorder) +#endif + { + if(pool.getPoolSize() > 0){ // Parallel computation handler = [&] (SimulatedPulse const &sp) -> void { handlePulseComputationParallel(sp); }; @@ -43,12 +49,20 @@ void WarehouseScanningPulseProcess::onLegComplete(){ randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); // Assist thread pool with pending tasks (on WarehouseThreadPool::join atm) shared_ptr task; while( (task=pool.get()) != nullptr){ - (*task)(apMatrix, randGen1, randGen2, intersectionHandlingNoiseSource); + (*task)( + apMatrix, randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif + ); } // Wait for threads to finish @@ -70,6 +84,9 @@ void WarehouseScanningPulseProcess::handlePulseComputationSequential( randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } @@ -93,6 +110,9 @@ void WarehouseScanningPulseProcess::handlePulseComputationParallel( randGen1, randGen2, intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } } diff --git a/src/scanner/WarehouseScanningPulseProcess.h b/src/scanner/WarehouseScanningPulseProcess.h index a7ebfe27c..d30b4c18f 100644 --- a/src/scanner/WarehouseScanningPulseProcess.h +++ b/src/scanner/WarehouseScanningPulseProcess.h @@ -3,6 +3,9 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +#endif #include @@ -20,6 +23,14 @@ class WarehouseScanningPulseProcess : public ScanningPulseProcess { public: #else protected: +#endif +#if DATA_ANALYTICS >= 2 + /** + * @brief The helios::analytics::PulseRecorder to be used to handle the + * records representing the computed pulse tasks. + * @see helios::analytics::PulseRecorder + */ + std::shared_ptr pulseRecorder; #endif // *** ATTRIBUTES *** // // ******************** // @@ -82,6 +93,9 @@ class WarehouseScanningPulseProcess : public ScanningPulseProcess { RandomnessGenerator &randGen1, RandomnessGenerator &randGen2, UniformNoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ); virtual ~WarehouseScanningPulseProcess() = default; diff --git a/src/scanner/detector/AbstractPulseRunnable.cpp b/src/scanner/detector/AbstractPulseRunnable.cpp index 6bae9fd0b..278504588 100644 --- a/src/scanner/detector/AbstractPulseRunnable.cpp +++ b/src/scanner/detector/AbstractPulseRunnable.cpp @@ -62,6 +62,10 @@ void AbstractPulseRunnable::capturePoint( std::mutex *allMeasurementsMutex, std::vector *cycleMeasurements, std::mutex *cycleMeasurementsMutex +#if DATA_ANALYTICS >= 2 + ,std::vector &calcIntensityRecord, + std::shared_ptr pulseRecorder +#endif ) { // Abort if point distance is below mininum scanner range: // TODO Pending : This check is already done in FullWaveformPulseRunnable @@ -77,6 +81,10 @@ void AbstractPulseRunnable::capturePoint( // TODO Pending : Is it necessary to compute position again? Notice it is // known from ray intersection point at FullWaveformPulseRunnable m.position = m.beamOrigin + m.beamDirection * m.distance; +#if DATA_ANALYTICS >= 2 + calcIntensityRecord[10] = 1; + pulseRecorder->recordIntensityCalculation(calcIntensityRecord); +#endif if(allMeasurements != nullptr){ std::unique_lock lock(*allMeasurementsMutex); allMeasurements->push_back(m); diff --git a/src/scanner/detector/AbstractPulseRunnable.h b/src/scanner/detector/AbstractPulseRunnable.h index 90727e758..91052c094 100644 --- a/src/scanner/detector/AbstractPulseRunnable.h +++ b/src/scanner/detector/AbstractPulseRunnable.h @@ -7,6 +7,9 @@ class Measurement; #include "LasSpecification.h" class Scanner; #include +#if DATA_ANALYTICS >= 2 +#include +#endif #include @@ -90,6 +93,10 @@ class AbstractPulseRunnable : public PulseTask{ std::mutex *allMeasurementsMutex, std::vector *cycleMeasurements, std::mutex *cycleMeasurementsMutex +#if DATA_ANALYTICS >= 2 + ,std::vector &calcIntensityRecord, + std::shared_ptr pulseRecorder +#endif ); /** diff --git a/src/scanner/detector/FullWaveformPulseRunnable.cpp b/src/scanner/detector/FullWaveformPulseRunnable.cpp index dc799aa3e..be11e3c12 100644 --- a/src/scanner/detector/FullWaveformPulseRunnable.cpp +++ b/src/scanner/detector/FullWaveformPulseRunnable.cpp @@ -9,13 +9,20 @@ #include "maths/Rotation.h" #include -#include "MarquardtFitter.h" #include #include #include #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_GV; + +#include +#include +#endif + using namespace std; // *** CONSTANTS *** // @@ -30,6 +37,9 @@ void FullWaveformPulseRunnable::operator()( RandomnessGenerator &randGen, RandomnessGenerator &randGen2, NoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ){ // Deferred/lazy initialization initialize(); @@ -42,24 +52,36 @@ void FullWaveformPulseRunnable::operator()( // performance optimization, so we should keep it nevertheless. sbecht 2016-04-24 // Early abort if central axis of the beam does not intersect with the scene: - vector tMinMax = scene.getAABB()->getRayIntersection( + vector const tMinMax = scene.getAABB()->getRayIntersection( pulse.getOriginRef(), beamDir ); - if (tMinMax.empty()) { +#if DATA_ANALYTICS >= 2 + HDA_GV.incrementGeneratedRaysBeforeEarlyAbortCount(); +#endif + if (checkEarlyAbort(tMinMax)) { logging::DEBUG("Early abort - beam does not intersect with the scene"); scanner->setLastPulseWasHit(false, pulse.getDeviceIndex()); return; } +#if DATA_ANALYTICS >= 2 + HDA_GV.incrementGeneratedRaysAfterEarlyAbortCount(); +#endif // Ray casting (find intersections) map reflections; vector intersects; +#if DATA_ANALYTICS >= 2 + std::vector> calcIntensityRecords; +#endif computeSubrays( - tMinMax, intersectionHandlingNoiseSource, reflections, intersects +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords, + pulseRecorder +#endif ); // Digest intersections @@ -70,6 +92,10 @@ void FullWaveformPulseRunnable::operator()( beamDir, reflections, intersects +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords, + pulseRecorder +#endif ); } @@ -85,14 +111,16 @@ void FullWaveformPulseRunnable::initialize(){ ); } void FullWaveformPulseRunnable::computeSubrays( - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ){ scanner->computeSubrays( [&] ( - vector const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -100,9 +128,12 @@ void FullWaveformPulseRunnable::computeSubrays( NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects - ) -> void{ +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord +#endif + ) -> void { handleSubray( - _tMinMax, circleStep, circleStep_rad, r1, @@ -110,18 +141,24 @@ void FullWaveformPulseRunnable::computeSubrays( intersectionHandlingNoiseSource, reflections, intersects +#if DATA_ANALYTICS >= 2 + ,subrayHit, + subraySimRecord, + calcIntensityRecords +#endif ); }, - tMinMax, intersectionHandlingNoiseSource, reflections, intersects, pulse.getDeviceIndex() +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } void FullWaveformPulseRunnable::handleSubray( - vector const &_tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -129,14 +166,44 @@ void FullWaveformPulseRunnable::handleSubray( NoiseSource &intersectionHandlingNoiseSource, map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord, + std::vector> &calcIntensityRecords +#endif ){ +#if DATA_ANALYTICS >= 2 + subrayHit = false; +#endif // Rotate around the circle: - vector tMinMax = _tMinMax; Rotation r2 = Rotation(Directions::forward, circleStep_rad * circleStep); r2 = r2.applyTo(r1); glm::dvec3 subrayDirection = pulse.getAttitude().applyTo(r2) .applyTo(Directions::forward); + vector tMinMax = scene.getAABB()->getRayIntersection( + pulse.getOriginRef(), + subrayDirection + ); +#if DATA_ANALYTICS >= 2 + glm::dvec3 rayDirection = pulse.computeDirection(); + subraySimRecord[5] = glm::l2Norm(rayDirection); // Ray norm + subraySimRecord[6] = glm::l2Norm(subrayDirection); // Subray norm + subraySimRecord[7] = glm::angle( // Angle between ray and subray + rayDirection, + subrayDirection + ); + subraySimRecord[8] = (rayDirection[0] < 0) == (subrayDirection[0] < 0); + subraySimRecord[9] = tMinMax[0]; + subraySimRecord[10] = tMinMax[1]; + subraySimRecord[11] = subrayDirection.x; + subraySimRecord[12] = subrayDirection.y; + subraySimRecord[13] = subrayDirection.z; + subraySimRecord[14] = rayDirection.x; + subraySimRecord[15] = rayDirection.y; + subraySimRecord[16] = rayDirection.z; +#endif + if(checkEarlyAbort(tMinMax)) return; glm::dvec3 subrayOrigin(pulse.getOrigin()); bool rayContinues = true; @@ -150,8 +217,12 @@ void FullWaveformPulseRunnable::handleSubray( ); if (intersect != nullptr && intersect->prim != nullptr) { +#if DATA_ANALYTICS >= 2 + HDA_GV.incrementSubrayIntersectionCount(); + subrayHit = true; +#endif // Incidence angle: - if(!scanner->isFixedIncidenceAngle()) { + if (!scanner->isFixedIncidenceAngle()) { incidenceAngle = intersect->prim->getIncidenceAngle_rad( pulse.getOriginRef(), @@ -167,7 +238,7 @@ void FullWaveformPulseRunnable::handleSubray( ); // Distance must be inside [rangeMin, rangeMax] interval - if(detector->isDistanceNotInRange(distance)) continue; + if (detector->isDistanceNotInRange(distance)) continue; // Distance between beam's center line and intersection point: double const radius = sin(divergenceAngle) * distance; @@ -175,7 +246,7 @@ void FullWaveformPulseRunnable::handleSubray( distance, pulse.getDeviceIndex() ); double intensity = 0.0; - if(intersect->prim->canComputeSigmaWithLadLut()){ + if (intersect->prim->canComputeSigmaWithLadLut()) { // LadLut based intensity computation double sigma = intersect->prim->computeSigmaWithLadLut( subrayDirection @@ -183,9 +254,11 @@ void FullWaveformPulseRunnable::handleSubray( intensity = scanner->calcIntensity( distance, radius, sigma, pulse.getDeviceIndex() ); - } - else{ + } else { // Lighting-based intensity computation +#if DATA_ANALYTICS >= 2 + HDA_GV.incrementIntensityComputationsCount(); +#endif intensity = scanner->calcIntensity( incidenceAngle, distance, @@ -193,12 +266,15 @@ void FullWaveformPulseRunnable::handleSubray( targetArea, radius, pulse.getDeviceIndex() +#if DATA_ANALYTICS >= 2 + , calcIntensityRecords +#endif ); } // Intersection handling - if(intersect->prim->canHandleIntersections()) { + if (intersect->prim->canHandleIntersections()) { glm::dvec3 outsideIntersectionPoint = RayUtils::obtainPointAfterTraversing( *intersect->prim->getAABB(), @@ -223,32 +299,49 @@ void FullWaveformPulseRunnable::handleSubray( subrayDirection ); rayContinues = true; - } - else{ // Update distance considering noise + } else { // Update distance considering noise distance = glm::distance( ihr.getIntersectionPoint(), pulse.getOriginRef() ); } } - if(!rayContinues) { // If ray is not continuing + if (!rayContinues) { // If ray is not continuing // Then register hit by default reflections.insert( pair(distance, intensity) ); intersects.push_back(*intersect); } +#if DATA_ANALYTICS >= 2 + std::vector &calcIntensityRecord = + calcIntensityRecords.back(); + calcIntensityRecord[0] = intersect->point.x; + calcIntensityRecord[1] = intersect->point.y; + calcIntensityRecord[2] = intersect->point.z; + } + else { + HDA_GV.incrementSubrayNonIntersectionCount(); +#endif } } +#if DATA_ANALYTICS >= 2 + if(subrayHit) HDA_GV.incrementIntersectiveSubraysCount(); + else HDA_GV.incrementNonIntersectiveSubraysCount(); +#endif } void FullWaveformPulseRunnable::digestIntersections( - std::vector>& apMatrix, + std::vector> &apMatrix, RandomnessGenerator &randGen, RandomnessGenerator &randGen2, glm::dvec3 &beamDir, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ){ // Find maximum hit distances double minHitDist_m, maxHitDist_m; @@ -257,6 +350,11 @@ void FullWaveformPulseRunnable::digestIntersections( // If nothing was hit, get out of here if (maxHitDist_m < 0) { scanner->setLastPulseWasHit(false, pulse.getDeviceIndex()); +#if DATA_ANALYTICS >= 2 + if(!calcIntensityRecords.empty()) { + pulseRecorder->recordIntensityCalculation(calcIntensityRecords); + } +#endif return; } @@ -272,7 +370,12 @@ void FullWaveformPulseRunnable::digestIntersections( distanceThreshold, peakIntensityIndex, numFullwaveBins - )) return; + )){ +#if DATA_ANALYTICS >= 2 + pulseRecorder->recordIntensityCalculation(calcIntensityRecords); +#endif + return; + } std::vector fullwave(numFullwaveBins, 0.0); // Populate full waveform @@ -299,6 +402,10 @@ void FullWaveformPulseRunnable::digestIntersections( numFullwaveBins, peakIntensityIndex, minHitTime_ns +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords, + pulseRecorder +#endif ); // Export measurements and full waveform data @@ -311,6 +418,10 @@ void FullWaveformPulseRunnable::digestIntersections( maxHitTime_ns, randGen, randGen2 +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords, + pulseRecorder +#endif ); } @@ -404,6 +515,10 @@ void FullWaveformPulseRunnable::digestFullWaveform( int const numFullwaveBins, int const peakIntensityIndex, double const minHitTime_ns +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ){ // Extract points from waveform data via Gaussian decomposition size_t const devIdx = pulse.getDeviceIndex(); @@ -419,31 +534,20 @@ void FullWaveformPulseRunnable::digestFullWaveform( double echo_width = 0.0; - for (int i = 0; i < numFullwaveBins; ++i) { - if (fullwave[i] < eps) continue; +#if DATA_ANALYTICS >= 2 + std::unordered_set capturedIndices; +#endif - // Continue with next iteration if there is no peak - if(!detectPeak(i, win_size, fullwave)) continue; - - // Gaussian model fitting - if (scanner->isCalcEchowidth()) { - try { - fit.setParameters( - vector{0, fullwave[i], (double) i, 1} - ); - fit.fitData(); - } - catch (std::exception &e) { - logging::ERR(e.what()); - continue; - } - echo_width = fit.getParameters()[3]; - echo_width = echo_width * nsPerBin; - - if (echo_width < 0.1) { // TODO Pending : 0.1 to threshold variable - continue; - } - } + for (int i = 0; i < numFullwaveBins; ++i) { + // Skip iteration if the handling of the bin_i does not accept the hit + if(handleFullWaveformBin( + fullwave, + fit, + echo_width, + i, + win_size, + nsPerBin + )) continue; // Compute distance double distance = SPEEDofLIGHT_mPerNanosec * @@ -453,6 +557,10 @@ void FullWaveformPulseRunnable::digestFullWaveform( double minDifference = numeric_limits::max(); shared_ptr closestIntersection = nullptr; +#ifdef DATA_ANALYTICS + size_t intersectionIdx = 0; + size_t closestIntersectionIdx = 0; +#endif for (RaySceneIntersection intersect : intersects) { double const intersectDist = glm::distance( intersect.point, pulse.getOriginRef()); @@ -462,7 +570,13 @@ void FullWaveformPulseRunnable::digestFullWaveform( minDifference = absDistDiff; closestIntersection = make_shared( intersect); +#ifdef DATA_ANALYTICS + closestIntersectionIdx = intersectionIdx; +#endif } +#ifdef DATA_ANALYTICS + ++intersectionIdx; +#endif } string hitObject; @@ -501,11 +615,65 @@ void FullWaveformPulseRunnable::digestFullWaveform( pointsMeasurement.push_back(tmp); ++numReturns; +#if DATA_ANALYTICS >= 2 + capturedIndices.insert(closestIntersectionIdx); +#endif // Check if maximum number of returns per pulse has been reached if(!scanner->checkMaxNOR(numReturns, devIdx)) break; } +#if DATA_ANALYTICS >= 2 + // Record all non captured points and remove them from records + size_t const numRecords = calcIntensityRecords.size(); + size_t nonCapturedCount = 0; + for(size_t i = 0 ; i < numRecords ; ++i){ + if(capturedIndices.find(i) == capturedIndices.end()){ + pulseRecorder->recordIntensityCalculation( + calcIntensityRecords[i-nonCapturedCount] + ); + calcIntensityRecords.erase( + calcIntensityRecords.begin()+i-nonCapturedCount + ); + ++nonCapturedCount; + } + } +#endif +} + +bool FullWaveformPulseRunnable::handleFullWaveformBin( + std::vector const &fullwave, + MarquardtFitter &fit, + double &echoWidth, + int const binIndex, + int const winSize, + double const nsPerBin +){ + if (fullwave[binIndex] < eps) return true; + + // Continue with next iteration if there is no peak + if(!detectPeak(binIndex, winSize, fullwave)) return true; + // Gaussian model fitting + if (scanner->isCalcEchowidth()) { + try { + fit.setParameters( + vector{0, fullwave[binIndex], (double) binIndex, 1} + ); + fit.fitData(); + } + catch (std::exception &e) { + logging::ERR(e.what()); + return true; + } + echoWidth = fit.getParameters()[3]; + echoWidth = echoWidth * nsPerBin; + + if (echoWidth < 0.1) { // TODO Pending : 0.1 to threshold variable + return true; + } + } + // Return false as the bin corresponds to an accepted hit + return false; } void FullWaveformPulseRunnable::exportOutput( @@ -517,9 +685,16 @@ void FullWaveformPulseRunnable::exportOutput( double const maxHitTime_ns, RandomnessGenerator &randGen, RandomnessGenerator &randGen2 +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ){ // ############ END Extract points from waveform data ################ if (numReturns > 0) { +#if DATA_ANALYTICS >= 2 + size_t i = 0; +#endif for (Measurement & pm : pointsMeasurement) { pm.pulseReturnNumber = numReturns; capturePoint( @@ -529,7 +704,14 @@ void FullWaveformPulseRunnable::exportOutput( scanner->allMeasurementsMutex.get(), scanner->cycleMeasurements.get(), scanner->cycleMeasurementsMutex.get() +#if DATA_ANALYTICS >= 2 + ,calcIntensityRecords[i], + pulseRecorder +#endif ); +#if DATA_ANALYTICS >= 2 + ++i; +#endif } if(scanner->isWriteWaveform()) { scanner->setLastPulseWasHit(true, pulse.getDeviceIndex()); @@ -606,7 +788,7 @@ bool FullWaveformPulseRunnable::detectPeak( } } - int size = fullwave.size(); + int const size = fullwave.size(); for (int j = std::min(size, i + 1); j < std::min(size, i + win_size); j++){ if (fullwave[j] < eps || fullwave[j] >= fullwave[i]) { return false; diff --git a/src/scanner/detector/FullWaveformPulseRunnable.h b/src/scanner/detector/FullWaveformPulseRunnable.h index ab5ea00a3..3dcba2158 100644 --- a/src/scanner/detector/FullWaveformPulseRunnable.h +++ b/src/scanner/detector/FullWaveformPulseRunnable.h @@ -1,13 +1,17 @@ #pragma once #include -#include "AbstractDetector.h" -#include "AbstractPulseRunnable.h" -#include "FullWaveformPulseDetector.h" -#include "RaySceneIntersection.h" -#include "ScenePart.h" +#include +#include +#include +#include +#include +#include #include #include +#if DATA_ANALYTICS >= 2 +#include +#endif #include @@ -79,10 +83,13 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { * @see FullWaveformPulseRunnable::handleSubray */ void computeSubrays( - vector const &tMinMax, NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ); /** * @brief Handle sub-rays along the circle @@ -96,7 +103,6 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { * @see FullWaveformPulseRunnable::computeSubrays */ void handleSubray( - vector const &tMinMax, int const circleStep, double const circleStep_rad, Rotation &r1, @@ -104,6 +110,11 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { NoiseSource &intersectionHandlingNoiseSource, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,bool &subrayHit, + std::vector &subraySimRecord, + std::vector> &calcIntensityRecords +#endif ); /** * @brief Digest intersections found through ray casting @@ -126,6 +137,10 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { glm::dvec3 &beamDir, std::map &reflections, vector &intersects +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ); /** * @brief Find min and max hit distances in meters @@ -197,6 +212,30 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { int const numFullwaveBins, int const peakIntensityIndex, double const minHitTime_ns +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif + ); + /** + * @brief Handle the bin of the fullwave that corresponds to the given + * bin index. + * @param fullwave The fullwave vector. + * @param fit Reference to the MarquardFitter object used to handle the + * fullwave. + * @param echoWidth Output variable to store the echo width. + * @param binIndex Index of the bin to be handled. + * @return True when the handled bin does not correspond to an accepted + * hit, i.e., true to skip the current iteration. False to proceed with + * the computation of the current iteration. + */ + bool handleFullWaveformBin( + std::vector const &fullwave, + MarquardtFitter &fit, + double &echoWidth, + int const binIndex, + int const winSize, + double const nsPerBin ); /** * @brief Export measurements and full waveform data @@ -213,6 +252,10 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { double const maxHitTime_ns, RandomnessGenerator &randGen, RandomnessGenerator &randGen2 +#if DATA_ANALYTICS >= 2 + ,std::vector> &calcIntensityRecords, + std::shared_ptr pulseRecorder +#endif ); // *** ASSISTANCE METHODS *** // @@ -265,6 +308,28 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { RandomnessGenerator &rg2 ); + /** + * @brief Check whether the ray/subray must be aborted or not depending + * on its intersection times. + * + * NOTE that typically the tMinMax vector is expected to come from a + * intersection computation on the minimum axis aligned bounding box + * containing the scene, i.e., the Scene::bbox attribute. + * + * @param tMinMax The first component is the minimum intersection time + * \f$t_*\f$, the second component is the maximum intersection time + * \f$t^*\f$, and no components means no intersection. See + * AABB::getRayIntersection for more details because an early abort check + * is a ray-AABB (Axis-Aligned Bounding-Box) check. + * @return True if the ray fails to intersects with the axis-aligned + * bounding box, False otherwise. + * @see Scene + * @see AABB + */ + inline bool checkEarlyAbort(std::vector const &tMinMax){ + return tMinMax.empty() || (tMinMax[0] < 0 && tMinMax[1] <= 0); + } + public: // *** O P E R A T O R *** // // ************************* // @@ -293,5 +358,8 @@ class FullWaveformPulseRunnable : public AbstractPulseRunnable { RandomnessGenerator &randGen, RandomnessGenerator &randGen2, NoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) override; }; \ No newline at end of file diff --git a/src/scanner/detector/PulseTask.h b/src/scanner/detector/PulseTask.h index 182c30c94..23762d215 100644 --- a/src/scanner/detector/PulseTask.h +++ b/src/scanner/detector/PulseTask.h @@ -2,6 +2,10 @@ #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_PulseRecorder; +#endif #include @@ -33,5 +37,8 @@ class PulseTask{ RandomnessGenerator &randGen, RandomnessGenerator &randGen2, NoiseSource &intersectionHandlingNoiseSource +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr pulseRecorder +#endif ) = 0; }; \ No newline at end of file diff --git a/src/scanner/detector/PulseTaskDropper.h b/src/scanner/detector/PulseTaskDropper.h index a6a1cb3ff..444bbb0ff 100644 --- a/src/scanner/detector/PulseTaskDropper.h +++ b/src/scanner/detector/PulseTaskDropper.h @@ -6,6 +6,10 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_PulseRecorder; +#endif /** @@ -24,6 +28,9 @@ class PulseTaskDropper : public BuddingTaskDropper< RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >{ public: // *** CONSTRUCTION / DESTRUCTION *** // @@ -49,6 +56,9 @@ class PulseTaskDropper : public BuddingTaskDropper< RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >(maxTasks, delta1, initDelta1, delta2, lastSign) {} virtual ~PulseTaskDropper() = default; @@ -62,6 +72,9 @@ class PulseTaskDropper : public BuddingTaskDropper< RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >::drop; // To avoid overriding hides drop overloads using TaskDropper< PulseTask, @@ -70,6 +83,9 @@ class PulseTaskDropper : public BuddingTaskDropper< RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >::tryDrop; // To avoid override hides tryDrop overloads /** * @brief Like TaskDropper::drop but dropping all pulse tasks through a diff --git a/src/scanner/detector/PulseThreadPool.h b/src/scanner/detector/PulseThreadPool.h index f541aae0c..979aa3265 100644 --- a/src/scanner/detector/PulseThreadPool.h +++ b/src/scanner/detector/PulseThreadPool.h @@ -5,6 +5,9 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +#endif /** * @version 1.0 @@ -17,6 +20,9 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >, public PulseThreadPoolInterface { @@ -24,6 +30,9 @@ class PulseThreadPool : public: #else protected: +#endif +#if DATA_ANALYTICS >= 2 + std::shared_ptr pulseRecorder; #endif // *** ATTRIBUTES *** // // ******************** // @@ -86,10 +95,18 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >(_pool_size), dynamic(dynamic) { // Allocate +#if DATA_ANALYTICS >= 2 + pulseRecorder = std::make_shared( + "helios_pulse_records" + ); +#endif apMatrices = new std::vector>[this->pool_size]; randGens = new RandomnessGenerator[this->pool_size]; randGens2 = new RandomnessGenerator[this->pool_size]; @@ -110,6 +127,10 @@ class PulseThreadPool : } virtual ~PulseThreadPool(){ +#if DATA_ANALYTICS >= 2 + // Flush and close pulse recorder + this->pulseRecorder->closeBuffers(); +#endif // Release memory delete[] apMatrices; delete[] randGens; @@ -130,6 +151,9 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) override { run_res_task(dropper); @@ -145,6 +169,9 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) override { return try_run_res_task(dropper); @@ -158,9 +185,21 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif >::join(); } +#if DATA_ANALYTICS >= 2 + /** + * @see PulseThreadPoolInterface::getPulseRecorder + */ + std::shared_ptr getPulseRecorder() override{ + return pulseRecorder; + } +#endif + protected: // *** M E T H O D S *** // // *********************** // @@ -175,6 +214,9 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif )> &task, int const resourceIdx ) override{ @@ -183,6 +225,9 @@ class PulseThreadPool : randGens[resourceIdx], randGens2[resourceIdx], intersectionHandlingNoiseSources[resourceIdx] +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } /** @@ -196,6 +241,9 @@ class PulseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif )> &task, int const resourceIdx ) override { diff --git a/src/scanner/detector/PulseThreadPoolInterface.h b/src/scanner/detector/PulseThreadPoolInterface.h index e3d7a6b76..a42148baa 100644 --- a/src/scanner/detector/PulseThreadPoolInterface.h +++ b/src/scanner/detector/PulseThreadPoolInterface.h @@ -4,6 +4,11 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_PulseRecorder; +#endif + #include @@ -37,6 +42,9 @@ class PulseThreadPoolInterface{ RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) = 0; /** @@ -59,6 +67,9 @@ class PulseThreadPoolInterface{ RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) = 0; /** @@ -66,4 +77,13 @@ class PulseThreadPoolInterface{ */ virtual void join() = 0; +#if DATA_ANALYTICS >= 2 + /** + * @brief Obtain the pointer to the pulse recorder being used by the + * pulse thread pool. + * @return The pulse recorder used by the pulse thread pool. + * @see helios::analytics::HDA_PulseRecorder + */ + virtual std::shared_ptr getPulseRecorder() = 0; +#endif }; \ No newline at end of file diff --git a/src/scanner/detector/PulseWarehouseThreadPool.h b/src/scanner/detector/PulseWarehouseThreadPool.h index 24a601475..b72803521 100644 --- a/src/scanner/detector/PulseWarehouseThreadPool.h +++ b/src/scanner/detector/PulseWarehouseThreadPool.h @@ -5,6 +5,9 @@ #include #include #include +#if DATA_ANALYTICS >= 2 +#include +#endif #include @@ -25,6 +28,14 @@ class PulseWarehouseThreadPool : public: #else protected: +#endif +#if DATA_ANALYTICS >= 2 + /** + * @brief The helios::analytics::PulseRecorder to be used to handle the + * records representing the computed pulse tasks. + * @see helios::analytics::PulseRecorder + */ + std::shared_ptr pulseRecorder; #endif // *** ATTRIBUTES *** // // ******************** // @@ -70,6 +81,11 @@ class PulseWarehouseThreadPool : randGens2 = new RandomnessGenerator[this->pool_size]; intersectionHandlingNoiseSources = new UniformNoiseSource[this->pool_size]; +#if DATA_ANALYTICS >= 2 + pulseRecorder = std::make_shared( + "helios_pulse_records" + ); +#endif // Initialize for (std::size_t i = 0; i < this->pool_size; ++i){ @@ -85,6 +101,10 @@ class PulseWarehouseThreadPool : } virtual ~PulseWarehouseThreadPool(){ +#if DATA_ANALYTICS >= 2 + // Flush and close pulse recorder + this->pulseRecorder->closeBuffers(); +#endif // Release memory delete[] apMatrices; delete[] randGens; @@ -105,6 +125,9 @@ class PulseWarehouseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) override { throw HeliosException( @@ -123,6 +146,9 @@ class PulseWarehouseThreadPool : RandomnessGenerator&, RandomnessGenerator&, NoiseSource& +#if DATA_ANALYTICS >= 2 + ,std::shared_ptr +#endif > &dropper ) override { return post(make_shared( @@ -135,6 +161,14 @@ class PulseWarehouseThreadPool : inline void join() override{ WarehouseThreadPool::join(); } +#if DATA_ANALYTICS >= 2 + /** + * @see PulseThreadPoolInterface::getPulseRecorder + */ + std::shared_ptr getPulseRecorder() override{ + return pulseRecorder; + } +#endif protected: // *** WAREHOUSE THREADPOOL *** // @@ -154,6 +188,9 @@ class PulseWarehouseThreadPool : randGens[tid], randGens2[tid], intersectionHandlingNoiseSources[tid] +#if DATA_ANALYTICS >= 2 + ,pulseRecorder +#endif ); } }; \ No newline at end of file diff --git a/src/scene/Scene.cpp b/src/scene/Scene.cpp index b7d943d39..f9eee51c5 100644 --- a/src/scene/Scene.cpp +++ b/src/scene/Scene.cpp @@ -23,6 +23,11 @@ using namespace glm; #include #include +#if DATA_ANALYTICS >= 2 +#include +using helios::analytics::HDA_GV; +#endif + using SurfaceInspector::maths::Plane; using SurfaceInspector::maths::PlaneFitter; @@ -182,6 +187,9 @@ std::shared_ptr Scene::getIntersection( ) const{ if (tMinMax.empty()) { logging::DEBUG("tMinMax is empty"); +#if DATA_ANALYTICS >= 2 + HDA_GV.incrementNonIntersectiveSubraysDueToNullTimeCount(); +#endif return nullptr; } return shared_ptr(raycaster->search( diff --git a/src/scene/primitives/AABB.cpp b/src/scene/primitives/AABB.cpp index 56d2bb250..8154d7a14 100644 --- a/src/scene/primitives/AABB.cpp +++ b/src/scene/primitives/AABB.cpp @@ -163,14 +163,15 @@ std::vector AABB::getRayIntersection( const int ysign = dir.y < 0; const int zsign = dir.z < 0; + // As long as divide by zero yields +-infinity, this logic should work double tmin, tmax; tmin = (bounds[xsign].x - orig.x) / dir.x; tmax = (bounds[1 - xsign].x - orig.x) / dir.x; const double tymin = (bounds[ysign].y - orig.y) / dir.y; const double tymax = (bounds[1 - ysign].y - orig.y) / dir.y; - if (tmin > tymax || tymin > tmax) return std::vector(); - if (tymin > tmin)tmin = tymin; + if (tmin > tymax || tymin > tmax) return {}; + if (tymin > tmin) tmin = tymin; if (tymax < tmax) tmax = tymax; const double tzmin = (bounds[zsign].z - orig.z) / dir.z; diff --git a/src/scene/primitives/AABB.h b/src/scene/primitives/AABB.h index 7fc6478ca..c1e06ddd8 100644 --- a/src/scene/primitives/AABB.h +++ b/src/scene/primitives/AABB.h @@ -128,6 +128,27 @@ class AABB : public Primitive { const glm::dvec3& intersectionPoint ) override; /** + * WARNING! The returned + * intersection times should be interpreted as follows, where \f$t_*\f$ + * is the minimum time and \f$t^*\f$ the maximum: + * + * If \f$t_* \geq 0, t^* \geq 0\f$ then the ray intersects the box twice, + * one time to enter the box, another time to leave it. + * + * If \f$t_* < 0, t^* < 0\f$ then the ray does not intersect the + * box, the line will do, but not the ray which is not allowed to + * go backward (see Convex Optimization by Stephen Boyd and Lieven + * Vandenberghe, where a ray is defined as + * \f$\left\{o + tv : t \geq 0\right\}\f$). + * + * If either \f$t_* < 0\f$ or \f$t^* < 0\f$ then the ray intersects the + * box once, it starts inside the box and the intersection occurs when + * the ray leaves the box. + * + * If the returned vector is empty, then there is no line-box intersection + * at all. + * + * * @see Primitive::getRayIntersection */ std::vector getRayIntersection(