WO2022132079A2 - A computer applied simulation system for acoustic/electromagnetic imaging and/or tomography - Google Patents

Abstract

This invention relates to a computer applied simulation method (100) enabling realistic wave propagation simulations particularly in heterogeneous media in all angular sectors for acoustic/electromagnetic imaging and/or tomography by use of low computer sources.

Description

A COMPUTER APPLIED SIMULATION SYSTEM FOR ACOUSTIC/ELECTROMAGNETIC IMAGING AND/OR TOMOGRAPHY

Technical Field of the Invention

This invention relates to a computer applied simulation method enabling realistic wave propagation simulations particularly in heterogeneous media in all angular sectors for acoustic/electromagnetic imaging and/or tomography by use of low computer sources.

Background of the Invention

Ultrasound tomography (UST) method was firstly proposed in 1977 for imaging tissue ultrasound parameters (J. F. Greenleaf, S. A. Johnson, and R. C. Bahn, ’’Quantitative Cross-Sectional Imaging of Ultrasound Parameters,” 1977 Ultrasonics Symposium, Phoenix, AZ, USA, 1977, pp. 989-995).

The initial algorithms used for UST image reconstruction were adopted from X-ray tomography and based on the assumption that ultrasound waves propagate as rays (Duric N et al 2013 Breast imaging with the softvue imaging system: first results, SPIE Medical Imaging Int. Society for Optics and Photonics p 86750K). Although fast imaging is possible by use of this assumption, since said assumption is not true for ultrasound waves used in clinical frequency range, image quality is degraded. As a more realistic model, a filtered back propagation method has been developed wherein wave diffraction effects are included (A J Devaney. A filtered backpropagation algorithm for diffraction tomography. Ultrasonic Imaging, 4:336-350, 1982.).

With the developments in computational power, more realistic wave propagation models can be resolved in a relatively shorter period. As disclosed under the US patent numbered US2015/0025388 in the related art, waveform tomography provides a high resolution UST image by taking the propagation effects such as multiple dispersion and diffraction into account. In waveform tomography, imaging scenario is simulated with a computer model, and simulation results are compared with the signals obtained from the transducers (experimental data) (Tarantola, A. (1987). Inverse Problem Theory: Methods for Data Fitting and Parameter Estimation (Elsevier, Amsterdam), pp. 1-601 ). As also disclosed under another US patent document numbered US2016/0030000, sound and/or attenuation coefficient in the field of view is updated according to difference between simulations and the experimental data. The simulations are repeated until difference between experimental data and the data from the computer simulations is smaller than a predetermined level. In order to apply this method, highly accurate simulation means to solve the wave propagation equation in heterogonous media in a practically applicable duration are needed.

Since image generation process depends on simulation results, simulation time should be minimized. Realistic simulation models use full-wave simulation methods such as finite differences (Perez-Liva, M., Herraiz, J. L., Udias, J. M., Miller, E., Cox, B. T., Treeby, B. E. (2017). Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion. The Journal of the Acoustical Society of America, 141 (3),1595-1604). Since the computational power and memory needed by these methods are high, two-dimensional models are used in practice in order to minimize the simulation time. Therefore, not only the image resolution is low along one axis due to two- dimensional tissue assumption, but also scattering in three dimensions can not be modelled.

In order to solve this problem, a 3-dimensional UST method providing fast solution for speed of sound reconstruction, in poor heterogonous media was proposed based on travel-time (Martiartu, Naiara Korta, Christian Boehm, and Andreas Fichtner. ”3d wave-equation-based finite-frequency tomography for ultrasound computed tomography.” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control (2020)). Although this method provides a three-dimensional solution for the reconstruction of sound speed, a reconstruction cannot be made for attenuation coefficient inside the medium. In addition, suggested model is valid for the case when tissue heterogeneity is low (where sound speed changes at the most +10%). In particular, speed of sound in bones can change by about 100% compared to soft tissues. For that reason, said method has limited usage.

In another US patent document numbered US8568318 in the related art, a method for angular spectrum based on Fourier transform in UST simulations is disclosed. In the method disclosed under said document, low heterogeneity assumption is also required, which limits the applicability of this method for imaging soft tissues. In addition, as multiple reflections are not taken into account, image quality can not be improved adequately.

