Skip to content

Commit

Permalink
Doc and API updates (#20)
Browse files Browse the repository at this point in the history
* Add a new page about other modelling software

* add more modelling software

* add more software

* Revert "add more software"

This reverts commit 0b9bf6c.

* Resources dir for docs, better coordinate system figure

* Simplify API for BenchmarkData and ReferenceModels

* Update to work with API changes
  • Loading branch information
gavinmacaulay authored Aug 20, 2024
1 parent ed61748 commit bb5f2ab
Show file tree
Hide file tree
Showing 10 changed files with 235 additions and 50 deletions.
11 changes: 10 additions & 1 deletion docs/contributing.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,16 @@ A single coordinate system should be used for all models provided by echoSMs. Th

The right-handed spherical coordinate system as defined by ISO 80000-2[1] is to be used, where organisms are oriented such that θ corresponds to pitch and ɸ to roll, as illustrated below, where the organism lies along the z-axis and the positive x-axis extends above the organism:

![coordinate system](coordinate_system.jpg)

<!--- This code will include an html file, originally used to
include a live 3D view of the coordinate system, but there are
issues with the html so for the moment a 2D image is used.
<p align="center">
<iframe src="../coordinate_system2.html" title="Coordinate system" width="100%" height="500" frameborder="0"></iframe>
</p>
--->

![The coordinate system](resources/coordinate_system.svg){:style="height:400px;width400px"}

The definitions are such that θ values of 0°, 90°, and 180° correspond to organism pitches of head on, dorsal, and tail on, respectively, and positive values of ɸ indicate a roll to starboard. All model code should accept angles and produce results in this coordinate system. If the model calculations use a different coordinate system, the code should internally convert between the system given above and the version used in the code.

Expand Down
Binary file removed docs/coordinate_system.jpg
Binary file not shown.
23 changes: 23 additions & 0 deletions docs/resources/coordinate_system.svg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
File renamed without changes
Binary file added docs/resources/herring.stl
Binary file not shown.
167 changes: 167 additions & 0 deletions docs/src/make_coordinate_system_figure.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
# %%
"""Create the coordinate system figure for the documentation."""
import pyvista as pv
import numpy as np
from pathlib import Path

# These two functions came from https://github.com/pyvista/pyvista/discussions/5023 and are
# used to create arced arrows.


def rotate_to(a, b):
"""Return a rotation matrix taking the vector A to B."""
# Naive approach is unstable if a and b are near parallel
theta = np.arccos(np.dot(a, b))
eps = 1e-3
if np.absolute(theta) < eps:
return np.eye(4)
elif np.absolute(np.pi - theta) < eps: # Close to 180 degrees
# Choose the coordinate axis most orthogonal to A
x = np.zeros(3)
x[np.argmin(np.absolute(a))] = 1.0
axis = np.cross(a, x)
axis /= np.linalg.norm(axis)
return pv.utilities.transformations.axis_angle_rotation(axis, theta, deg=False)
else:
axis = np.cross(a, b)
return pv.utilities.transformations.axis_angle_rotation(axis, theta, deg=False)


def semi_circular_arrow(
start_angle=0,
circ_frac=0.75,
body_axial_res=100,
body_radial_res=50,
head_radial_res=100,
circ_radius=10,
body_radius=1,
head_length=3,
head_radius_frac=1.5,
normal=None,
center=None):
"""Create a semi circular arrow."""
t = np.linspace(0, circ_frac * 2 * np.pi, body_axial_res) + start_angle
x = circ_radius * np.cos(t)
y = circ_radius * np.sin(t)
z = np.zeros(body_axial_res)
body_pts = np.column_stack([x, y, z])
body = pv.MultipleLines(body_pts).tube(body_radius, n_sides=body_radial_res)

# Direction the head points
dhead = body_pts[-1] - body_pts[-2]
dhead /= np.linalg.norm(dhead)

head = pv.Cone(
center=body_pts[-1] + dhead * (head_length / 2),
direction=dhead,
height=head_length,
radius=head_radius_frac * body_radius,
resolution=head_radial_res,
)

arrow = body.merge(head, merge_points=False)

if normal is not None:
arrow = arrow.transform(rotate_to((0, 0, 1), normal), inplace=False)

if center is not None:
arrow = arrow.translate(center, inplace=False)

return arrow

# %%
# This code generates the figure


resourcesDir = Path(r'../resources')
al = 10.0 # axes length
sr = 0.01 # shaft radius for axes and arrows

# Axes arrows
arrow_x = pv.Arrow(start=(0, 0, 0), direction=(al, 0, 0), shaft_radius=sr, tip_radius=0.02,
tip_length=0.1, scale='auto')
arrow_y = pv.Arrow(start=(0, 0, 0), direction=(0, al, 0), shaft_radius=sr, tip_radius=0.02,
tip_length=0.1, scale='auto')
arrow_z = pv.Arrow(start=(0, 0, 0), direction=(0, 0, al), shaft_radius=sr, tip_radius=0.02,
tip_length=0.1, scale='auto')

# Angle arrows and arcs
arrow_theta = pv.Arrow(start=(0, 0, 0), direction=(al/2, 0, al/2),
shaft_radius=sr/2, tip_radius=0.02, tip_length=0.1, scale='auto')
arc_theta = semi_circular_arrow(center=(0, 0, 0), circ_frac=38/360, start_angle=-np.pi/2,
circ_radius=al/3, normal=(0, 1, 0),
body_radius=0.05, head_length=0.5)

arrow_phi = pv.Arrow(start=(0, 0, 0), direction=(al/2, al/2, 0),
shaft_radius=sr/2, tip_radius=0.02, tip_length=0.1, scale='auto')
arc_phi = semi_circular_arrow(center=(0, 0, 0), circ_frac=38/360, start_angle=0,
circ_radius=al/3, normal=(0, 0, 1),
body_radius=0.05, head_length=0.5)

# axes labels
axes_label_pts = [[al, 0., 0.], [0., al, 0.], [0., 0., al]]
axes_label_txt = ['x', 'y', 'z']

# angle labels
angles_label_pts = [[al/3*np.tan(22.5/180*np.pi), 0., al/3],
[al/3, al/3*np.tan(22.5/180*np.pi), 0.]]
angles_label_txt = ['θ', 'φ']

# the fish model came from: https://3dmag.org/en/market/download/item/6255/
reader = pv.get_reader('../resources/herring.stl')
fish = reader.read()

# centre the fish on the origin
b = fish.bounds
offset = (-((b[1]-b[0])/2 + b[0]), (b[3]-b[2])/2 + b[2], (b[5]-b[4])/2 + b[4])
fish.translate(offset, inplace=True)

# rotate the fish to fit our coordinate system
fish.rotate_x(-90, inplace=True)
fish.rotate_z(-90, inplace=True)

# scale the fish to length (along the z axis)
length = fish.bounds[5] - fish.bounds[4]
fish.scale(0.7*al/length, inplace=True)

# Assemble the 3D scene
p = pv.Plotter(window_size=[1600, 1200]) # to get a good size for raster outputs
p.add_mesh(fish, opacity=.9)
p.add_mesh(arrow_x, color='gray')
p.add_mesh(arrow_y, color='gray')
p.add_mesh(arrow_z, color='gray')
p.add_mesh(arrow_theta, color='red')
p.add_mesh(arc_theta, color='red')
p.add_mesh(arrow_phi, color='green')
p.add_mesh(arc_phi, color='green')

# the angle labels
p.add_point_labels(angles_label_pts[0], [angles_label_txt[0]], font_family='times',
bold=False, shape=None, always_visible=True, show_points=False,
font_size=50)
p.add_point_labels(angles_label_pts[1], [angles_label_txt[1]], font_family='times',
bold=False, shape=None, always_visible=True, show_points=False,
justification_horizontal='right', font_size=50)
# the axes labels
p.add_point_labels(axes_label_pts, axes_label_txt,
font_size=50, italic=True, bold=False, shape=None,
always_visible=True, show_points=False, font_family='times')

p.camera_position = [(10.0, 22.8, 15.4),
(4.7, 3.9, 3.3),
(0.97, -0.14, -0.19)]

# Unfortunately, all exports have issues...
# p.export_html('coordinate_system2.html') # loses all text
# p.export_obj('coordinate_system2.obj') # no colour?
# p.export_gltf('coordinate_system2.gltf') # loses text
# p.export_vrml('coordinate_system2.vrml') # no text??

# no greek symbols and scale > 1 loses some of the text
# p.screenshot('coordinate_system2.png', transparent_background=True, scale=1)

# Best option for the moment. Html would be preferrable if the text showed up
p.save_graphic(resourcesDir/'coordinate_system.svg') # raster, otherwise good

# this generates an on-screen version. But doesn't show the greek symbols
p.show()
7 changes: 5 additions & 2 deletions mkdocs.yml
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@ nav:

theme:
name: material
logo: echoSMs_logo.png
favicon: echoSMs_logo.png
logo: resources/echoSMs_logo.png
favicon: resources/echoSMs_logo.png
icon:
repo: fontawesome/brands/github
features:
Expand All @@ -21,6 +21,9 @@ theme:
repo_url: https://github.com/ices-tools-dev/echoSMs
edit_uri: edit/main/docs/

markdown_extensions:
- attr_list

plugins:
- search
- mkdocstrings:
Expand Down
38 changes: 10 additions & 28 deletions src/echosms/benchmarkdata.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,14 @@ class BenchmarkData:
Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015.
Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research.
Journal of the Acoustical Society of America 138, 3742–3764. <https://doi.org/10.1121/1.4937607>
Attributes
----------
angle_dataset : Pandas DataFrame
The angle dataset from the benchmark model runs.
freq_dataset : Pandas DataFrame
The frequency dataset from the benchmark model runs.
"""

def __init__(self):
Expand All @@ -22,31 +30,5 @@ def __init__(self):
angle_data_file = data_directory/'Benchmark_Angle_TS.csv'
freq_data_file = data_directory/'Benchmark_Frequency_TS.csv'

self.angle_data = pd.read_csv(angle_data_file)
self.freq_data = pd.read_csv(freq_data_file)

def dataset_angle(self):
"""Provide the angle benchmark data.
Returns
-------
: Pandas DataFrame
The TS as a function of incidence angle.
"""
return self.angle_data

def dataset_freq(self):
"""Provide the frequency benchmark data.
Returns
-------
: Pandas DataFrame
The TS as a function of frequency.
"""
return self.freq_data

# def dataset_angle(self, names):
# pass

# def dataset_freq(self, names):
# pass
self.angle_dataset = pd.read_csv(angle_data_file)
self.freq_dataset = pd.read_csv(freq_data_file)
12 changes: 6 additions & 6 deletions src/echosms/referencemodels.py
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ def models(self):
"""
return self.defs

def model_names(self):
def names(self):
"""Names of all model definitions.
Returns
Expand All @@ -46,8 +46,8 @@ def model_names(self):
"""
return [n['name'] for n in self.defs['target']]

def get_model_specification(self, name):
"""Get the model defintions for a particular model.
def specification(self, name):
"""Model defintions for a particular model.
Parameters
----------
Expand All @@ -67,8 +67,8 @@ def get_model_specification(self, name):
else:
return None

def get_model_parameters(self, name):
"""Get the model parameters for a particular model.
def parameters(self, name):
"""Model parameters for a particular model.
Model parameters are a subset of the model specification where the non-numerical
items have been removed.
Expand All @@ -84,7 +84,7 @@ def get_model_parameters(self, name):
The model parameters for the requested model or ``None`` if no model with that name.
"""
s = self.get_model_specification(name)
s = self.specification(name)

if s is None:
return s
Expand Down
27 changes: 14 additions & 13 deletions src/example_code.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# %%
"""Examples of using the echoSMs code to estimate scatter from objects."""

import matplotlib.pyplot as plt
Expand All @@ -11,13 +12,13 @@
# Load the reference model defintiions
rm = ReferenceModels()
print('Available reference models are:')
for n in rm.model_names():
for n in rm.names():
print('\t' + n)

# Load the benchmark data (from Jech et al., 2015)
bm = BenchmarkData()
bmf = bm.dataset_freq()
bm_theta = bm.dataset_angle()
bmf = bm.freq_dataset
bmt = bm.angle_dataset

# %% ###############################################################################################
# Run the benchmark models and compare to the frequency-varying benchmark results.
Expand Down Expand Up @@ -58,11 +59,11 @@

for name in names:
# Get the model parameters used in Jech et al. (2015) for a particular model.
s = rm.get_model_specification(name[0])
m = rm.get_model_parameters(name[0]) # the subset of s with string items removed
s = rm.specification(name[0])
m = rm.parameters(name[0]) # the subset of s with string items removed

# Add frequencies and angle to the model parameters
m['f'] = bmf['Frequency_kHz']*1e3 # [Hz]
m['f'] = bm.freq_dataset['Frequency_kHz']*1e3 # [Hz]
m['theta'] = 90.0

# and run these
Expand Down Expand Up @@ -106,27 +107,27 @@

for name in names:
# Get the model parameters used in Jech et al. (2015) for a particular model.
s = rm.get_model_specification(name[0])
m = rm.get_model_parameters(name[0]) # the subset of s with string items removed
s = rm.specification(name[0])
m = rm.parameters(name[0]) # the subset of s with string items removed

# Add frequencies and angle to the model parameters
m['f'] = 38000 # [Hz]
m['theta'] = bm_theta['Angle_deg']
m['theta'] = bmt['Angle_deg']

# and run these
ts = mod.calculate_ts(m, model_type=s['model_type'])

jech_index = np.mean(np.abs(ts - bm_theta[name[1]]))
jech_index = np.mean(np.abs(ts - bmt[name[1]]))

# Plot the mss model and benchmark results
fig, axs = plt.subplots(2, 1, sharex=True)
axs[0].plot(m['theta'], ts, label='echoSMs')
axs[0].plot(bm_theta['Angle_deg'], bm_theta[name[1]], label='Benchmark')
axs[0].plot(bmt['Angle_deg'], bmt[name[1]], label='Benchmark')
axs[0].set_ylabel('TS re 1 m$^2$ [dB]')
axs[0].legend(frameon=False, fontsize=6)

# Plot difference between benchmark values and newly calculated mss model values
axs[1].plot(m['theta'], ts-bm_theta[name[1]], color='black')
axs[1].plot(m['theta'], ts-bmt[name[1]], color='black')
axs[1].set_xlabel('Angle (°)')
axs[1].set_ylabel(r'$\Delta$ TS [dB]')
axs[1].annotate(f'{jech_index:.2f} dB', (0.05, 0.80), xycoords='axes fraction',
Expand All @@ -140,7 +141,7 @@
mss = MSSModel()

# Can add more variations to the model parameters
m = rm.get_model_parameters('weakly scattering sphere')
m = rm.parameters('weakly scattering sphere')
m['f'] = np.linspace(12, 200, num=800) * 1e3 # [Hz]
m['target_rho'] = np.arange(1020, 1030, 1) # [kg/m^3]
m['theta'] = [0, 90.0, 180.0]
Expand Down

0 comments on commit bb5f2ab

Please sign in to comment.