Method and system for estimating a velocity vector field of uterine strain

Abstract

Some embodiments are directed to estimating a velocity vector field of a strain wave in a uterus For example, distances and delays in a strain wave between pairs of selected defined locations may be estimated over a time period. Local velocity vector may be computed from the estimated delays and distances. The velocity vector field may be visualized and / or characterized.

Description

[0001] METHOD AND SYSTEM FOR ESTIMATING A VELOCITY VECTOR FIELD OF

[0002] UTERINE STRAIN

[0003] TECHNICAL FIELD

[0004] The presently disclosed subject matter relates to a method for estimating a velocity vector field of a strain wave in a uterus, a system comprising for estimating a velocity vector field of a strain wave in a uterus, and a computer storage medium.

[0005] BACKGROUND

[0006] Uterine mechanical behavior in the form of uterine contraction plays a key role in reproduction. Irregular patterns of uterine contractility associate with dysfunctions of the uterus.

[0007] The diagnosis of uterine dysfunctions, such as adenomyosis, today mainly based on the assessment of some anatomical measures by MRI, based on the MUSA criteria.

[0008] Uterine activity, characterized primarily by periodic contractions within the sub endometrial layer of the myometrium, is believed to correlate with various uterine dysfunctions during pregnancy and natural menstrual cycles. These dysfunctions, including dysmenorrhea, endometriosis, and adenomyosis [3, 6, 8], may ultimately contribute to infertility. An accurate and objective assessment of uterine contractility outside pregnancy may provide valuable insights for improving clinical diagnosis and treatment of these conditions.

[0009] Several studies have attempted to quantify and characterize uterine contractions to reveal the complex mechanism of the uterus. Transabdominal electrohysterogram (EHG) has been proposed as a noninvasive tool to investigate the presence and propagation patterns of uterine electrical activity that triggers and drives the mechanical contraction of the uterus during [10, 11] and outside pregnancy [7, 14], However, the EHG signals fail to precisely localize the origin of uterine contractions within the uterus

[0012] , which limits the accurate characterization of the uterine electrophysiological activity.

[0010] More recently, two-dimensional transvaginal ultrasound (2D TVUS) has been employed as a non-invasive potential tool for quantitative assessment of uterine motion. Dedicated speckle tracking algorithms were developed to track the propagation of uterine contraction throughout the acquired TVUS recording [4, 15], Following this, strain analysis was applied to quantify the deformation of the subendometrial layer resulting from these contractions. Subsequently, a set of features, including contraction amplitude, frequency, energy, velocity, and coordination were extracted to either distinguish four menstrual phases [4, 15] or predict the success of IVF treatment prior the embryo transfer [4, 13],

[0011] Unfortunately, several disadvantages are still associated with the detection and display of strain fields in uteri. Strain fields provide a less intuitive and immediate visual representation of uterine motion. The direction and magnitude of uterus movement is not directly observable, making it hard to understand contraction dynamics. Strain fields do not allow a detailed quantitative analysis of the speed and direction of uterus displacement. This limitation makes it challenging to identify patterns such as areas of high or low activity, which may correlate with specific clinical conditions.

[0012] Strain measurements do not provide information on the evolution of uterine strain over time. This hinders understanding the temporal dynamics of uterine contractions, which is important for identifying abnormalities. Using strain fields alone makes it harder for clinicians to assess the efficiency and effectiveness of uterine contractions.

[0013] References to known systems are collected near the end of the description.

[0014] SUMMARY

[0015] A method for estimating a velocity vector field of a strain wave in a uterus, a system comprising for estimating a velocity vector field of a strain wave in a uterus, and a computer storage medium are provided as defined in the claims. Specific embodiments of the invention are set forth in the dependent claims.

[0016] Deriving a velocity vector field for a uterus is advantageous, and in particular, especially so for the study of uterine contractions. Embodiments may be used for quantitative ultrasound analysis of the uterine contractions (peristalsis), especially for diagnostic purposes and for the prediction of in-vitro fertilization success.

[0017] Velocity vector fields offer a more intuitive and immediate visual representation of uterine motion compared to strain fields. The direction and magnitude of uterus movement can be directly observed. Accordingly, contraction dynamics can be understood from a visualization in a manner not possible with a strain map.

[0018] Velocity vector fields allow for a qualitative as well as quantitative analysis of the speed and direction of uterus strain, which in turn indicates how the strain evolves over time. This is particularly useful for identifying patterns, such as areas of high or low activity, which may correlate with clinical conditions.

[0019] By converting strain fields to velocity vector fields, clinicians can better assess the efficiency and effectiveness of uterine contractions. The velocity vector field gives an impression of coordinated movement vs disorder in the uterus. This provides help in diagnosing issues such as adenomyosis, cesarean niches, fibromas, myomas, and endometriosis, and predicting the need for interventions, and monitoring the response to treatments. Velocity vector fields furthermore may be used to predict the success of in- vitro-fertilization success.

[0020] Estimating a velocity vector field may use 2D or 3D ultrasound measurements of the uterus. For example, a strain map may first be derived from the ultrasound measurements from which a velocity vector field is derived in turn.

[0021] The inventors found that using 3D ultrasound measurements was advantageous. Quantifying uterine contraction with 2D TVUS has some important limitations. Out-of-plane (OOP) motion is a frequent phenomenon during in-vivo 2D TVUS acquisition. It is primarily caused by the movement of the patient and the probe, but it may also be influenced by uterine contractions resulting into uterus displacement that is perpendicular to the acquisition plane [4], Once OOP motion occurs, the speckle pattern being tracked will undergo significant decorrelation or may even entirely disappear, ultimately resulting in tracking failure. Thus, the recordings with frequent OOP motion have to be excluded from the quantitative analysis, leading to a high rejection rate of the acquired recordings. Furthermore, uterine contractions propagate along the endometrium in three dimensions. Analyzing uterine contractions from a single 2D plane may limit our comprehensive understanding of uterine motion.

[0022] 3D TVUS shows the potential to improve diagnostic accuracy for conditions like adenomyosis and endometrial cancer.

[0023] A system for deriving a velocity vector field is an electronic device, e g., a suitably adapted computer, ultrasound device, etc. A method for deriving a velocity vector field may be executed on an electronic device, e.g., a computer, an ultrasound device, etc. A person skilled in the art will appreciate that the method may be applied to multidimensional image data, e.g., to two-dimensional (2D), three-dimensional (3D) or four-dimensional (4D) images, acquired by various acquisition modalities such as, but not limited to, standard X-ray Imaging, Computed Tomography (CT), Magnetic Resonance Imaging (MRI), Ultrasound (US), Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT), and Nuclear Medicine (NM).

[0024] An embodiment of the method may be implemented on a computer as a computer implemented method, or in dedicated hardware, or in a combination of both. Executable code for an embodiment of the method may be stored on a computer program product. Examples of computer program products include memory devices, optical storage devices, integrated circuits, servers, online software, etc. Preferably, the computer program product comprises non-transitory program code stored on a computer readable medium for performing an embodiment of the method when said program product is executed on a computer.

[0025] In an embodiment, the computer program comprises computer program code adapted to perform all or part of the steps of an embodiment of the method when the computer program is run on a computer. Preferably, the computer program is embodied on a computer readable medium.

[0026] BRIEF DESCRIPTION OF DRAWINGS

[0027] Further details, aspects, and embodiments will be described, by way of example only, with reference to the drawings. Elements in the figures are illustrated for simplicity and clarity and have not necessarily been drawn to scale. In the figures, elements which correspond to elements already described may have the same reference numerals. In the drawings,

[0028] Figure 1 schematically shows an example of an embodiment of a device configured for estimating a velocity vector field of a strain wave in a uterus,

[0029] Figure 2 schematically shows an example of an embodiment of a plurality of defined locations distributed over a uterus,

[0030] Figures 3a.1 and 3a.2 schematically show examples of an embodiment of a visualization of a velocity vector field,

[0031] Figures 3b.1 and 3b 2 schematically show examples of an embodiment of a visualization of a velocity vector field, Figure 4a-4c schematically show examples of an embodiment of transversal slices selected along the longitudinal direction,

[0032] Figure 4d schematically shows an example of an embodiment of evenly spaced lines starting from a centroid of the endometrium border,

[0033] Figure 5a schematically shows an example of an embodiment of transverse, sagittal, and coronal planes of endometrial cavity,

[0034] Figure 5b schematically shows an example of an embodiment of a 3D endometrium models and grid positioning,

[0035] Figure 5c schematically shows an example of an embodiment of a 3D endometrium models and grid positioning.

[0036] Figures 6a and 6b show examples of an embodiment of a 2D contraction propagation map.

[0037] Figures 7a and 7b show an example of an embodiment of detected propagation paths,

[0038] Figures 8a and 8b show an example of an embodiment of detected propagation paths,

[0039] Figure 9 schematically shows an example of an embodiment of a method for estimating a velocity vector field of a strain wave in a uterus,

[0040] Figure 10a schematically shows a computer readable medium having a writable part comprising a computer program according to an embodiment,

