Enabling Augmented Segmentation and Registration in Ultrasound-Guided Spinal Surgery via Realistic Ultrasound Synthesis from Diagnostic CT Volume

Since the goal of the proposed US simulation system is to train doctors or artificial intelligence agents for surgical robots, the moving space of probe motion cannot be the entire 3D space. The probe in the system can move in the provided scanning space, but the probe is a curved surface, it is the probe surface, and the scanning space will not be fitted perfectly. The scanning space will wrap around and fit the probe if the pressure is simulated. However, the CT volume will not deform with the surface. Hence, an algorithm needs to be proposed to make the CT data deform with the probe surface along the scanning space. Using the spring-mass model to simulate deformations in scanning space by using the moving least squares algorithm to fit the CT deformation is certainly good, but also has two demerits. The first is solving the spring-mass simulation and the moving least squares equation will consume a lot of computational resources. Secondly, the algorithm requires a lot of extra work to mark the anchor points on the images. Under this situation, we changed the local translation image warping algorithm proposed by A. Gustafsson for CT data and apply it with the shape of the probe to the image. The algorithm only needs to know the shape of the probe and the HU value of the CT data of that slice to simulate the motion of the tissues. If further acceleration is desired, the weight value of the HU parameter can be set to 1, then the UV coordinate mapping only needs to be calculated once, without affecting the imaging efficiency at all.

The equation is formulated as Equ. 1

In the above equation, $f$ controls the ratio of deformation. In our case, we want to push the intersection point of the center-line of the probe and the top line of the sampled image, as the blue dot is shown in Fig.4, to the top of the probe which is shown in the figure as the red dot. According to the mentioned algorithm, all the tissue between the blue and red dots needs to be pressed downwards. Assuming that the tissue is a rigid body, the thickness of the tissue does not change by pressing with a very big force. In this case, the original position at the red dot should be pushed to the top of the green curve. Then a green arc that represents the limit of tissue movement can be drawn. The squeezing of the tissue in the area between the probe arc and the green arc will cause high reflections or absorption in this area. Hence, in the process of simulating transmission, it is significant to make the sound waves attenuate less in this area, otherwise, the squeezed tissue generated by the algorithm will completely block the transmission of US. As you can see from Fig. 3, there are a few bright lines below the probe curve, which are the signal of the squeezed tissue and it is identical to the vivo US image.

As mentioned above, the probe cannot move in space arbitrarily. To automate the acquired US images, two types of probe moving regulation are designed in this paper using spinal US as an instance. Define $M$ as a mesh in 3d space, clipping a 3D sub-mesh $M_{r}$ of interest by the bounding box of the scanning target, the movement of the probe is only meaningful in this manifold. The first way is to move freely on the grid, for any point $p$ on $M_{r}$, the direction of movement of p will be curved along the grid normal vector. In this case, a very small $dx$ or $dy$ will not make a great $dz$, which would not make the image change dramatically and could not acquire a continuous image. The second method is to slice $M_{r}$ into a series of curves in 3D space along the forward direction $M_{r0}...M_{rn}$. In each polyline, the points can only move forward or backward along the tangent direction by a certain distance. It has the merit that keeping the components of the probe’s normal vector remain consistent across that curve, which is often used for automatic circular or flat scanning of a target and ensuring that the target being swept is on a line or a point. Besides, both of these movement methods can be employed by a reinforcement learning-based agent, and the first method is more suitable for continuous control, while the second method is more suitable for discrete control.

The first step in the simulation is to map the CT-HU image into an acoustic impedance image, with the mapping look-up table references to [21]. In the mapping table, each HU value will be mapped to an acoustic impedance Z. The image of this acoustic impedance will then be processed instead of the CT image. Based on the acoustic physical transfer properties, we first introduce the Fresnel equation to calculate the reflection. Dividing sound waves into those parallel and those perpendicular to the plane, the equations are:

When sound waves transmitted from one medium to another, it is not only be reflected but also be refracted. This refracted has the same physical properties as the refracted of light passing through different media which is given by Snell’s law:

In fact, when a sound wave is refracted or bent, the direction will be kept until the next transmission. It is possible to simulate multi-level refraction and reflection values using a recursive ray-tracing algorithm However, since the amplitude of the refracted sound waves is small and has little effect on the generated image quality, the proposed method does not consider its direction after refraction and only superimposes the amplitude on the incident sound wave of the next medium.

Irradiance: Lambert’s cosine theorem is introduced to take the irradiance of the sound wave should into consideration. This equation means that the more perpendicular the direction of the tissue gradient and the sound wave, the stronger the sound reflection and the brighter the tissue.