Microwave tomography is also a method providing imaging electromagnetic features of tissue in a similar way (conductivity, complex dielectric coefficient). The problems encountered in UST are also seen in microwave tomography. Vyas et al. proposed an angular spectrum method for fast ultrasound simulations, although not related to UST ( Vyas, Urvi, and Douglas Christensen. ’’Ultrasound beam simulationsin inhomogeneous tissue geometries using the hybrid angular spectrum method.” IEEE transactions on ultrasonics, ferroelectrics, and frequency control 59, no. 6 (2012): 1093- 1100). Even though the method suggested by Vyas et al. is applicable in highly heterogeneous tissues, it is not convenient for tomographic imaging as it can provide solution in a limited angular sector due to plane wave assumption. This method was proposed for the simulation of focused ultrasound beams, but not for imaging.

Below, the methods applied in the related art for wave propagation from acoustic sources are described. Similar methods can also be used for electromagnetic waves.

In the related art, acoustic pressure field in the linear medium is obtained by solution of the non-homogenous Helmholtz equation in frequency domain.

Here f(x, y, z) is acoustic source term. Acoustic pressure field is defined as the local pressure difference applied by the acoustic waves with respect to the initial pressure pf the medium.

Pressure field in a plane located in z position can be divided into plane wave components as follows:

Pressure field in another Zd plane in homogenous medium with wave number k can be calculated as follows:

Vector propagation coefficient k is sum of propagation coefficients in x, y, and z directions:

As equation (2) can be obtained by the Fast Fourier Transformation (FFT) method, plane wave propagation in homogenous medium can be calculated very fast using equation (3). However, equation (3) cannot be used for heterogeneous media.

In cases where medium parameters such as sound speed, object density, attenuation coefficient vary as a function of position, medium is defined as heterogonous. In case of acoustic and ultrasound imaging problems, acoustic waves are sent from the acoustic sources to the medium, and the waves reflected back from the medium and/or waves passed through medium are detected by acoustic sensors, and position dependent acoustic parameters are predicted.

In the related art, Hybrid Angular Spectrum (HAS) method was suggested for solution of pressure field in non-homogenous medium using plane wave assumption (Vyas, Urvi, and Douglas Christensen. ’’Ultrasound beam simulationsin inhomogeneous tissue geometries using the hybrid angular spectrum method.” IEEE transactions on ultrasonics, ferroelectrics, and frequency control 59, no. 6 (2012): 1093-1100). In HAS method, simulation region is divided into two or three dimensional rectangular cells depending on problem type, and acoustic features are assigned to each cell. In this method acoustic source must be outside simulation region. Firstly, acoustic pressure distribution on the simulation boundary (surface) that is nearest to the acoustic source is calculated. Then, plane wave spectrum of this pressure distribution is calculated by using equation (2). Pressure field in the next plane of the simulation domain is calculated in two steps. Firstly, average wave propagation coefficient in the area between the planes is calculated as follows:

In equation (6), sound speed in (i, j, k) cell in Ci , _k simulation region is / ultrasound wave frequency.

In the first step, the spatial pressure field, the difference between propagation constant of each cell in the plane and the average propagation constant, and the reflection coefficient are used to calculate the interim pressure of cells in the next plane:

Is the reflection coefficient, which is calculated as follows:

In equation (8), pij,_k and aij,_k are density and attenuation coefficient in cell (i, j, k), respectively. 6ki,j,_k is the deviation from the average wave number:

z’ is the distance of the plane wave fronts between grids (i + 1 , j) and (i, j). 0_d, is plane wave propagation direction, z” is total distance between successive cells, which is calculated as d_z / cos (0d). .

In the second step, the transmitted pressure field is converted into spectral domain by Fourier transform, and transmitted to the next plane in the direction of maximum amplitude using the average wave number:

Then, the pressure field is converted back to the spatial domain by inverse Fourier transform. This process is repeated until the pressure field is reached to the last plane in the simulation geometry. To calculate the back-scattered waves, reflected pressure distribution in each plane is recorded. Simulation is proceeded towards initial plane starting from the last plane by means adding the saved reflected pressure field in each plane.

