# kWaveDiffusion

Time-domain simulation of heat diffusion and perfusion.

## Syntax

```kdiff = kWaveDiffusion(kgrid, medium, source)
kdiff = kWaveDiffusion(kgrid, medium, source, sensor)
kdiff = kWaveDiffusion(kgrid, medium, source, sensor, ...)
kdiff = kWaveDiffusion(kgrid, medium, source, [], ...)
```

## Description

`kWaveDiffusion` is a class definition for the time-domain solution of the diffusion equation or Pennes' bioheat equation in 1D, 2D, and 3D. In addition to heat diffusion, Pennes' bioheat equation accounts for advective heat loss due to tissue perfusion (the flow of blood through tissue), and heat deposition (e.g., due to ultrasound absorption) . The computation is based on a k-space pseudospectral scheme in which spatial gradients are calculated using the Fourier collocation spectral method, and temporal gradients are calculated using a k-space corrected finite difference scheme. For a homogeneous medium, the formulation is exact and unconditionally stable. For a heterogeneous medium, the time scheme allows larger time-steps to be taken for the same level of accuracy compared to conventional pseudospectral time-domain methods.

The most general partial differential equation solved by the `kWaveDiffusion` class is given by

```A * dT/dt = div(Kt * grad(T)) - B * (T - Ta) + Q
```

where

```A  = density [kg/m^3] * specific heat capacity [J/(kg.K)]
Kt = thermal conductivity [W/(m.K)]
B  = blood density [kg/m^3] * blood specific heat [J/(kg.K)] * blood perfusion rate [1/s]
Ta = arterial temperature (blood ambient temperature) [degC]
Q  = volume rate of heat deposition [W/m^3]
```

In a homogeneous medium, the thermal coefficients are related to the diffusion coefficient (or thermal diffusivity) by

```diffusion coefficient [m^2/s] = Kt / A
```

Note, when the diffusion coefficient is specified instead of the individual thermal coefficients, the equation that is solved is

```dT/dt = div(D * grad(T))
```

For non-constant coefficients, this differs from the conventional heat equation (where the diffusion coefficient is taken outside the divergence operator). For convenience, the thermal coefficients related to perfusion can also be combined to give a single "perfusion coefficient" given by

```perfusion coefficient [1/s] = B / A
```

The input parameters are assigned as fields to four input structures (`kgrid`, `medium`, `source`, and `sensor`) in the same way as the other models in the k-Wave toolbox. These structures define the properties of the computational grid, the distribution of medium properties, source terms, and the locations of the sensor points used to record the evolution of the temperature field over time.

The medium parameters can each be specified as a single scalar values in SI units (for homogeneous coefficients), or as matrices the same size as the computational grid (for heterogeneous coefficients).

The initial temperature distribution is specified by assigning a single scalar value or a matrix (the same size as the computational grid) to `source.T0`. A heat source can also be specified in the same way by defining `source.Q` (the volume rate of heat deposition).

The time history of the temperature field can be recorded automatically by specifying a series of sensor locations using sensor.mask. This is defined as a binary matrix (i.e., a matrix of 1's and 0's with the same dimensions as the computational grid) representing the grid points within the computational grid that will record the temperature field at each time step. The current sensor data can be queried at any point using the property `kdiff.sensor_data` (where `kdiff` is the name of the `kWaveDiffusion` object). The `sensor_data` is returned using MATLAB's standard column-wise linear matrix index ordering, indexed as `sensor_data(sensor_point_index, time_index)`. If no time dependent output is required, the sensor input can be replaced with an empty array `[]`.

After an object of the `kWaveDiffusion` class is created, the simulation is run by calling the method `kdiff.takeTimeStep(Nt, dt)`, where `kdiff` is the object name, `Nt` is the number of time steps to take, and `dt` is the size of the time step. During the simulation, a visualisation of the temperature field is displayed. The current temperature can be queried (or modified) at any point using the property `kdiff.T`.

 Pennes, H. H. (1948). Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of Applied Physiology, 1(2), 93-122.

## Inputs

 `kgrid` grid object returned by `kWaveGrid` `medium.diffusion_coeff` diffusion coefficient [m^2/s] ` ` OR `medium.density` tissue mass density [kg/m^3] `medium.specific_heat` tissue specific heat capacity [J/(kg.K)] `medium.thermal_conductivity` tissue thermal conductivity [W/(m.K)] `medium.perfusion_coeff` perfusion coefficient [1/s] ` ` OR `medium.blood_density` blood mass density [kg/m^3] `medium.blood_specific_heat` blood specific heat capacity [J/(kg.K)] `medium.blood_perfusion_rate` blood perfusion rate [1/s] (volumetric flow rate per unit volume) `medium.blood_ambient_temperature` ambient blood temperature within perfused tissue regions [degC] `medium.diffusion_coeff_ref` reference diffusion coefficient used within the k-space operator (default = 'max') `medium.perfusion_coeff_ref` reference perfusion coefficient used within the k-space operator (default = 'max') `source.T0` initial temperature distribution [degC] `source.Q` volume rate of heat deposition [W/m^3](note, `medium.density` and `medium.specific_heat` must be defined when using `source.Q`) `sensor.mask` binary grid specifying where the temperature is recorded at each time step

## Optional Inputs

Optional 'string', value pairs that may be used to modify the default computational settings.

Input Valid Settings Default Description
`'DisplayUpdates'` (Boolean scalar) `true` Boolean controlling whether details of the simulation are printed to the MATLAB command line.
`'MovieArgs'` (string cell array) `{}` Settings for `VideoWriter`. Parameters must be given as {'param', value, ...} pairs within a cell array (default = {}), where 'param' corresponds to a writable property of a VideoWriter object.
`'MovieName'` (string) `'date-time-kWaveDiffusion'` Name of the movie produced when `'RecordMovie'` is set to `true`.
`'MovieProfile'` (string) `'Uncompressed AVI'` Profile input passed to `VideoWriter`.
`'PlotFreq'` (integer numeric scalar) `10` The number of iterations which must pass before the simulation plot is updated.
`'PlotScale'` (numeric two element vector) or
`'auto'`
`'auto'` [min, max] values used to control the scaling for `imagesc` (visualisation). If set to `'auto'`, a symmetric plot scale is chosen automatically for each plot frame.
`'PlotSim'` (Boolean scalar) `true` Boolean controlling whether the simulation iterations are progressively plotted.
`'RecordMovie'` (Boolean scalar) `false` Boolean controlling whether the displayed image frames are captured and stored as a movie using `VideoWriter`.

## Outputs

 `kdiff` `kWaveDiffusion` object which can be used to run thermal simulations using the diffusion equation or Pennes bioheat equation

## Dynamic Properties

Properties which can be queried or modified after the object is created.

Fieldname Description
`cem43` thermal dose given in cumulative equivalent minutes (cem) relative to T = 43 degC [mins]
`T` current temperature field [degC]
`Q` volume rate of heat deposition [W/m^3]

## Static Properties

Properties which can be queried, but not modified, after the object is created.

Fieldname Description
`dt_limit` maximum time step for which the simulation remains stable [s]
`lesion_map` binary matrix of `cem43 >= 240` mins
`lesion_size` total size of `lesion_map` (distance in 1D [m], area in 2D [m^2], volume in 3D [m^3])
`sensor_data` time varying temperature recorded at the sensor positions given by `sensor.mask` [degC]

## Methods

Methodname Description
`plotTemp` plot current temperature field in current figure window
`setOptionalInputs('string', value, ...)` modify the optional inputs after the object is created
`takeTimeStep(Nt, dt)` calculate the given number of time steps of the temperature field

`bioheatExact`