Attenuation in the propagation of US passing through tissue could be broadly divided into reflection, scattering, and absorption, where absorption accounts for the majority of it. The paper [22]states that US absorption of substances such as bovine and porcine livers accounts for $90\%$ or even $100\%$ of US propagation attenuation. Even in the tissues that exhibit anisotropy such as the bovine brain and leg muscle absorption still makes a major contribution to the attenuation. Therefore, mapping the CT scan to an absorption image is important. The mapping table is obtained by interpolating the dataI in the paper [23]. The formula for absorption is given by the Lambert-Beer law. Since different tissues may not have the same US absorption capacity in the same Hu value, in order to get a better result from the generated image, an adjustable parameter $\alpha$ is added to the original equation. The modified equation is as follows:

The propagation of the US is the most time-consuming part of the generation process. Particularly, the steps of finding receiver tissues of reflection and refraction would take a lot of time in detecting intersections that occur frequently.

At the same time, the direction of the adjacent pixels is discrete, whereas the direction of acoustic transmission is continuous. So it may happen that the line segments intersect while the pixels do not have an intersection. Each pixel among the eight neighborhoods should be regarded as the obstacle of the wave propagation and the normalized result of the cosine value of the direction of propagation and the direction of the pixel connection is used as the weight of the obstacle. That is, if the connection direction is in the opposite direction to the direction of propagation, the weights will be zero and the obstacle weight is maximized when the directions are exactly the same.

In order to calculate the amplitude enhancement due to reflection and refraction easily, each point also calculates its wave propagation to the neighborhood. In this case, a negative cosine value of the pixel connection direction and the propagation direction will consider as reflection. Cosine values less than 0.8 (a hyper-parameter) will be considered as contrast enhancement by refraction and scattering. The schematic diagram of propagation is shown in 1

When the transducer receives the reflected sound waves from the medium gap, a bright stripe appears in the US image. Because the reflection only occurs when the properties of the medium change, there should be only one reflection between the bone surface and the muscle tissue, but the US image often shows a very thick bright stripe between bone and muscles. hacihaliloglu et al. [24] propose that this bright stripe is produced by the thickness of the US in the elevational direction. Based on this theory, this paper sampled 3d gradient textures in front of and behind the current plane (As the blue line in figure 6). These textures are multiplied by a weight value $\alpha$ which is inversely proportional to the distance and then accumulated into the current generated US image. Simulating thick stripes is significant for neural network training; if the thick stripes are not present in the training data, the output of the model will deviate from the ground truth when segmenting the bone surface on the real US data.

In summary, the generated US image is obtained by calculating the propagation image with the reflection image and the absorption image and then blending these three images together with the radial noise.

To cover the global information, MobileViT block first utilizes a $3\times 3$ convolution to aggregate the neighboring features. The pixels are then grouped into patches, each patch contains $h\times w$ pixels in which $h$ and $w$ represent the height and the width of the patch. The $i_{th}$ token of each patch transverses its feature from each other by SA. When $h<=3\wedge w<=3$, each pixel can access information from all pixels. Details of the MobileViT block could be found in [25].

We find that the region of interest (RoI) in US images has similar textures, but is weak in amplitude to other tissues. Basic $3\times 3$ convolutions are capable of modeling the neighboring texture but fail to learn long-range textures. As the local texture of the to-be-segmented region is similar to the other tissues, we have to design a module that covers the texture within and outside the to-be-segmented region.

In vision transformers, long-range texture modeling is achieved by self-attention (SA) of all tokens as in Equ. 6:

in which $x\in\mathbb{R}^{N\times d}$ represents the $N$ input tokens with dimension $d$, $W_{q}\in\mathbb{R}^{d\times d_{q}}$, $W_{k}\in\mathbb{R}^{d\times d_{q}}$, and $W_{v}\in\mathbb{R}^{d\times d_{v}}$ are learnable parameters. We rewrite the self-attention to the formula as in Equ. 7:

where $\mathcal{A}_{i}$ represents the set of accessible tokens of token $x_{i}$. As $W_{q}$ and $W_{k}$ only attribute in modeling the relationship of queries, once $W_{v}$ is fixed, the space generated by $x$ is therefore limited. Meanwhile, as shown in Fig. 8(a), SA tends to aggregate “similar” features from the tokens, which helps the model to learn long-range dependency. However, because of the nature of US images, long-range contrast is also important to be learned, but SA fails to do so.

Convolution operation could be formulated in the same formulation as SA, as shown in Equ. 8:

In this case, if $rank(x_{j}W_{v})=d$, the space generated by $x$ is $\mathbb{R}^{d}$. We randomize $1000$ different pairs of $W_{q}$ and $W_{v}$, $1000$ different $\lambda$ and fixed $W_{v}$. The outputs are shown in Fig. 8. As demonstrated in the figure, the outputs of convolution fulfill the space, while the outputs of SA are gathered.

A recent work [26] also demonstrates that CNNs with super large kernels learn feature representations from the super large receptive field (sub-global), and achieve comparable performance with state-of-the-art transformers with fewer parameters.