HAS method is convenient for cases where plane wave acoustic pressure field propagates inside heterogeneous media, but not convenient for point sources. In this method, simulation grid should be arranged in a manner perpendicular to the wave propagation direction. For this reason, in the case when there is more than one main propagation direction, or when transducer (acoustic source) axis changes, grid should be rearranged. In addition, this argument is also true for back scattered waves. Scattering in the direction perpendicular to the propagation vector cannot be taken into account. Consequently, HAS method is not a suitable for imaging problems and also, it does not support the use of acoustic sources inside the simulation region.

Brief Description of the Invention

Purpose of the invention is to provide to a computer applied simulation method enabling realistic wave propagation simulations particularly in heterogeneous media in all angular sectors for acoustic/electromagnetic imaging and/or tomography by use of low computer sources.

In order to achieve the purpose of the invention, computer applied simulation method described under claim one and claims dependant on this claim and particularly enabling modelling of wave propagation in heterogonous medium, discloses a simulation region where simulation process is made and a simulation region bounded by a first plane from at least one end and a final plane at the other end in a manner to have planes more than one between them. Waves emanating from the acoustic/electromagnetic sources are decomposed according to their propagation direction by preferably taking spatial Fourier transform. Simulation region is divided into cubical cells formed by means of preferably a first plane, a last plane, and interim planes. Acoustic/electromagnetic field propagated by each source plane wave component is proceeded toward the last plane starting from the initial plane of the simulation region. Back scattered and transmitted wave components due to non-heterogonous feature of medium are calculated during propagation. Back scattered component in each cell is recorded for the related plane. When the wave propagates from one plane to another, firstly pressure weighted average propagation constant of all planes is calculated. In this way cells having high pressure contributes to general propagation coefficient more. Propagation directions of the waves with an amplitude above a predetermined threshold calculated. Relative propagation constant is calculated by subtracting the absolute propagation constant of each cell in the plane from average propagation constant. Pressure progress from each cell in a plane to the next plane, is calculated using the phase and amplitude, relative propagation constant and wave propagation direction. Plane wave spectrum is calculated preferably by Fourier transform of the pressure distribution in the plane. Average propagation constant of the plane wave spectrum is used to propagate the wave to next plane. Then, spatial pressure distribution is calculated by taking Inverse Fourier Transform. At this stage, back scattered and transmitted pressure is also calculated and thus pressure in next cell is found. This operation is conducted for all propagation directions. Then, pressure in the next plane is calculated in a similar way. Same operations repetitively proceeded until the last plane of simulation region. Fourier transform of the back scattered components are taken and decomposed into plane wave directions as done for the waves emanating from acoustic/electromagnetic sources. These compositions are used for wave propagation simulation in the respective direction. Then same operations are applied in the reverse direction. Components of the acoustic/electromagnetic sources propagating in reverse direction and back scattered waves calculated in the initial iteration are summed to provide the source for simulating the wave propagation from the last plane to the first plane. Wave propagation in the back and forward directions along the same axis is repeated and all reflections are included in the calculation. Then wave propagation in other axes are made in a similar way and pressure distribution in all directions are collected and total pressure distribution is obtained. This operation is repeated once or several times depending on the back scattered acoustic/electromagnetic density.

Detailed Description of the Invention

Computer applied simulation system for acoustic/electromagnetic imaging and/or tomography to achieve purpose of this invention is shown in the drawings and the drawings are as follows;

Figure 1 shows a flow diagram of an application of the simulation method of the invention.

Figure 2 shows a schematic view of the simulation region that is separated into planes in an application of the method of the invention for wave propagation in the ± x directions in Cartesian coordinate system.

Figure 3 shows a schematic view of the simulation region that is separated into planes in an application of the method of the invention for wave propagation in the ± z directions in Cartesian coordinate system.

Figure 4 shows a schematic view of the simulation region that is separated into planes in an application of the method of the invention for wave propagation in the ± y directions in Cartesian coordinate system.

Figure 5 shows a schematic view of the spectral filter applied to separate the pressure components propagating mainly in the + x direction in a three-dimensional frequency space in an embodiment the invention. Figure 6 shows a schematic view of the spectral filter applied separate the pressure components propagating mainly in the + x, - x, + y and - y directions in a two-dimensional frequency space in an embodiment of the invention.