[0041] Figure 10b schematically shows a representation of a processor system according to an embodiment.

[0042] Reference signs list

[0043] The following list of references and abbreviations corresponds to figures 1, 2, and 10, and is provided for facilitating the interpretation of the drawings and shall not be construed as limiting the claims.

[0044] 110 a velocity vector field estimation device

[0045] 111 a processor system

[0046] 112 storage

[0047] 113 communication interface

[0048] 200 a uterus 210 a kernel

[0049] 1000, 1001 a computer readable medium

[0050] 1010 a writable part

[0051] 1020 a computer program

[0052] 1110 integrated circuit(s)

[0053] 1120 a processing unit

[0054] 1122 a memory

[0055] 1124 a dedicated integrated circuit

[0056] 1126 a communication element

[0057] 1130 an interconnect

[0058] 1140 a processor system

[0059] DESCRIPTION OF EMBODIMENTS

[0060] While the presently disclosed subject matter is susceptible of embodiment in many different forms, there are shown in the drawings and will herein be described in detail one or more specific embodiments, with the understanding that the present disclosure is to be considered as exemplary of the principles of the presently disclosed subject matter and not intended to limit it to the specific embodiments shown and described.

[0061] In the following, for the sake of understanding, elements of embodiments are described in operation. However, it will be apparent that the respective elements are arranged to perform the functions being described as performed by them.

[0062] Further, the subject matter that is presently disclosed is not limited to the embodiments only, but also includes every other combination of features described herein or recited in mutually different dependent claims.

[0063] Figure 1 schematically shows an example of an embodiment of a velocity vector field estimation device 110. Velocity vector field estimation device 110 is configured to obtain a time-dependent strain field indicating an estimate of strain at a plurality of defined locations distributed over a uterus during a time period, and to compute a velocity vector field therefore. The velocity vector field may indicate the velocity with which a stain wave progresses through the uterus.

[0064] Velocity vector field estimation device 110 may comprise a processor system 111, a storage 112, and a communication interface 113.

[0065] Velocity vector field estimation devices 110 is represented here as a single device, though it could just as well be implemented as systems, e.g., a geographically distributed system, e.g., a cloud computing system, e.g., a system comprising multiple computers.

[0066] Storage 112 may be, e.g., electronic storage, magnetic storage, etc. The storage may comprise local storage, e.g., a local hard drive or electronic memory. Storage 112 may comprise non-local storage, e.g., cloud storage. In the latter case, storage 112 may comprise a storage interface to the non-local storage. Storage may comprise multiple discrete sub-storages together making up storage 112

[0067] Storage 112 may be non-transitory storage. For example, storage 112 may store data in the presence of power such as a volatile memory device, e.g., a Random Access Memory (RAM). For example, storage 112 may store data in the presence of power as well as outside the presence of power such as a non-volatile memory device, e.g., Flash memory. Storage 112 may comprise a volatile writable part, say a RAM, a non-volatile writable part, e.g., Flash. Storage may comprise a non-volatile non-writable part, e.g., ROM, e.g., storing part of the software.

[0068] In the various embodiments of communication interface 113, the communication interface may be selected from various alternatives. For example, the interface may be a network interface to a local or wide area network, e.g., the Internet, a storage interface to an internal or external data storage, an application interface (API), etc.

[0069] The device 110 may communicate internally, with other systems, with other devices, external storage, input devices, output devices, and / or one or more sensors over a computer network. The computer network may be an internet, an intranet, a LAN, a WLAN, a WAN, etc. The computer network may be the Internet. The device 110 comprise a connection interface which is arranged to communicate within device 110 or outside of device 110 as needed. For example, the connection interface may comprise a connector, e.g., a wired connector, e.g., an Ethernet connector, an optical connector, etc., or a wireless connector, e.g., an antenna, e.g., a Wi-Fi, 4G or 5G antenna. The communication interface 113 may be used to receive digital data, e.g., a time-dependent strain field, or ultrasound measurements from which the time-dependent strain field may be derived. The communication interface 113 may be used to send digital data, e.g., the velocity vector field, a visualization thereof, a characterization thereof, etc.

[0070] Velocity vector field estimation device 110 may have a user interface, which may include well-known elements such as one or more buttons, a keyboard, display, touch screen, etc. The user interface may be arranged for accommodating user interaction for performing actions such as, initiating the computation of a velocity vector field, to show and / or manipulate a visualization, to initiating the computation of a characterization, and / or to show and / or manipulate the characterization.

[0071] The execution of device 110 may be implemented in a processor system 111. The device 110 may comprise functional units to implement aspects of embodiments. The functional units may be part of processor system 111. For example, functional units shown herein may be wholly or partially implemented in computer instructions that are stored in a storage of the device and executable by processor system 111.

[0072] Processor system 111 may comprise one or more processor circuits, e.g., microprocessors, CPUs, GPUs, etc. Device 110 may comprise multiple processors. A processor circuit may be implemented in a distributed fashion, e g., as multiple subprocessor circuits. For example, device 110 may use cloud computing.

[0073] Typically, the velocity vector field estimation device 110 comprises one or more microprocessors which executes appropriate software stored at the device; for example, that software may have been downloaded and / or stored in a corresponding memory, e.g., a volatile memory such as RAM or a non-volatile memory such as Flash.

[0074] Instead of using software to implement a function, device 110 may, in whole or in part, be implemented in programmable logic, e g., as field-programmable gate array (FPGA). The device may be implemented, in whole or in part, as a so-called applicationspecific integrated circuit (ASIC), e.g., an integrated circuit (IC) customized for their particular use. For example, the circuits may be implemented in CMOS, e.g., using a hardware description language such as Verilog, VHDL, etc. In particular, velocity vector field estimation device 110 may comprise circuits, e g., for graphical processing, and / or arithmetic processing.

[0075] In hybrid embodiments, functional units are implemented partially in hardware, e.g., as coprocessors, and partially in software stored and executed on the device. Figure 9 schematically shows an example of an embodiment of a method 500 for estimating a velocity vector field of a strain wave in a uterus. Method 500 comprises obtaining (510) a time-dependent strain field indicating an estimate of strain at a plurality of defined locations distributed over a uterus during a time period. Method 500 may be executed, e.g., on a computer, on an ultrasound device, a workstation, etc. Method 500 may be executed on device 110 described herein.

[0076] There are various ways to obtain the time-dependent strain field. For example, the time-dependent strain field may be received from a further device configured for computing and / or measuring the time-dependent strain field. The time-dependent strain field may also be computed locally, e.g., on a device such as device 110.

[0077] One way to derive a time-dependent strain field is described in the International Patent Application publication WO2019053249A1, which explains a method for the derivation of uterine strain in 2D and 3D such that the estimated strain is relative to the uterine anatomy.

[0078] The plurality of defined locations may be distributed over one or more layers of the uterus, e.g., the subendometrial layer of the uterus.

[0079] Strain can be estimated by applying speckle tracking techniques over the acquired ultrasound B-mode video loops. Multiple locations can be tracked that are distributed over specific uterine layers and / or structures, such as the endometrium, the sub -endometrial layer or the myometrium.

[0080] The derivation of the velocity vectors can be based on any different strain derived according to the uterine anatomy and orientation, with the longitudinal direction defined along the endometrial line and the transversal or radial directions defined perpendicular to it. The different derived strains along these directions are named as longitudinal, transversal (or radial), and circular strain, with the latter being perpendicular to both the longitudinal and radial directions. Directionless strain such as area and volume strains can be derived by measuring the variations in the area (2D) or volume (3D) determined by the tracked markers.

[0081] It proved not possible to derive the velocity vector field directly from the observed relative movement, e.g., movement of the speckle patterns. It was an insight of the inventor to obtain the velocity vector field, from the time-dependent strain field. An advantageous method is to first derive a time-dependent strain field, and based on the estimated local strains over time and over multiple locations, derive the velocity vector field of the uterus, in particular at some point during a uterine contraction. The velocity vector field is itself time dependent. For example, in an embodiment, the velocity vector field is computed for a particular point in time, e.g., as indicated by a user, e.g., at a predefined point during a uterine contraction, e.g., at 25% or 50%, etc., of the contraction. Computing the velocity vector field takes into account the evolution of the strain field in a time segment or time window around the particular point in time. For example, given a particular point in time, a time segment may be defined around it. The size of the time segment may be determined empirically, and depends, inter alia on the temporal resolution of the measurements, and the desired resolution of the velocity vector field.

[0082] Alternatively options are to compute the velocity vector field for multiple points in time, e.g., multiple points in time during a contraction, and to show an animation of the changing velocity vector field, or to compute therefrom an average velocity vector field, which may also be showed.

[0083] One strain wave lasts for about 20 seconds. An embodiment may show how the velocity vector field changes over that time, or part thereof. For example, a standard acquisition may last longer, e.g., 4 minutes, and my capture several contractions. Start and end of contractions can be identified by a medical segmentation model, e.g., a machine learning model, or may be indicated by a human user of the system.

