How do I calculate the full width at half maximum (FWHM) volume data in acoustic ultrasound simulation?

[Hank_H]

Hi

I am using SEFT for ultrasound simulation.
In the process of 'Analysis' results. I know how to view and get the result of maximum intensity from 2D slice view.
Could anyone point out how to calculate the volume data of full width at half maximum?

[montanaro]

The easiest way to do it (that I know of) is rather awkward. You export the whole field results to a vtk file and then use vtk to calculate the isosurface / isocontour of the Full Width at Half Maximum, then import it back into the modeling of Sim4Life and go to Tools ->Geometry to get the data you need (Surface Area, Volume, Bounding box size, etc). I'll attach a Python script below to: Export Field Data to a vtk file and use python vtk to extract the Full Width at Half Maximum, and to import it back into S4L. Then, in the modelling tab you can select the object and go to Tools->Geometry to get the data you want.

Notes:

Make sure to update 'simulation_name' and that the correct field is selected ("Overall Field" and "p(x,y,z,f)" in this case)

Select appropriate temporary file name and path in "vtk_exporter.FileName"

Do appropriate postprocessing of data (in my case, pressure field is complex data, I took the absolute value of this data in this line:
pabs = np.abs(pdata_wrapper.PointData[0] + 1j * pdata_wrapper.PointData[1])

Select appropriate threshold ... Keep this if you want FWHM
contour.SetValue(0, pabs_max / 2.0) # Half max

Careful that data when imported back into s4l is scaled correctly (conversion from meters to mm)
tr.Scale(1e3, 1e3, 1e3)

import s4l_v1 as s4l
import s4l_v1.document as document
import s4l_v1.analysis as analysis
import s4l_v1.model as model
import s4l_v1.units as units
import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs
import numpy as np

simulation_name = "Ac"

# Select Simulation
simulation = document.AllSimulations[simulation_name]
simulation_extractor = simulation.Results()
acoustic_sensor_extractor = simulation_extractor["Overall Field"]
# document.AllAlgorithms.Add(acoustic_sensor_extractor)

inputs = [acoustic_sensor_extractor.Outputs["p(x,y,z,f)"]]

vtk_exporter = analysis.exporters.VtkExporter(inputs=inputs)
vtk_exporter.FileName = u"C:\\Users\\montanaro\\Desktop\\ExportedData1.vtk"
vtk_exporter.Update()

# Read in vtk file containing field
reader = vtk.vtkDataSetReader()
reader.SetFileName(vtk_exporter.FileName)

# Convert cell centered data to node centered data 
# Needed for contour (isosurface) filter
cell2point = vtk.vtkCellDataToPointData()
cell2point.SetInputConnection(reader.GetOutputPort())
cell2point.Update()
# Turn into numpy array
pdata = cell2point.GetOutput()
pdata_wrapper = dsa.WrapDataObject(pdata)
# Get absolute value of real and imaginary pressure
pabs = np.abs(pdata_wrapper.PointData[0] + 1j * pdata_wrapper.PointData[1])
# Create new field with this calculated data
pdata_wrapper.PointData.append(pabs, 'pabs')
pdata.GetPointData().SetScalars(pdata.GetPointData().GetArray('pabs'))

# Create Contour Filter used to calculate isosurface
contour = vtk.vtkContourFilter()
contour.SetInputData(pdata)
contour.SetValue(0, pabs_max / 2.0) # Half max 

# Make sure we use correct scale between export and import
tr = vtk.vtkTransform()
tr.Scale(1e3, 1e3, 1e3)
transform = vtk.vtkTransformPolyDataFilter()
transform.SetInputConnection(contour.GetOutputPort())
transform.SetTransform(tr)

# Write vtk file containing contour filter as an stl file
writer = vtk.vtkSTLWriter()
writer.SetInputConnection(transform.GetOutputPort())
writer.SetFileName(vtk_exporter.FileName.replace(".vtk", ".stl"))
writer.SetFileTypeToBinary()
writer.Update()

# Import back into S4L
ents = model.Import(writer.GetFileName())

[Hank_H]

@montanaro

Hi montanaro

Thanks for all your help.
I got your point in your above script and it works on python.