Figure 7 shows comparison of resultant pressure field distribution obtained using the method of the invention and methods in the related art in an illustrative two-dimensional simulation case.

The parts indicated in the figures have been designated separate numbers and the said numbers are given below:

SR: Simulation Region

Pi: First plane limiting simulation region from one end

P2: Last plane limiting simulation region from the end not limited by the first plane

P_n: Planes between the first plane and last plane x: Ordinary positive x axis in Cartesian coordinate system y: Ordinary positive y axis in Cartesian coordinate system z: Ordinary positive z axis in Cartesian coordinate system

+x: Proceeding direction of wave on positive x axis

-x: Proceeding direction of wave on negative x axis

+y: Proceeding direction of wave on positive y axis

-y: Proceeding direction of wave on negative y axis

PS-FDTD: Pseudo spectral Finite Difference Time Domain method

HAS: Hybrid Angular Spectrum method

Particularly, computer applied simulation method (100) of the invention allowing simulation of wave propagation in heterogonous media comprises process steps of: definition of at least one simulation region (SR) extending between at least one first plane (Pi) and at least one last plane (P₂) having plurality of planes (P_n) extending between them;

• conversion of the field emanating from acoustic/electromagnetic sources inside the simulation region (SR) into spectral field distribution (102);

• separation of spectral field distribution into components extending in at least two directions on each axis of at least two axis making angle with one another other than 0 and 180 degrees (103);

• propagating each component by calculating the transmitted and reflected fields for each plane inside the simulation region (Pi, P_n, P₂) starting from the first plane (Pi) or from last plane (P₂) to the last plane (P₂) or to the first plane (Pi) (104);

• calculation of one primary total field distribution in the simulation region (SR) when the propagation process is ended, and calculating spectral field distribution of the reflected field to find reflected field propagation directions (105);

• propagating the field in the opposite direction of step (104) by adding reflected field components in the related direction for each plane (Pi, P_n, P₂) in the simulation region (SR);(106);

• calculating the updated total field distribution inside the simulation region (SR) when propagation ends (107);

• checking whether or not at least one rule defined in advance is verified in order to ensure the convergence of the simulation is above a desired level (108);

• returning to step 104 if the rule is not verified as a result of the checking operation (108);

• summing up the total field calculated for all components if the rule is verified as a result of checking process (108) (109).

In one embodiment, the predefined rule for the application of the computer applied simulation method (100) of the invention comprises checking if primary total pressure density and the updated total pressure density in the next iteration equals to or smaller than a predefined value in order to ensure the convergence of the simulation is above a desired level. In this application when checking process (108) is initiated after calculation of the updated total pressure density (107), the difference between primary total pressure density and updated total pressure density is calculated and if the calculated difference is equal to or smaller than the predetermined value, total pressure densities calculated for all components are summed up (10) and the simulation (100) is terminated. If the difference between primary total pressure density and updated total pressure density as a result of the checking process (108) is larger than the predetermined value, processes between steps 104 to 108 are repeated until the difference between primary total pressure density and updated total pressure density is equal to or smaller than the predetermined value.

In one embodiment, the predefined rule for the application of the computer applied simulation method (100) of the invention comprises process of checking if the ratio between the primary total pressure density and updated total pressure density equals to or smaller than a predefined value in order to ensure the convergence of the simulation is above a desired level. In this application when checking process (108) is initiated after calculation of update total pressure density (107), ratio of the primary total pressure density and the updated total pressure density is calculated and if the calculated ratio is equal to or smaller than the predetermined value, total pressure densities calculated for all components are summed up (10) and the simulation (100) is terminated. If the ratio between the primary total pressure density and the updated total pressure density as a result of checking process (108) is larger than predetermined value, processes between steps 104 to 108 are repeated until the ratio of the primary total pressure density and the updated total pressure density is equal to or smaller than the predetermined value.