[0084] In practice, the peristaltic waves are often continuous and have a sinusoidal behavior, and as such do not have a strictly defined start and end point. Accordingly, the system and / or its user may define a start and end of contractions so as to capture part of a contraction, a full contraction, or more than one full contraction, as desired. If needed, a start and end point may be defined with respect to the wave like behavior of the contractions, e.g., a start and end point may be defined as a minimum in strain activity, or the like.

[0085] In an embodiment, the strain field is relative to the uterus anatomy and orientation. For example, the strain data is aligned with the specific anatomical structure and natural orientation of the uterus. In particular, strain directions may be relative to the outer contour of the uterus or the endometrium and may be presented in radial, longitudinal, and circular directions. As mentioned the strain field may be computed locally; the segmentation of the uterus and / or endometrium may be received or may be computed locally.

[0086] For example, the strain field may be determined in a longitudinal direction defined along a longitudinal axis of the uterus, and / or a transversal direction and / or a circular direction. The longitudinal direction is along the length of the uterus, e.g., the endometrial line. The transversal and circular directions are perpendicular to the longitudinal direction. Strain in a particular direction can be derived from movement of markers, e.g., speckle patterns, along that particular direction.

[0087] Strain may also be directionless, e.g., area strain and / or volume strain. Directionless strains such as area and volume strains measure changes in the overall size Which strain type to use may be determined in the method. For example, a strain field may be computed for multiple strain types, e.g., for some or all of the mentioned strain types, and a strain type may be selected which shows a high, e.g., highest, signal to noise ratio. Which strain type to use may be selected by a user, or may be fixed in the software.

[0088] In an embodiment, the time-dependent strain field may be computed from 2D or 3D ultrasound measurements of the uterus. The ultrasound measurements have multiple frames. For example, an embodiment may comprise:

[0089] - determining an orientation of the uterus,

[0090] - generating a plurality of tracking points inside a region of interest of the uterus,

[0091] - tracking and estimating the displacement of the plurality of tracking points between frames of the ultrasound measurements;

[0092] - calculating, from a varying distance between subsets of the tracking points, at least one of: a transversal strain, a longitudinal strain, a circular strain, an area strain and a volume strain.

[0093] Figure 2 schematically shows an example of an embodiment of a plurality of defined locations distributed over a uterus 200. The defined locations may be the same as the tracking markers, but this is not necessary. The defined locations and / or tracking markers may be regularly distributed over the uterus, e.g., one or more layers of the uterus. The defined locations and / or tracking markers may be randomly distributed over the uterus, e.g., one or more layers of the uterus.

[0094] For example, the defined locations may be the locations over the transversal planes, with regular radial spacing around the endometrial line, where the local radial strains are estimated. Figures 4a-4d provide an example of such spacing.

[0095] Figure 2 provides an example of tracking markers distributed over the subendometrial layer of the uterus. The markers are determined by defining a number of parallel planes perpendicular to the endometrial line. In a possible application, for each plane, the interception with the uterine layer of interest (e.g., subendometrial layer) of the radii originating from the endometrial line, with a step of 10 degrees, determines the markers that are tracked for calculating the strain. When radial strain is to be calculated, two subsequent markers are tracked along each radius and their relative distance variation used to calculate the radial strain. The endometrial line may be defined as a longitudinal axis of the uterus. Alternatively or in addition, the tracking markers may be distributed over other layers of the uterus, e.g., more into the muscle, in the myometrium.

[0096] At reference 210, an example of 3x3 kernel is shown, which may be used for the estimation of the local velocity vector.

[0097] In the embodiment shown in figure 2, radial strains are measured. Note that these tracking markers are not sufficient to calculate radial strain as you need two markers along each radius to calculate the strain, e.g., from their relative distance variation. Therefore, to calculate the strain at these locations two radial markers per location were used. The location may correspond to one of the two markers, or may be associated to a location between the two markers used to calculate the strain.

[0098] More generally, although an association between defined locations and tracking points is typically used, this is not necessary. For example, strain may be interpolated at any location from strain computed for the tracking locations.

[0099] Given the time-dependent strain field, whether computed locally or received, a velocity vector field may be derived therefor. Method 500 comprises performing (520) a spatiotemporal analysis of the obtained strain field to determine a velocity vector field of a strain wave in the uterus. In particular, the velocity vector field of uterine contractions may be computed. In an embodiment, estimations are made of local propagation delays from which velocity vectors are computed. Method 500 comprises determining (530) local velocity vectors for the plurality of defined locations, the local velocity vectors together forming the velocity vector field.

[0100] For example, the local velocity vectors may be computed at all or part of the multiple defined locations and / or corresponding markers. For example, for any particular location, a local velocity vector may be computed. This allows for the estimation of the velocity vector field over the full uterine surface. Several approaches can be taken to derive the local velocity vectors. Determining a local velocity vector for a particular defined location may proceed as follows. First multiple defined locations are selected from the plurality of defined locations in a neighborhood of the particular defined location;

[0101] This can be done in various ways. For example, in an embodiment, the defined locations may be distributed in a grid, and the neighborhood of the particular defined location comprises a kernel of defined locations centered on the particular defined location. This is option is illustrated in figure 2.

[0102] Alternatively, the neighborhood of the particular defined location may comprise defined locations within a threshold distance to the particular defined location. In some situations, the two approaches may coincide. For example, a in gird taking the threshold distance at V2 times the unit grid length, will produce the same kernel as shown in figure 2. However, other types of grid, and in particular for randomized locations, the distance approach may be more flexible.

[0103] Using randomized locations is an advantage because it can mitigate potential artifacts that arise from regular grid patterns, leading to a more accurate representation of the strain and velocity fields. Randomized points can provide a more robust sampling of the uterine surface, reducing the risk of artefacts due to grid alignment issues. Randomization also allows the selection of only well-defined speckle patterns, enhancing the accuracy of the strain and velocity estimations. Additionally, randomization enables an increase in the number of points in areas of interest, thereby increasing the resolution and providing a more detailed analysis of those specific regions.

[0104] Once neighboring points are selected, the method may comprise estimating (532) delays in the strain wave between pairs of selected defined locations over all or part of the time period. It is possible to include each pair of selected defined locations, but this is typically not necessary, and a subset will suffice.

[0105] To estimate a delay between two given points, e.g., point 1 and point 2 neighboring locations, e.g., within the kernel. The following approach may be used. We assume that the same strain wave is recorded at the two points, such that a couple of strains yx(t) and y2(t) may be modelled as yx(t) = s(t) + nx(t) and y2(t) = s(t - r) + n2(t), with s(t) being the radial strain wave propagating through the kernel, T being the delay of the strain wave between the two locations, and nx(t) and n2(t) being additive noise. Based on the time-strain evolution at each location, the delay f between the strain measured for each pair of locations in the kernel (yx(t) and y2(t)) can be calculated by maximizing the expectation of the two signals as flj2= argmax£[y1(t)y2(t + t)]

[0106] Wherein t1;2denotes the estimated delay of the strain wave between points 1 and 2.

[0107] For example, estimating a delay in the strain wave between a pair of selected defined locations over all or part of the time period may comprise maximizing a utility function for a time delay. The utility function comprises a correlation between a strain over time at a first defined location of the pair according to the strain field, and a strain over time shifted over the time delay at a second defined location of the pair. For example, the utility function may be E[y1(t)y2(t + 1)], computed over a time period. Here E( ) is the expectation operator.

[0108] The time period may be the full contraction period, but may also be a part thereof. For example, delays can be computed for multiple consecutive, possibly partially overlapping time periods, to support determining the velocity vector field at multiple locations.

[0109] This operation can be performed by maximizing the cross-correlation between the two signals in the time domain or more efficiently by using the Fourier transform as

[0110] T1 2= arg where Vj (m) and K2(m) are the Fourier transforms of yx(t) and y2(t), respectively, and (•)* denotes the complex conjugate.

[0111] Interestingly, this approach also allows an estimation of the delay that is not limited by the sampling rate in time, e.g., frame rate, as you would have by crosscorrelation in time.

[0112] The method may further comprise estimating (533) distances between the pairs of selected defined locations. The distances may be expressed as vectors, e.g., indicating diction as well as magnitude. Distance may be expressed as Euclidian distance, but this is not needed. For example, the distance may be inter-pixel distance vectors. The distance vectors may be obtained from various sources, e.g., the ultrasound measurements, a modelling of the uterus, from the strain field, and so on.

[0113] The method may further comprise computing (534) the local velocity vector for the particular defined location from the estimated delays and distances. For example, given a distance vector and a delay, the product between an estimate for the local velocity vector should be close to the distance vector. For example, in an embodiment, the local velocity vector may be obtained by minimizing a loss function for the local velocity vector. The loss function comprising differences between the estimated distances between the pairs of selected defined locations and products of the local velocity vector and the estimated delays.