As aforementioned, we need a module to gather long-range textures. As SA fails to learn long-range contrast, and a super-large convolution kernel leads to large computational cost, we need to design a light yet effective convolution module for long-range contrast learning, namely Long-range Contrast Learning Module (LCLM).

To cover long-range contrast with fewer parameters, a simple way is to utilize dilation convolution. In this work, we adopt 3 cascaded dilation convolution layers with dilation {3,5,11}. Along with one basic convolution layer which appears at the end of the MobileViT block, LCLM covers up to $41\times 41$ pixels densely. We show the receptive field in Fig. 9. Each convolution layer is followed by a normalization layer and an activation layer. A depth-wise LCLM only needs $3\times 3\times 3\times d$ parameters, which is even smaller than a standard convolution layer ($3\times 3\times d\times d$, $d>>3$).

To validate the proposed US simulation method and the lightweight segmentation model, we train the proposed model on synthetic US images and inference on real US images. We set the batch size to 32, utilize the Adam optimizer, and set the learning rate to 1e-4.

We first compare the proposed model with other US segmentation models in Tab. II. All performances reported in the table are re-implemented by ourselves. Compared with UNet-based models, the proposed model is much smaller and more effective. Compared with MViT-FPN, the proposed approach achieves much better performance.

To utilize the model, in this case, the segmentation is to achieve better registration, therefore we remove the cases with obvious segmentation errors to obtain our (selected) result, which is used in the subsequent registration. In this case, the effect of the network trained with generated data can be equivalent to or better than the effect of the SOTA model trained by a large number of real US images with an accurately labeled. Our label is acquired by ct segmentation, which is more challenging for the convergence of lightweight neural networks since he can label bone surfaces which is invisible or ambiguous in some situation in US.

To demonstrate the effectiveness of synthetic US over initial CT scans, long-range dependency, and long-range contrast, we validate the aforementioned settings and show the results in Tab. III. Training the model with CT scans leads to a dramatic IoU decrease, which demonstrates that the proposed method can effectively synthesize US images. Without LCLM, the model fails to learn the long-range contrast of US images and results in a performance decrease. Without MViT, the IoU slightly drops, which represents that the model benefits from long-range dependency, but is less significant compared with long-range contrast.

The US simulation system is proposed to better obtain labeled data from multiple poses, patients, and environments, enabling US-based surgical navigation to be accomplished. In this paper, we validated whether the output of the model trained on synthetic data can align the preoperative CT and the planned trajectory into the intra-operative phase, and the spinal bone surface registration for pedicle screw placement is the instance.

The validation is divided into two parts, the first is to validate the accuracy of the rigid body registration, i.e. translation and rotation. The second part is to validate whether the pre-operatively planned trajectory with the aforementioned transforms will cause complications for the surgery or not, i.e. whether the intra-operative trajectory will touch vital organs.

The experiments were performed on three phantoms, a spine phantom in water, a 3D-printed human spine in water, and a bovine spine in agar gel. The registration steps were divided into a coarse alignment and an individual registration for each segment. The coarse registration used all point clouds to calculate the approximate position and orientation, and the final registration used a de-noised Iterative Closest Point(ICP) method. The results of the registration are shown in the table below.

In the registration phase, the model inputs the US image outputs the segmentation label, and extracts the upper contour of its maximum connectivity, hence the point cloud is created with the calibration matrix of the US probe. The calibration of the US is done by two cross-wire phantom [27], with 0.97mm MSE(min square loss) loss. Since the coordinate system of each segment does not coincide, the rotation evaluation in angular will be represented by the rotation of pedicle screws in the next part. To evaluate the error of registration, the point-to-point MSE is used, which is also the objective function of the ICP algorithm. The result of registration is shown in Fig. 15, which shows the ground truth in the intra-operative space(white model) and the results of registration with US (green model). In summary, the segmentation algorithm trained by the generated US image works in the registration with the error $0.133\sim 3.366$mm. The maximum appears in the z-axis of L5, which is around 3 mm, which is probably caused by the shape of L5 being more different from the shape of other lumbar.

To verify the registration effectiveness of the system in the target tasks, we compared the translation and rotation of the tip of the screw with the US point cloud registration and ground truth using an expert-planned intra-operative plan for pedicle screw placement. The rotation is determined by the angle between two pedicle axis vectors in 3D space which is expressed in degrees. As shown in Fig. 14, the error in translation for the tip of the screw ranges from 0.027 mm to 4.031 mm, with the maximum value still occurring on the z-axis of L5, which is consistent with the above assessment. The error in the rotation is between 3.05 degrees and 4.47 degrees, therefore, it can be concluded that the pre-trained model with generated US image can perform the alignment with small errors. The results of the pedicle screw placement are shown in figure 16, which shows this registration does not put the patient at risk of complications.