In one embodiment, the predefined rule for the application of the computer applied simulation method (100) of the invention comprises process of checking whether or not processes in steps 104 to 108 are repeated in a predefined number of times to ensure the convergence of the simulation is above a desired level. In this application when checking process (108) is initiated after calculation of the updated total pressure density (107), the number of repetitions of the processes from step 104 to step 108 is calculated and if that calculated value is equal to the predetermined value, total pressure densities calculated for all components are summed up (109) and the simulation (100) is terminated. If the number of repetitions of the process from steps 104 to 108 is smaller than the predetermined value after checking process (108), processes from step 104 to step 108 are repeated until the number of repetitions of the processes in steps 104 to 108 is equal to predetermined value. In an embodiment of the invention, in the step of separation of spectral field distribution into directional components (103), the spectral field distribution is separated into at least four components extending in directions opposite to each other on each of at least two axes.

In a special embodiment of the invention, in the step of separation of spectral field distribution into directional components (103), spectral field distribution is separated into at least six components extending in directions opposite to each other on each of at least three axes.

In the computer applied simulation method (100) of the invention, in an embodiment of the process step of separation of spectral field distribution into components extending in at least two directions on each axis of at least two axis making angle with one another other than 0 and 180 degrees (103), spectral field distribution is separated into components in propagation directions of ordinary ± x and ± y in Cartesian coordinate system by preferably two dimensional Fast Fourier Transform method.

In the computer applied simulation method (100) of the invention, in an embodiment of the process step of separation of spectral field distribution into components extending in at least two directions on each axis of at least two axis making angle with one another other than 0 and 180 degrees (103), spectral field distribution is separated into components in propagation directions of ordinary ± x, ± y and ± z in Cartesian coordinate system by preferably three dimensional Fast Fourier Transform method.

An illustrative embodiment of the computer applied simulation method (100) of the invention is described below, wherein the acoustic sources are used and the spectral pressure field is separated into components mainly in propagation directions of ordinary ± x, ± y and ± z in a Cartesian coordinate system.

After defining the simulation region (SR) where it is desired to perform the simulation (101 ), pressure field distribution (P_s (x, y, z)) remaining inside simulation region (SR) generated by the acoustic sources is converted into spectral pressure field by preferably Fast Fourier Transformation (FFT) method.

Spectral pressure field obtained thereafter is separated into ordinary ± x, ± y and ± z direction components in the Cartesian coordinate system by a filtering operation performed in simulation region (SR) (103). In a preferred embodiment, each component is adapted in a manner to support wave propagation between ± 45 degrees around a main direction. While 4 components extending in ± x, ± y directions are used for a two-dimensional simulation application, 6 components extending in ± x, ± y and ± z directions are used for a three-dimensional application:

The spectral filter S in the equation numbered (12) takes value of one for the related main propagation direction region and zero for other direction regions.

In the described illustrative embodiment, component of spectral pressure region in +x direction is proceeded from an initial plane (Pi) of simulation region (SR) until last plane (P₂) and pressure propagation in each plane (Pi, P_n, P₂) for +x propagation of spectral pressure is calculated by use of following equation (104).

Here ‘ ^; represents propagation transfer operation defined in equations (7) and (10).

During propagation of the spectral pressure field from initial plane (Pi) to last plane (P₂) of the simulation region (SR) of component in +x direction, reflected pressure field (P_r ^{+ X} (x, y, z)) in each plane is also calculated (104).

When the last plane (P₂) boundary of the simulation region (SR) is reached, in other words, when propagation reaches the end, primary total pressure density in simulation region (SR) is calculated (105). When the last plane (P₂) boundary of the simulation region (SR) is reached, in other words, when propagation reaches the end, reflected pressure field (P_r ^{+ X} (x, y, z)) is separated into subcomponents (- x, ± y and ± z) in 5 different propagation directions by use of preferably Fast Fourier Transformation (FFT) method.

Then reflected pressure field propagating in the -x direction calculated in step 105 is added to the component of spectral pressure field propagating in the -x direction and pressure propagation from the last plane (P₂) to initial plane (Pi) of simulation region (SR) is calculated by use of the following equation (106):

During propagating spectral pressure field with added spectral field component in the -x direction and the reflected pressure field propagating in -x direction, from the last plane (P₂) to initial plane (Pi) of simulation region (SR), reflected pressure field (P_r ^{- X} (x, y, z)) at each plane (Pi, P_n, P₂) is also calculated (106).