[0114] For example, by estimating the delay between each pair of locations surrounding the reference location at the center of the kernel, we obtain a vector f of collecting all the estimated delays in the kernel. The vector f may be represented as a row vector, e.g., as an array. For example, the vector may have length N(N-l) / 2, assuming all pairs are computed once for a set of N points. For example, in an embodiment, all points on the edge of the kernel are used, which in case of a 3x3 kernel, means N=8.

[0115] Given a 3D vector v defining the average local propagation velocity of the strain waves across the neighborhood of a particular location, e.g., with the kernel, and accounting for the distances between the locations in the three directions, collected in a matrix D of size 3xN(N-l) / 2, we can write

[0116] VT = D,

[0117] The velocity vector can then be calculated by minimizing an error function, defined in terms of L2 norm, between the actual distances and those estimated by multiplying the previously estimated delays with the velocity vector to be estimated as where Dtis the i-th column of D and f£the corresponding time delay.

[0118] The optimization can be further improved by introducing a weighted optimization in the error function, with weights wtthat are larger for larger si nal-to-noise ratio (SNR) of the measured signals y£;1(t) and y£;2(t) (estimation confidence). For example, one may take E[yi,i(t)yi,2(t + T)] WI = -

[0119] <JVYi.,i <TVYI.,2 The weight wtmay be taken as the value of the normalized cross correlation function at lag i, reflecting our confidence in the time-delay estimate. Herey.iandy. 2depict the standard deviations of the strain field at locations 1 and 2 of pair i, respectively.

[0120] The same approach may be used for locations that do not form a grid. For example, locations in a neighborhood of a particular define location may be selected, e.g., those within a threshold distance. For part of the pairs of points in the neighborhood, distances and delays may be computed. For example, all distances and delays for points different from the particular point at the center of the neighborhood may be taken. However, a smaller sampling may be used. For example, a set n of pairs may be selected, possibly randomly, e.g., at least 10 pairs, at least 20 pairs, etc. Distances between the selected pairs may be recorded in a matrix D and delays in an array f, the local velocity vector may then be obtained by minimizing the loss function as above.

[0121] Computing the velocity vector field has several applications. For example, the velocity vector field or a derivative thereof may be visualized. The visualization may be displayed on a display. The visualization may be displayed to a user.

[0122] Such velocity vector field can be used for visualization as well as for further characterization of uterine contractions. Average velocity or dominant propagation directions can be estimated and shown to the user.

[0123] For example, the velocity vector field characterized using various values or maps. The characterization may be displayed on a display. The characterization may be displayed to a user. The characterization may be reported.

[0124] In an embodiment, the velocity vector field may be determined for multiple subsequent parts of the time period. Possibly the subsequent parts are overlapping. For example, 5 second time periods, overlapping in periods of 0.5 seconds may be used to derive 40 velocity vector field for a typical 20 second contraction. The multiple velocity vector fields may be visualized in an animation. Different parameters may be used, e.g., time periods that are shorter than 5 seconds, e.g., 1 second, or longer, e.g. 10 seconds. The segments may overlap more, e.g., 1 second, or less, e.g., 0.1 second.

[0125] Such measures, which can be measured both globally or locally, and can also be estimated over a running time window to assess their dynamics, can serve as a valuable tool for the diagnosis of uterine dysfunctions, such as adenomyosis, cesarean niches, fibromas, myomas, and endometriosis, as well as for the prediction of in-vitro- fertilization success Various visualization are possible. For example, the velocity vector field visualization may comprise an overlay, which may be shown overlayed over a visualization of the uterus, e.g., a 3D representation. This is not needed, the velocity vector field may be shown as a map, e.g., a rectangular map.

[0126] Figures 3a.l and 3a.2 schematically show examples of an embodiment of a visualization of a velocity vector field. The velocity vector field of the uterine contractions are derived from spatiotemporal analysis of local strains. The estimation of a velocity vector represents the propagation of uterine strain waves at each location. The velocity vector field visualization comprises a plurality of arrows indicating the velocity vector field, e.g., at the defined locations or a subset thereof, and / or the magnitude of the velocity.

[0127] Figures 3b.1 and 3b.2 schematically show examples of an embodiment of a visualization of a velocity vector field. We also propose a representation of the propagation pattern of the uterine contractions in terms of main streamlines connecting the local velocity vectors. Streamlines give an indication of dominant direction(s) to the user. Example of streamlines derived from the velocity vector fields are shown in figures 3b.1 and 3b.2.

[0128] The streamlines connecting the estimated vectors provide a clear rendering of the propagation path of the uterine contractions, e.g., peristaltic waves. Rendering the propagation path may use of colors, e.g., to indicate a magnitude of velocity, strain, or the like.

[0129] The visualization comprises a rendering of a propagation path of contractions in the uterus. For example, the propagation path may be computed from the velocity vector field. An algorithm is presented herein to derive propagation path, even without first explicitly deriving the velocity vector field. For example, a propagation path may be started from one or more defined points, e.g., defined by the software, the user, and the like.

[0130] The derivation of the vector field and related streamlines can be performed on the full ultrasound video loop as well as over a time window of, e.g., 20 s. Other time periods are also possible. Such time window can then be shifted over time, with a given overlap, so as to capture and render how the vector field varies over time. This allows for investigating the dynamic variations of the propagation pattern of uterine contractions and peristaltic waves, providing additional relevant information for assessing the uterine condition. Several characterizations are possible. For example, it was found that useful characterizations include derivation of entropy, mutual information, vorticity, and dyssynchrony measures.

[0131] For example, in an embodiment, characterization of the velocity vector field comprises computing a global uterine activity parameter from the velocity vector field quantifying the velocity vector field. For example, parameter may be one or more of: average velocity, average velocity magnitude, and / or average velocity direction, a vorticity, vorticity map, a dyssynchrony index, where the latter can be derived by calculating the spread of the timing (phase) of the peristaltic waves from the estimated mutual delays over a certain subset of locations (region) distributed e.g. in the circular direction (around the uterine main longitudinal).

[0132] For example, in an embodiment, the characterization may comprise

[0133] - computing an entropy of the velocity vector field globally and / or within a kernel and / or comprising determining a probability distribution of velocity vectors occurring globally and / or within a kernel, computing an entropy of the probability distribution, and / or

[0134] - computing a conditional entropy and / or mutual information between one or more velocity vectors and between one or more neighboring velocity vectors, and / or

[0135] - computing an entropy map and / or mutual information map, and visualizing the map, e.g., as an overlay.

[0136] The local estimations of the propagation delays and velocity vectors enable the derivation of relevant parameters describing the uterine contraction activity. In particular, we can analyze the entropy of the vector field within the kernel (or globally) as a measure of its irregularity and complexity. Shannon’s entropy is defined as with P(v = [vx, vy, vz]) being the three-dimensional probability distribution of a certain velocity vector occurring. This can be calculated by counting the number of vectors within the kernel that fall within a certain bin i, j, and k in the three directions. As a result, the estimated Shannon entropy H(V) is given as

