waltsims / k-wave-python

A Python interface to k-Wave GPU accelerated binaries
https://k-wave-python.readthedocs.io/en/latest/
GNU General Public License v3.0
102 stars 30 forks source link

[Help Wanted] #459

Closed TullyMonster closed 3 weeks ago

TullyMonster commented 3 weeks ago

Description I have collected a dataset of vascular images along with vascular semantic segmentation. Previously, based on the manually segmented images, I completed binarization of the vessels, marking vessel pixels as 1 and non-vessel pixels as 0, represented by the variable vessel_matrix.

Next, I plan to simulate ultrasound signals based on vessel_matrix.

1. Import:

from pathlib import Path

import numpy as np
from kwave.data import Vector
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksensor import kSensor
from kwave.ksource import kSource
from kwave.kspaceFirstOrder2D import kspaceFirstOrder2DC
from kwave.options.simulation_execution_options import SimulationExecutionOptions
from kwave.options.simulation_options import SimulationOptions
from kwave.utils.mapgen import make_cart_circle

2. kgrid, medium and sensor:

kgrid = kWaveGrid(
    N=Vector(list(vessel_matrix.shape)),  # Number of grid points
    spacing=Vector([1e-4, 1e-4])  # Physical size of each grid point
)
medium = kWaveMedium(sound_speed=1540,  # Sound speed in water or tissue
                     density=1000,
                     alpha_coeff=0.75,  # Absorption coefficient
                     alpha_power=1.5,  # Power index of absorption coefficient
                     BonA=6)  # Nonlinear absorption parameter in the medium

kgrid.makeTime(c=medium.sound_speed,
               t_end=(kgrid.spacing * kgrid.N * 2.2 / medium.sound_speed).max())  # 2.2 times the spatial length to ensure enough reflection paths

sensor = kSensor(mask=make_cart_circle(radius=1e-5, num_points=20))

3. source:

# Segmentation mask of the vessels. Each pixel is marked as 1 if it belongs to a vessel, and 0 otherwise
vessel_matrix: np.ndarray = ...  # shape: (584, 565)
source = kSource()
source.p0 = 1e6 * vessel_matrix

4. simulate:

def simulate(**opts):
    simul_opts = SimulationOptions(save_to_disk=True, **opts)
    exec_opts = SimulationExecutionOptions(is_gpu_simulation=False)

    return kspaceFirstOrder2DC(kgrid=kgrid, source=source, sensor=sensor, medium=medium,
                               simulation_options=simul_opts,
                               execution_options=exec_opts)

result = simulate(pml_inside=False, data_cast='single')

5. result:

Nt = result['Nt'].item()
p: np.ndarray = result['p']

print(f'{Nt = }')
print(f'{p.shape = }')
Nt = 4283
p.shape = (4283,)

I thought, (Nx, Ny, Nt) was my target shape, but the shape I saw was (Nt,).

So, I guess this is just the result of getting the change time of the first or last sensor unit.

What I want to ask is, how to traverse each unit?

djps commented 3 weeks ago

I think you need to check the sensor information: make_cart_circle has a subvoxel radius. When I try and use cart2grid it gives a single point at the centre of the grid, which I assume is what is being recorded.

Also, it is probably best to continue this on the discord channel.

TullyMonster commented 3 weeks ago

Thank you very much! I will continue to address this on discord, thanks for the inspiration!

waltsims commented 3 weeks ago

Cross-referencing response on discord here.