When the initial plane (Pi) boundary of the simulation region (SR) is reached, in other words, when propagation reaches to the end, updated total pressure density in the simulation region (SR) is calculated (107).

Then in order to find out whether or not simulation is conducted at the desired convergence level, it is checked whether or not a predefined rule is executed. In an embodiment defined in a manner comprising checking whether or not the difference between the primary total pressure density and the updated total pressure density is equal to or smaller than a predetermined value, if it is discovered that the difference between primary total pressure density and updated total pressure density is larger than the predetermined value as a result of checking process (108), calculation of the total pressure densities in the + x and -x directions is continued repetitively in similar manner until the said two values are equal to one another, in other words, it is returned to step 104 and processes from step 104 to step 108 are repeated. If as a result of control process (108), the difference between primary total pressure and updated total pressure density is found to be equal or smaller than the predetermined value, above mentioned steps are applied to ±y and ±z direction components after above mentioned operations, or simultaneous with the above mentioned processes, and the total pressure densities calculated for each ±x, ±y and ±z direction components are summed up to find the resultant pressure density.

In Figure 7, results obtained from computer applied method (100) of the invention and methods in the related art in an illustrative two-dimension simulation operation are shown. Here an object having a sound speed 33 % higher than background plan is placed in a homogenous medium. It is seen that Pressure propagation calculated by full-wave Pseudo spectral finite difference time domain method (PS-FDTD) and proposed method give similar results. Said simulation process lasts for 19 seconds with the computer applied method of the invention and 370 seconds with PS-FDTD method. In addition, it is also seen that HAS method of prior art shown in said Figure 7 can not model reflections adequately.

The invention also relates to a data processing apparatus comprising instruments adapted to realize computer applied simulation method (100) described above.

The invention also relates to a computer program comprising instructions of computer for realization of the steps of computer applied simulation method (100) described above when applied by computer.

This invention also relates to a data carrier that can be readable by computer wherein above described computer program can be stored.

The invention also relates to method of generating an image wherein computer applied simulation method (100) described above is used. More specifically, said image generation method comprises process steps of transformation into frequency space of at least one signal incident into a medium from at least one acoustic/electromagnetic source and reflecting from the medium and/or passing through medium and received by at least a sensor, obtaining at least a signal received by at least a sensor by applying above described computer applied simulation method (100) by means of simulation method, updating material parameters in simulation region (SR) used by simulation method (100) subject to one or more than one features of signal by means of simulation if it is discovered that the difference between them is larger than a predetermined threshold value as a result of comparison of the signal received by at least a sensor to the signal obtained by means of a simulation, obtaining updated related signal by applying the computer applied simulation method (100) by use of material parameters updated in simulation region (SR) and repeating the steps described above until the difference between them is equal to or less than the threshold value determined before as a result of comparison of updated signal by simulation method to signal received by at least a sensor. Material parameters can be sound speed and/or attenuation coefficient and/or density of tissue or displayed object in acoustic and ultrasound imaging. For instance, since different tissues such as fat tissue, mammary glands, cysts, benign or malign tumours have different material parameters, they can be distinguished from one another by imaging such parameters. Thus, pathologic tissues can be diagnosed. In electromagnetic imaging, material parameters can be dielectric coefficient and/or conductivity and/or attenuation coefficient, and different tissues and materials can be distinguished by said parameters in a similar way as described above for acoustic and ultrasound imaging.

Instead of selecting main propagation direction in HAS method, computer applied method (100) of the invention propagates waves in more than one direction simultaneously. Thus, simulation of isotropic sources can also be conducted. In order to decrease simulation time, number of propagation directions can be restricted according to a certain threshold normalized to maximum direction component. Reflected pressure field (P_r ^{+ X} (x, y, z)) in related planes (Pi, P_n, P₂) of the defined simulation region (SR) can also be calculated and stored.

Computer applied simulation method (100) of the invention is a frequency based method and thus can be transformed to signal for broad band signals, to frequency domain by Fourier Transformation and resolved for each frequency component. The components can be parallelized by use of CPU and/or GPUs and thus fast solution can be achieved. In light of this basic concept, it is possible to develop various applications of a computer applied method (100) of the invention and the invention can not be limited to the examples described under the invention and is fundamentally as described under the claims.