[0137] H(F) = - i,j,k P(i,j, fc) log P (i,j, k

[0138] The use of the information theory for assessing the irregularity of the velocity vector field can be further extended by deriving, e.g., the conditional entropy or mutual information and, therefore, looking at the statistical dependency of the velocity vectors between neighboring locations. Such neighboring locations can be defined as the center of the kernel compared to its periphery, but different subdivisions, such as in sectors, are also possible.

[0139] The standard definitions of conditional entropy H(V1,V2) and mutual information I(V1,V2) are given as with P(.,.) defining the joint probability. All probability distributions can again be estimated from the vector field, similar to Shannon entropy, by defining a number of bins, i,j, k, in the three directions and counting the number of vectors falling within such bins. Exploiting the dynamic estimation of the vector field over a running temporal window, the estimation of entropy, conditional entropy and mutual information can also benefit from accounting for multiple estimations over time, increasing the statistics of the estimates.

[0140] Next to parameters derived from information theory, the irregularity and complexity of the vector field can also be assessed by other metrics such as vorticity.

[0141] Vorticity, e.g., the curl of the velocity field, quantifies the circular, swirling motion characterizing a vortex. The formation of abnormal vortices is associated with dysfunction of the uterus. Mathematically, flow vorticity can be defined as the curl of the velocity vectors. The overall degree of vorticity may be estimated as the mean vorticity.

[0142] When we measure vorticity over the 3D surface where the velocity vector is estimated, this can be estimated as . . dvi dvcvorticity

[0143] J= - - — , de dl with vtand vcrepresenting the longitudinal and circular components of the contraction propagation, respectively.

[0144] Vorticity can be assessed globally over the entire surface or locally over a predefined kernel, allowing for generating a map of vorticity over the full 3D surface of the uterus or part thereof.

[0145] Additionally, it is possible to assess the degree of dyssynchrony of the contraction by, e g., evaluating the symmetry of the contraction at each transversal plane for varying angle around the endometrial line. The maximum delay between waves can be used to provide an index of dyssynchrony in the overall propagation from fundus to cervix and vice versa. Such analysis can also be extended to multiple transversal planes up to covering the full surface.

[0146] The visualization and / or characterization may be used to diagnose uterine dysfunctions, e.g., one or more of: adenomyosis, cesarean niches, fibromas, myomas, and endometriosis. The visualization and / or characterization may be used for the prediction of in-vitro-fertilization success.

[0147] By providing an analysis of the uterine contraction patterns, users can identify abnormal behaviors indicative of specific conditions.

[0148] For example, fibromas and myomas, as well as caesarean scar niches, may cause localized disruptions in the contraction patterns. This can also be the case for adenomyosis, although adenomyosis can often be more diffused over tissue rather than localized. These disruptions can be characterized by analyzing the velocity vector fields for anomalies in specific areas. The entropy map can also reveal areas with abnormal propagation of the waves, which may be caused by tissue consistency, providing further evidence for the presence of abnormalities, etc. Endometriosis can also lead to irregular contractions and abnormal tissue structures, which can be detected as deviations in these visualizations and / or characterizations.

[0149] For predicting in-vitro fertilization (IVF) success, the characterization of uterine contractions can provide valuable insights into the readiness of the uterus (uterine receptivity) for embryo implantation. By analyzing the synchronization and regularity of contractions, as well as the overall uterine activity, clinicians can assess the likelihood of successful embryo implantation.

[0150] Additionally, the visualization and / or characterization of the velocity vector fields can aid in monitoring the progress of treatments for these conditions. By comparing pre- and post-treatment data, clinicians can evaluate the effectiveness of interventions.

[0151] Below several further optional refinements, details, and embodiments are illustrated.

[0152] Data Acquisition

[0153] 2D and 3D TVUS data were collected at the Catharina Hospital Eindhoven from a healthy volunteer in the late follicular (LF) phase and a patient undergoing in-vitro fertilization (IVF) treatment using a Samsung-Medison ultrasound scanner WS80A equipped with a transvaginal V5-9 probe. The center frequency of the probe was 5.6 MHz. 2D TVUS recordings were acquired in the sagittal plane. The acquisition frame rate of the 2D TVUS recording was 30 frames per second and the volume rate of the 3D acquisition was 1.0 volume per second. The time length for the 2D acquisition was 4 minutes. The 3D acquisition time length was 41 s and 38 s for the healthy volunteer and IVF patient, respectively.

[0154] All 3D volumes were exported in MVL (Medison 3D Volume Data) format from the ultrasound scanner, and then converted to RAW format with a dedicated Cartesian Export tool provided by Samsung. After conversion, the size of the individual 3D volume was 600 X 600 X 600 pixels, with pixel size equal to 0.240 mm. After that, volumes in RAW format were converted to MAT (MATLAB fde) format by MATLAB for further analysis.

[0155] Propagation detection with 3D TVUS

[0156] In order to fully utilize the advantages of 3D data, the geometrical model of the endometrium was extracted from the 3D volume, and speckle tracking was then applied to investigate uterine contractions across the entire model. Eventually, propagation paths were detected based on strain analysis.

[0157] 3D model construction

[0158] To position tracking markers (TMs) along the endometrial cavity, a 3D model of the endometrial cavity was first created as follows,

[0159] • The endometrium cavity was divided into slices with a 10-pixel interval in different transverse planes (see Fig. 4a-4c, N slices with visible endometrium were selected.

[0160] Figure 4a-4c schematically shows an example of an embodiment of a multiple transversal slices selected along the longitudinal direction from the healthy volunteer.

[0161] On each selected slice, the endometrium border was manually drawn, (b) On each selected slice, 36 lines with 10 degrees interval starting from the centroid of the endometrium border were generated. TMs were placed at the cross-section points of the lines and border.

[0162] Figures 4a-4c show slices 290, 300 and 310, which progress in the transversal direction.

[0163] Figure 4d schematically shows an example of an embodiment of a 36 lines with 10 degrees interval starting from the centroid of the endometrium border On each selected slice, 36 lines with 10 degrees interval starting from the centroid of the endometrium border were generated. TMs were placed at the cross-section points of the lines and border.

[0164] • At each slice, the contour of the endometrium was manually segmented. Starting from the centroid of the contour, 36 lines were generated with a 10-degree interval (see Fig. 4d).

[0165] • The TMs were determined at the cross-section of these lines and the endometrium contour.

[0166] In the figures, the endometrium contour was manually drawn. Automatic segmentation of the endometrium may alternatively be used.

[0167] Figure 5a schematically shows an example of an embodiment of transverse, sagittal, and coronal planes of endometrial cavity.

[0168] Figure 5b schematically shows an example of an embodiment of a 3D endometrium models and grid positioning (black dots) for the healthy volunteer. Figure 5c schematically shows an example of an embodiment of a 3D endometrium models and grid positioning (black dots) for an IVF patient.

[0169] C and F represent cervix and fundus, respectively. A, P, RL, and LL represent anterior, posterior, right lateral, and left lateral sides of the endometrium, respectively.

[0170] 3D Ultrasound Speckle Tracking

[0171] 3D iterative OF based on Lukas-Kanade method was implemented to track the movement of all TMs [9], The fundamental hypothesis of OF is that the intensity, I, of a certain voxel located at (x, y, z, t) (x, y, z representing the space coordinates and t the time instant) in the reference volume does not change after a displacement (A x,A y,A z) in a period A t, i.e.,

[0172] I(x,y,z, t) = I(x + Ax,y + Ay, z + Az, t + At). (1)

[0173] Then, the velocity (v) of the selected voxel in x, y,z directions can be estimated by least squares, under the assumption that the flow was constant within a cubic window centered around the voxel: di . di di di _ y x ~ — F vv- — F v dx y dyz7- — F — — 0. dz dt (2) Strain Analysis

[0174] Radial strain (RS) was computed over time as the change in the distance between each TM on the endometrial contour and the corresponding centroid point normalized with respect to their initial distance. Other choices are possible, for example, a second marker along the radius could be used.

[0175] A fourth-order Butterworth bandpass filter was applied to the RS signals to minimize the interference from the respiration and probe motion [4], From the literature, the contraction frequency during a normal menstrual cycle varies from 0.5 - 4.1 contractions per minute [2], thus the cut-off frequencies of the filter were determined as fc= [0.0083,0.0683] Hz.

[0176] Velocity vector field

[0177] The obtained data may be used to derive one or more streamlines, possibly even a single streamline. Interestingly, this can be done with or without computing the velocity vector field.

[0178] If desired, the strain field, in this case, radial strain field may be used for estimating a velocity vector field, e.g., as indicated herein. The velocity vector field may be used to compute streamlines if desired. This approach is advantageous if many streamlines will be computed, e.g., for a visualization.

[0179] For example, in an embodiment, local velocity vectors are determined for the plurality of defined locations, by estimating both delays in the strain wave arriving between pairs of selected locations, and distances between the selected locations.

[0180] The velocity vector field can be used to determine a propagation path of the uterine contractions (peristaltic waves). For example, propagation paths may be obtained by generating streamlines, e.g., starting from a selected starting point. For example, one or more starting points may be selected by user, and / or may be automatically identified. For example, a plurality of starting points for a plurality of propagation path may be selected, e.g., distributed in an area of the uterus, e.g., in the cervix.

[0181] Propagation detection

[0182] As an alternative to computing the velocity vector field, a propagation path may be detected directly. The algorithm does not first compute the velocity vector field, though it may also use a cross-correlation.

[0183] Positive strain values indicate uterine relaxation, whereas negative values indicate contraction. Subsequently, the contraction frequency (CF) was estimated based on the number of zero-crossings per minute, representing the transition from relaxation to contraction or vice versa. The averaged CF was determined by calculating the mean of the CF estimated from all RS signals.

[0184] Firstly, we identify a point (xn, yn, zn) as the origin of contraction propagation based on preliminary inspection of RS signal variation. In an ideal scenario, the wavefront pattern remains consistent during propagation. Next, our objective is to determine a point (xn+1, yn+i>zn+i) likely to be the next location of propagation within a predefined search area that covers the current slice and slices below and above, and spans 10 degrees on both sides. This is determined by searching for the maximum normalized crosscorrelation between the RS signals at the current point (xn,yn,zn) and a potential next point (xn+1,yn+1,zn+1), given as where t0represents the initial time of the starting point, T is the time interval to evaluate cross-correlation, and tn— 0,1,2, ... , T denotes the delay time between two points on the path. The propagation path of a specific contraction ends when the maximal cross-correlation drops below 0.95.

[0185] For example, a method for determining a propagation path is provided, the method comprising

[0186] - selecting a point (xn,yn, zn) as the origin of contraction propagation

[0187] - repeatedly determining a next point (xn+1, yn+1, zn+1) likely to be a next location of propagation within a predefined search area.

[0188] The predefined search area may for example comprise the current slice and slices below and above, and spans 10 degrees on both sides.

[0189] Determining a next point may, for example, comprise searching for the maximum normalized cross-correlation between the RS signals at the current point (xn,yn, znand a potential next point (xn+1, yn+i,zn+i -

[0190] Results

[0191] Figures 6a and 6b show an example of an embodiment of a 2D contraction propagation map based on the RS signals extracted from the anterior wall of the endometrium from the healthy volunteer, and the IVF patient, respectively.

[0192] This plot shows a clear propagation from C2F for the healthy volunteer, while the propagation direction keeps changing from C2F to F2C for the IVF patient. Figures 7a and 7b show an example of an embodiment of detected propagation paths for the healthy volunteer. Figure 7a starts from Px(right lateral in cervical area). Figure 7b Starts from P2(anterior side in mid-cavity area). A clear propagation direction from C2F can be seen in the two paths.

[0193] Figures 8a and 8b show an example of an embodiment of detected propagation paths for the IVF patient. Figure 8a starts from P3(anterior side in mid-cavity area). Figure 8b starts from P4(posterior side in cervical area). Similarly, clear propagation direction from C2F can be seen in two paths.

[0194] As shown in Fig. 6a-b, the 2D contraction propagation maps are derived from the RS signals extracted from the anterior wall of the endometrium from the healthy volunteer (see Fig. 6(a)) and the IVF patient (see Fig. 6(b)). Positive RS, shown in red, indicates uterine relaxation, and negative RS, shown in blue, indicates uterine contraction. Clear cervix-to-fundus (C2F) propagation can be observed in the healthy volunteer. This finding is consistent with the anticipated uterine propagation direction in healthy women during the late follicular phase [5], However, a more complex propagation pattern, characterized by standing waves and alterations in the propagation direction, was observed in the IVF patient.

[0195] The endometrium cavity of the healthy volunteer and the IVF patient were divided into 36 and 18 slices, respectively. After the manual segmentation, the transverse, sagittal, and coronal planes of the endometrial cavity and models of the uterus are shown in Fig. 5a, and b-c, respectively. The average contraction duration of the healthy volunteer and the IVF patient were 12 s and 14 s, respectively. To visualize the propagation paths in 3D, we preliminarily inspected the start time of an obvious contraction (t0= 5 s for both subjects). Two starting points were selected within the contraction region in different anatomical orientations. As shown in Fig. 7a-b, the detected propagation lasted 15 s starting from P4(right lateral in cervical area, see Fig. 7(a)) and 10 s from P2(anterior side in mid-cavity area, see Fig. 7(b)). Similarly, in Fig. 8, the detected propagation lasted 8 s starting from P3(anterior side in mid-cavity area, see Fig. 8(a)) and 6 s from P4(posterior side in cervical area, see Fig. 8(b)).

[0196] A reliable evaluation of uterine motion is essential for delving into its underlying mechanisms. Traditional propagation detection methods using 2D TVUS fail to fully capture the complex contraction patterns of the uterus in three dimensions, thereby limiting the derived insight. Hence, we propose a 3D propagation detection method to enhance our understanding. As proof of concept, 2D and 3D TVUS recordings were obtained from a healthy volunteer and an IVF patient. 2D / 3D OF was employed to track the uterine motion and radial strain analysis was derived from motion signals to measure regional deformation around the endometrium.

[0197] As shown in Fig. 6(a) and Fig. 7a-b, the proposed 3D propagation detection method identified clear C2F propagation, especially along the anterior side, in the healthy volunteer during the late follicular phase, a finding confirmed by the 2D detection method. While 2D TVUS provided limited insights due to its reliance on a single plane, primarily from the mid-sagittal plane, 3D TVUS revealed a more detailed picture. Specifically, we observed that the contraction not only propagated from C2F but also exhibited a twisting motion from the right-lateral to the anterior side (Fig. 7a-b).

[0198] However, uterine contractions do not always propagate along the anterior or posterior wall of the endometrium, as can be captured by 2D TVUS. OOP motion may lead to false detection of the propagation pattern, thereby hindering the analysis of the propagation direction. For example, the IVF patient included in the study exhibited a complex uterine propagation pattern with the 2D detection method (see Fig. 6(b)), but in 3D, primary C2F propagation, even though with a twisting fashion, was again observed (see Fig. 8a-b).

[0199] References

[0200] [1] Marina Paula Andres, Giuliano Moyses Borrelli, Juliana Ribeiro, Edmund Chada Baracat, Mauricio Simoes Abrao, and Rosanne M Kho. Transvaginal ultrasound for the diagnosis of adenomyosis: systematic review and meta-analysis. Journal of minimally invasive gynecology, 2 (2): 257-264, 2018.

[0201] [2] Carlo Bulletti, Dominique de Ziegler, Valeria Polli, Lidia Diotallevi, Elena Del Ferro, and Carlo Flamigni. Uterine contractility during the menstrual cycle. Human reproduction, 15(suppl_l):81-89, 2000.

[0202] [3] Carlo Bulletti and Dominique De Ziegler. Uterine contractility and embryo implantation. Current Opinion in Obstetrics and Gynecology, 17(3):265-276, 2005.

[0203] [4] Yizhou Huang, Connie Rees, Federica Sammali, Celine Blank, Dick Schoot, and Massimo Mischi. Characterization of uterine peristaltic waves by ultrasound strain analysis. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 69(6):2050-2060, 2022. [5] Marga M Ijland, Johannes LH Evers, Gerard AJ Dunselman, Cornells van Katwijk, Cornelia R Lo, and Henk J Hoogland. Endometrial wavelike movements during the menstrual cycle. Fertility and sterility, 65(4):746-749, 1996.

[0204] [6] Marga M IJland, Johannes LH Evers, and Henk J Hoogland. Velocity of endometrial wavelike activity in spontaneous cycles. Fertility and sterility, 68(l):72-75, 1997.

[0205] [7] Nienke PM Kuijsters, Federica Sammali, Xin Ye, Celine Blank, Lin Xu, Massimo Mischi, Benedictus C Schoot, and Chiara Rabotti. Propagation of spontaneous electrical activity in the ex vivo human uterus. Pfliigers Archiv-European Journal of Physiology, 472(8): 1065-1078, 2020.

[0206] [8] G Leyendecker, G Kunz, L Wildt, D Beil, and H Deininger. Uterine hyperperistalsis and dysperistalsis as dysfunctions of the mechanism of rapid sperm transport in patients with endometriosis and infertility. Human reproduction, 11(7): 1542- 1551, 1996.

[0207] [9] Bruce D Lucas and Takeo Kanade. An iterative image registration technique with an application to stereo vision. In IJCAI’81: 7th international joint conference on Artificial intelligence, volume 2, pages 674-679, 1981.

[0208]

[0010] Massimo Mischi, Chuan Chen, Tanya Ignatenko, Hinke de Lau, Beijing Ding, SG Guid Oei, and Chiara Rabotti. Dedicated entropy measures for early assessment of pregnancy progression from single-channel electrohysterography. IEEE Transactions on Biomedical Engineering, 65(4):875-884, 2017.

[0209]

[0011] Chiara Rabotti, S Guid Oei, Janneke van’t Hooft, and Massimo Mischi. Electrohysterographic propagation velocity for preterm delivery prediction. American Journal of Obstetrics & Gynecology, 205(6):e9-el0, 2011.

[0210]

[0012] Chiara Rabotti. Characterization of uterine activity by electrohysterography. 2010.

[0211]

[0013] Federica Sammali, Celine Blank, Tom GH Bakkes, Yizhou Huang, Chiara Rabotti, Benedictus C Schoot, and Massimo Mischi. Multi-modal uterine-activity measurements for prediction of embryo implantation by machine learning. IEEE Access, 9:47096-47111, 2021.

[0212]

[0014] Federica Sammali, Nienke Pertronella Maria Kuijsters, Benedictus Christiaan Schoot, Massimo Mischi, and Chiara Rabotti. Feasibility of transabdominal electrohysterography for analysis of uterine activity in nonpregnant women. Reproductive Sciences, 25(7): 1124-1133, 2018.

[0015] Federica Sammali, Nienke PM Kuijsters, Yizhou Huang, Celine Blank, Chiara Rabotti, Benedictus C School, and Massimo Mischi. Dedicated ultrasound speckle tracking for quantitative analysis of uterine motion outside pregnancy. IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 66(3):581-590, 2018.

[0213]

[0016] Apostolos Ziogas, Emmanouil Xydias, Sofia Kalantzi, Despoina Papageorgouli, Polyxeni-Natalia Liasidi, Ioanna Lamari, and Alexandros Daponte. The diagnostic accuracy of 3d ultrasound compared to 2d ultrasound and mri in the assessment of deep myometrial invasion in endometrial cancer patients: A systematic review. Taiwanese Journal of Obstetrics and Gynecology, 61(5):746-754, 2022.

[0214] Furthermore, reference is made to

[0215] [Pl] M. Mischi and B.C. School, “Two-dimensional and three- dimensional strain mapping for uterine contractions,” WO2019053249A1, issued patent in the USA and pending in EU. TU / e: 17-09-2017.

[0216] [P2] M. Mischi and B.C. Schoot, “Quantitative analysis of uterine spatiotemporal motion patterns and Coordination,” WO2021255283A1, pending patent. TU / e: 19-6-2020.

[0217] • Van Sloun, R. J., Demi, L., Postema, A. W., De La Rosette, J. J., Wijkstra, H , & Mischi, M. (2016). Entropy of ultrasound-contrast-agent velocity fields for angiogenesis imaging in prostate cancer. IEEE transactions on medical imaging, 36(3), 826-837.

[0218] • Chen, P., van Sloun, R. J., Turco, S., Wijkstra, H., Filomena, D., Agati, L., ... & Mischi, M. (2021). Blood flow patterns estimation in the left ventricle with low-rate 2D and 3D dynamic contrast-enhanced ultrasound. Computer Methods and Programs in Biomedicine, 198, 105810.

[0219] • Saporito, S., van Assen, H. C., Houthuizen, P., Aben, J. P. M., Strik, M., van Middendorp, L. B , ... & Mischi, M. (2016). Assessment of left ventricular mechanical dyssynchrony in left bundle branch block canine model : Comparison between cine and tagged MRI. Journal of Magnetic Resonance Imaging, 44(4), 956-963.

[0220] • Rabotti, C., Mischi, M., Oei, S. G., & Bergmans, J. W. (2010). Noninvasive estimation of the electrohysterographic action-potential conduction velocity. IEEE Transactions on Biomedical Engineering, 57(9), 2178-2187.

[0221] • Huang, Y , Rees, C , Sammali, F , Blank, C , Schoot, D., & Mischi, M. (2022). Characterization of uterine peristaltic waves by ultrasound strain analysis. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 69(6), 2050-2060. The following numbered clauses define embodiments that have been considered.

[0222] Clause 1. A method (500) for estimating a velocity vector field of a strain wave in a uterus, the method comprising: obtaining (510) a time-dependent strain field indicating an estimate of strain at a plurality of defined locations distributed over a uterus during a time period; performing (520) a spatiotemporal analysis of the obtained strain field to determine a velocity vector field of a strain wave in the uterus, comprising: determining (530) local velocity vectors for the plurality of defined locations, the local velocity vectors together forming the velocity vector field, wherein determining a local velocity vector for a particular defined location comprises: selecting (531) multiple defined locations from the plurality of defined locations in a neighborhood of the particular defined location; estimating (532) delays in the strain wave between pairs of selected defined locations over all or part of the time period; estimating (533) distances between the pairs of selected defined locations; computing (534) the local velocity vector for the particular defined location from the estimated delays and distances; providing (540) the velocity vector field for visualization and / or for characterization of the velocity vector field.

[0223] Clause 2. A method for estimating a velocity vector field as in Clause 1, wherein the velocity vector field is determined for multiple subsequent parts of the time period.

[0224] Clause 3. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the plurality of defined locations are distributed over the subendometrial layer of the uterus.

[0225] Clause 4. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the defined locations are distributed in a grid, the neighborhood of the particular defined location comprises a kernel of defined locations centered on the particular defined location, and / or the neighborhood of the particular defined location comprises defined locations within a threshold distance to the particular defined location,

[0226] Clause 5. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein estimating a delay in the strain wave between a pair of selected defined locations over all or part of the time period comprises maximizing a utility function for a time delay, the utility function comprising a correlation between a strain over time at a first defined location of the pair according to the strain field, and a strain over time shifted over the time delay at a second defined location of the pair.

[0227] Clause 6. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the local velocity vector is obtained by minimizing a loss function for the local velocity vector, the loss function comprising differences between the estimated distances between the pairs of selected defined locations and products of the local velocity vector and the estimated delays.

[0228] Clause 7. A method for estimating a velocity vector field as in clause 5, wherein the loss function comprises weights depending on a signal -to-noise ratio (SNR) of the strain at the selected defined locations.

[0229] Clause 8. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the velocity vector field visualization comprises an overlay over a visualization of the uterus, and / or the velocity vector field visualization comprises a plurality of arrows indicating the velocity vector field, e.g., at the defined locations or a subset thereof, and / or the velocity vector field visualization comprises streamlines for the velocity vector field, and / or the velocity vector field visualization comprises a rendering of a propagation path of contractions in the uterus.

[0230] Clause 9. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the method comprises deriving the strain field relative to the uterine anatomy and orientation, and / or the strain field may be determined in a longitudinal direction defined along a longitudinal axis of the uterus, and / or a transversal direction and / or a circular direction, and / or a longitudinal direction is defined along the endometrial line of the uterus and the transversal or circular directions defined perpendicular to it, strain being computed relative to the longitudinal, transversal and / or radial direction, and / or wherein the strain is directionless, e.g., area strain and / or volume strain

[0231] Clause 10. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein obtaining a time-dependent strain field comprises obtaining 2D or 3D ultrasound measurements of the uterus, having multiple frames, determining an orientation of the uterus, generating a plurality of tracking points inside a region of interest of the uterus, tracking and estimating the displacement of the plurality of tracking points between frames of the ultrasound measurements; calculating, from a varying distance between subsets of the tracking points, at least one of: a transversal strain, a longitudinal strain, a circular strain, an area strain and a volume strain

[0232] Clause 11. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the visualization and / or characterization is used to diagnose uterine dysfunctions, e.g., one or more of: adenomyosis, cesarean niches, fibromas, myomas, and endometriosis, and / or wherein the visualization and / or characterization is used for the prediction of in-vitro-fertilization success.

[0233] Clause 12. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein characterization of the velocity vector field comprises computing a global uterine activity parameter from the velocity vector field quantifying the velocity vector field, e.g., average velocity, average velocity magnitude, and / or average velocity direction, a vorticity, vorticity map, a dyssynchrony index.

[0234] Clause 13. A method for estimating a velocity vector field as in any one of the preceding clauses, wherein the characterization comprises computing an entropy of the velocity vector field globally and / or within a kernel and / or comprising determining a probability distribution of velocity vectors occurring globally and / or within a kernel, computing an entropy of the probability distribution, and / or computing a conditional entropy and / or mutual information between one or more velocity vectors and between one or more neighboring velocity vectors, and / or computing an entropy map and / or mutual information map, and visualizing the map, e.g., as an overlay.

[0235] Clause 14. A method as in any one of the preceding clauses for determining a propagation path, the method comprising identify a point (xn, yn, zn~) as the origin of the propagation path, repeatedly determine a next point (xn+1,yn+1,zn+1) likely to be a next location of propagation within a predefined search area.

[0236] Clause 15. A system comprising one or more processors; and one or more storage devices storing instructions that, when executed by the one or more processors, cause the one or more processors to perform operations for a method according to any one of the preceding clauses.

[0237] Clause 16. A transitory or non-transitory computer storage medium encoded with instructions that, when executed by one or more computers, cause the one or more computers to perform a method according to any one of clauses 1-14.

[0238] Many different ways of executing the method, e.g., method 500, are possible, as will be apparent to a person skilled in the art. For example, the order of the steps can be performed in the shown order, but the order of the steps can be varied or some steps may be executed in parallel. Moreover, in between steps other method steps may be inserted. The inserted steps may represent refinements of the method such as described herein, or may be unrelated to the method. For example, some steps may be executed, at least partially, in parallel. Moreover, a given step may not have finished completely before a next step is started.

[0239] Embodiments of the method may be executed using software, which comprises instructions for causing a processor system to perform an embodiment of method 500. Software may only include those steps taken by a particular sub-entity of the system. The software and / or other data according to an embodiment may be stored in a non-transitory storage medium, such as a hard disk, a floppy, a memory, an optical disc, read only memory, random access memory, CD-ROMs, magnetic tape, optical data storage devices, etc. Transitory signals and carrier waves are excluded from non- transitory media.

[0240] The software may be sent as a transitory signal along a wire, or wireless, e.g., sent as a transitory signal over a data network, e.g., the Internet. For example, signals and / or carrier waves may serve as a transitory medium for carrying information. For example, a modulated electromagnetic wave may carry a signal bearing the software and / or other data according to an embodiment.

[0241] The software may be made available for download and / or for remote usage on a server. Embodiments of the method may be executed using a bitstream arranged to configure programmable logic, e.g., a field-programmable gate array (FPGA), to perform an embodiment of the method.

[0242] It will be appreciated that the presently disclosed subject matter also extends to computer programs, particularly computer programs on or in a carrier, adapted for putting the presently disclosed subject matter into practice. The program may be in the form of source code, object code, a code intermediate source, and object code such as partially compiled form, or in any other form suitable for use in the implementation of an embodiment of the method. An embodiment relating to a computer program product comprises computer executable instructions corresponding to each of the processing steps of at least one of the methods set forth. These instructions may be subdivided into subroutines and / or be stored in one or more files that may be linked statically or dynamically. Another embodiment relating to a computer program product comprises computer executable instructions corresponding to each of the devices, units and / or parts of at least one of the systems and / or products set forth. Figure 10a shows a computer readable medium 1000 having a writable part 1010, and a computer readable medium 1001 also having a writable part. Computer readable medium 1000 is shown in the form of an optically readable medium. Computer readable medium 1001 is shown in the form of an electronic memory, in this case a memory card. Computer readable medium 1000 and 1001 may store data 1020 wherein the data may indicate instructions, which when executed by a processor system, cause a processor system to perform an embodiment of a method for estimating a velocity vector field, according to an embodiment. The computer program 1020 may be embodied on the computer readable medium 1000 as physical marks or by magnetization of the computer readable medium 1000. However, any other suitable embodiment is conceivable as well. Furthermore, it will be appreciated that, although the computer readable medium 1000 is shown here as an optical disc, the computer readable medium 1000 may be any suitable computer readable medium, such as a hard disk, solid state memory, flash memory, etc., and may be non-recordable or recordable. The computer program 1020 comprises instructions for causing a processor system to perform an embodiment of said method for estimating a velocity vector field.

[0243] Figure 10b shows in a schematic representation of a processor system 1140 according to an embodiment. The processor system comprises one or more integrated circuits 1110. The architecture of the one or more integrated circuits 1110 is schematically shown in Figure 10b. Circuit 1110 comprises a processing unit 1120, e.g., a CPU, for running computer program components to execute a method according to an embodiment and / or implement its modules or units. Circuit 1110 comprises a memory 1122 for storing programming code, data, etc. Part of memory 1122 may be read-only. Circuit 1110 may comprise a communication element 1126, e.g., an antenna, connectors or both, and the like. Circuit 1110 may comprise a dedicated integrated circuit 1124 for performing part or all of the processing defined in the method. Processor 1120, memory 1122, dedicated IC 1124 and communication element 1126 may be connected to each other via an interconnect 1130, say a bus. The processor system 1140 may be arranged for contact and / or contact-less communication, using an antenna and / or connectors, respectively.

[0244] For example, in an embodiment, processor system 1140, e.g., a velocity vector field estimating device may comprise a processor circuit and a memory circuit, the processor being arranged to execute software stored in the memory circuit. For example, the processor circuit may be an Intel Core i7 processor, ARM Cortex-R8, etc. The memory circuit may be an ROM circuit, or a non-volatile memory, e.g., a flash memory. The memory circuit may be a volatile memory, e.g., an SRAM memory. In the latter case, the device may comprise a non-volatile software interface, e.g., a hard drive, a network interface, etc., arranged for providing the software.

[0245] While system 1140 is shown as including one of each described component, the various components may be duplicated in various embodiments. For example, the processing unit 1120 may include multiple microprocessors that are configured to independently execute the methods described herein or are configured to perform elements or subroutines of the methods described herein such that the multiple processors cooperate to achieve the functionality described herein. Further, where the system 1140 is implemented in a cloud computing system, the various hardware components may belong to separate physical systems. For example, the processor 1120 may include a first processor in a first server and a second processor in a second server.

[0246] It should be noted that the above-mentioned embodiments illustrate rather than limit the presently disclosed subject matter, and that those skilled in the art will be able to design many alternative embodiments.

[0247] In the claims, any reference signs placed between parentheses shall not be construed as limiting the claim. Use of the verb ‘comprise’ and its conjugations does not exclude the presence of elements or steps other than those stated in a claim. The article ‘a’ or ‘an’ preceding an element does not exclude the presence of a plurality of such elements. Expressions such as “at least one of’ when preceding a list of elements represent a selection of all or of any subset of elements from the list. For example, the expression, “at least one of A, B, and C” should be understood as including only A, only B, only C, both A and B, both A and C, both B and C, or all of A, B, and C. The presently disclosed subject matter may be implemented by hardware comprising several distinct elements, and by a suitably programmed computer. In the device claim enumerating several parts, several of these parts may be embodied by one and the same item of hardware. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measures cannot be used to advantage.

[0248] In the claims references in parentheses refer to reference signs in drawings of exemplifying embodiments or to formulas of embodiments, thus increasing the intelligibility of the claim. These references shall not be construed as limiting the claim.

Claims

CLAIMSClaim 1. A method (500) for estimating a velocity vector field of a strain wave in a uterus, the method comprising: obtaining (510) a time-dependent strain field indicating an estimate of strain at a plurality of defined locations distributed over a uterus during a time period; performing (520) a spatiotemporal analysis of the obtained strain field to determine a velocity vector field of a strain wave in the uterus, comprising: determining (530) local velocity vectors for the plurality of defined locations, the local velocity vectors together forming the velocity vector field, wherein determining a local velocity vector for a particular defined location comprises: selecting (531) multiple defined locations from the plurality of defined locations in a neighborhood of the particular defined location; estimating (532) delays in the strain wave between pairs of selected defined locations over all or part of the time period; estimating (533) distances between the pairs of selected defined locations; computing (534) the local velocity vector for the particular defined location from the estimated delays and distances; providing (540) the velocity vector field for visualization and / or for characterization of the velocity vector field.Claim 2. A method for estimating a velocity vector field as in Claim 1, wherein the velocity vector field is determined for multiple subsequent parts of the time period.Claim 3. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the plurality of defined locations are distributed over the subendometrial layer of the uterus.Claim 4. A method for estimating a velocity vector field as in any one of the preceding claims, whereinthe defined locations are distributed in a grid, the neighborhood of the particular defined location comprises a kernel of defined locations centered on the particular defined location, and / or the neighborhood of the particular defined location comprises defined locations within a threshold distance to the particular defined location,Claim 5. A method for estimating a velocity vector field as in any one of the preceding claims, wherein estimating a delay in the strain wave between a pair of selected defined locations over all or part of the time period comprises maximizing a utility function for a time delay, the utility function comprising a correlation between a strain over time at a first defined location of the pair according to the strain field, and a strain over time shifted over the time delay at a second defined location of the pair.Claim 6. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the local velocity vector is obtained by minimizing a loss function for the local velocity vector, the loss function comprising differences between the estimated distances between the pairs of selected defined locations and products of the local velocity vector and the estimated delays.Claim 7. A method for estimating a velocity vector field as in claim 5, wherein the loss function comprises weights depending on a signal -to-noise ratio (SNR) of the strain at the selected defined locations.Claim 8. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the velocity vector field visualization comprises an overlay over a visualization of the uterus, and / or the velocity vector field visualization comprises a plurality of arrows indicating the velocity vector field, e.g., at the defined locations or a subset thereof, and / or the velocity vector field visualization comprises streamlines for the velocity vector field, and / or the velocity vector field visualization comprises a rendering of a propagation path of contractions in the uterus.Claim 9. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the method comprises deriving the strain field relative to the uterine anatomy and orientation, and / or the strain field may be determined in a longitudinal direction defined along a longitudinal axis of the uterus, and / or a transversal direction and / or a circular direction, and / or a longitudinal direction is defined along the endometrial line of the uterus and the transversal or circular directions defined perpendicular to it, strain being computed relative to the longitudinal, transversal and / or radial direction, and / or wherein the strain is directionless, e.g., area strain and / or volume strainClaim 10. A method for estimating a velocity vector field as in any one of the preceding claims, wherein obtaining a time-dependent strain field comprises obtaining 2D or 3D ultrasound measurements of the uterus, having multiple frames, determining an orientation of the uterus, generating a plurality of tracking points inside a region of interest of the uterus, tracking and estimating the displacement of the plurality of tracking points between frames of the ultrasound measurements; calculating, from a varying distance between subsets of the tracking points, at least one of: a transversal strain, a longitudinal strain, a circular strain, an area strain and a volume strainClaim 11. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the visualization and / or characterization is used to diagnose uterine dysfunctions, e.g., one or more of: adenomyosis, cesarean niches, fibromas, myomas, and endometriosis, and / or wherein the visualization and / or characterization is used for the prediction of in-vitro-fertilization success.Claim 12. A method for estimating a velocity vector field as in any one of the preceding claims, wherein characterization of the velocity vector field comprises computing a globaluterine activity parameter from the velocity vector field quantifying the velocity vector field, e g., average velocity, average velocity magnitude, and / or average velocity direction, a vorticity, vorticity map, a dyssynchrony index.Claim 13. A method for estimating a velocity vector field as in any one of the preceding claims, wherein the characterization comprises computing an entropy of the velocity vector field globally and / or within a kernel and / or comprising determining a probability distribution of velocity vectors occurring globally and / or within a kernel, computing an entropy of the probability distribution, and / or computing a conditional entropy and / or mutual information between one or more velocity vectors and between one or more neighboring velocity vectors, and / or computing an entropy map and / or mutual information map, and visualizing the map, e.g., as an overlay.Claim 14. A method as in any one of the preceding claims for determining a propagation path, the method comprising identify a point (xn, yn, zn~) as the origin of the propagation path, repeatedly determine a next point (xn+1,yn+1, zn+1) likely to be a next location of propagation within a predefined search area.Claim 15. A system comprising one or more processors; and one or more storage devices storing instructions that, when executed by the one or more processors, cause the one or more processors to perform operations for a method according to any one of the preceding claims.Claim 16. A transitory or non-transitory computer storage medium encoded with instructions that, when executed by one or more computers, cause the one or more computers to perform a method according to any one of claims 1-14.

Images