Claims

CLAIMS A computer applied simulation method (100) allowing modelling of wave propagation in particularly heterogonous medium, characterized by comprising the process steps of

• definition (101 ) of at least one simulation region (SR) extending between at least one first plane (Pi) and at least one last plane (P₂) having plural number of planes (P_n) extending between them;

• transformation of the incident pressure field inside the simulation region (SR) into spectral pressure field (102);

• separation of spectral pressure field into components extending in at least two directions on each axis of at least two axes making angle with one another other than 0 and 180 degrees (103);

• propagating (104) each component by calculating transmitted and reflected pressure in each of the planes (Pi, PN, P2) in the simulation region starting from the first plane (Pi) to the last plane (P₂) or from last plane (P₂) or first plane (Pi) in the related direction of the simulation region (SR);

• calculation of one primary total pressure density inside the simulation region (SR) when propagation ends and detecting propagation directions for the reflecting pressure field (105);

• propagating the pressure field in the opposite direction of step 104 by adding the reflected (106) pressure field in the related direction onto each component and by calculating the pressure in each plane (Pi, PN, P2);

• calculation of the updated total pressure density inside simulation region (SR) when propagation reaches to a boundary plane (Pi or P₂) (107);

• checking whether or not at least one rule defined in advance is verified in order to ensure the convergence of the simulation is at a desired level (108);

• return to step 104 if rule is not verified as a result of the checking operation (108);

• summing (109) up total pressure densities calculated for all components if the rule is verified as a result of the checking process (108). The computer applied simulation method (100) according to claim 1 , characterized by comprising the process of checking the rule used in checking process (108), wherein the rule is that the difference between the primary total pressure density and the updated total pressure density is equal to or smaller than a predetermined value. The computer applied simulation method (100) according to claim 1 , characterized by comprising the process of checking the rule used in checking process (108), wherein the rule is that ratio of the primary total pressure density and the updated total pressure density is equal to or smaller than a predetermined value. The computer applied simulation method (100) according to claim 1 , characterized by comprising the process of checking the rule used in checking process (108), wherein the rule is that the processes of steps 104 to 108 are repeated predetermined number of times. The computer applied simulation method (100) according to any one of above claims, characterized by separating the spectral pressure distribution into at least four components wherein each of at least two axes extending in directions opposite to each other in the step of separation of spectral pressure region into components (103), The computer applied simulation method (100) according to claim 5, characterized by the spectral pressure distribution being separated into six components extending in directions opposite to each other on each of at least three axes. The computer applied simulation method (100) according to any one of above claims, characterized by spectral pressure field being separated into components of ordinary ± x and ± y propagation directions in the Cartesian coordinate system by two-dimensional Fast Fourier Transformation method. The computer applied simulation method (100) according to any one of claims 1 to 6, characterized by spectral pressure field being separated into components of ordinary ± x, ± y and ± z propagation directions in the Cartesian coordinate system by three- dimensional Fast Fourier Transformation method. A data processing apparatus comprising instruments adapted to achieve a computer applied simulation method (100) according to any one of above claims. A computer program comprising instructions providing realization of steps of a computer applied simulation method (100) according to any one of claims 1 to 8 where applied by a computer. A data carrier readable by a computer wherein computer program of claim 10 can be stored. An image generation method wherein a computer applied simulation method (100) according to claim 1 to claim 8 is used. The image generation method according to claim 12, characterized by comprising the process steps of

• transforming at least one signal incident into a medium from at least one acoustic/electromagnetic source which is reflected from the medium and/or passed through medium and received by at least one sensor into frequency domain,

• obtaining at least one signal received by at least one sensor by applying the computer applied simulation method (100) according to any of the claim 1 to claim 8,

• updating material parameters inside simulation region (SR) used in simulation method (100) according to one or more than one features of signal if the difference between the signal received by at least one sensor and the signal obtained by means of the simulation is larger than a predetermined threshold,

• obtaining updated signal by applying computer applied simulation method (100) wherein updated material parameters is used in the simulation region (SR) and, repeating the steps described above until the difference between the updated signal by simulation method and the signal received by at least a sensor is equal to or less than the predetermined threshold value.