Separation of perfusion phases in angiographies

1.   Introduction

Cerebral Angiography is an invasive diagnostic and therapeutic tool for neurovascular imaging. It is based on arterial vascular opacification using iodine contrast agents. It allows the neurointerventionnal radiologists to navigate along the arterial tree using catheters to diagnose and treat neurovascular diseases. It presents specificities related to acquisition and interpretation modalities.

The first specificity is the dynamical acquisition over time according to 2 planes in 2 dimensions. While 2D acquisition can be a source of potentially misleading superimpositions and constructed images, the use of two orthogonal planes reduces interpretation errors. The main advantage of this temporal dynamic acquisition in arteriography is the real-time evaluation of the evolution of the contrast agent within the traveled vascular structures: arteries, arterioles, capillaries, parenchyma, venules and then veins. From this evaluation, important therapeutic decisions are made in case of abnormalities: mechanical thrombectomy, preventive embolization, coiling, in situ thrombolysis, angioplasty, etc.

The second specificity is related to the interpretation of the images. This is carried out in immediate procedure. The main interpretation lies in the large number of anatomical variants of the cerebral arterial and venous vascularization. In the majority of cases, this is based on a qualitative or semi-quantitative dynamic reading procedure of the images. The qualitative analysis is based on the presence of stenosis, occlusion, outward bulging, the presence of early venous drainage or delayed arterial opacification in the venous phase. The exploration of the parenchymographic phase is a fundamental but often neglected step in cerebral arteriography. This phase represent the cerebral vascular dynamics counterpart. Yet, delayed parenchymogram will constitute the turning point for the radiologist when looking for vascular abnormality. Last, the parenchymogram is not subject to variation from one healthy subject to another. Such is not the case for the other angiographic phases.

In practice, semi-quantitative analysis can be carried out based on visual diagnostic scales: Capillary Index Score for predicting tissue viability before thrombectomy, TICI Score for vascular tree re- constructions after thrombectomy, Raymond Scale for aneurysmal recanalization etc... Until today, the only validated quantitative measure is the measurement the size of stenosis. Other studies have focused on quantitatively evaluating the flow within the structures of interest by studying perfusion in particular. Unfortunately, theses methods have not being used in clinical routine. Thus, despite the therapeutic consequences related to the interpretation of arteriographic images, interventional radiologists do not benefit from reliable quantitative evaluation tools. We will distinguish three sets of images corresponding to the decomposition of the angiogram in three phases (Arterial, Parenchymo- gram and Venous phase). These three phases rely on blood brain specific physiology : limited blood brain barrier permeability ^[3], polymorphic arterial vasculature and homogenized parenchymogram. The three phases decomposition were designed to fit the best to the radiologist reading method with the development of the dynamic three step reading grid.

Indeed, the angiogram reading relies on heuristic dynamic three-step analysis. The first step is the arterial phases. While analysis the Willis circle, the radiologists try to identify a slowed down or abnormally present artery. Due to the numerous anatomical variations ^[4], the reading is made simple by identifying 2 mains arteries: anterior cerebral artery and middle cerebral artery. From then, the collaterals and distals bifurcations are analyzed looking for lack of opacification and/or and hyperdense artery. Density of collaterals arteries is also analyzed to have an overview of the blood supply in the ischemic area. The second step is the parenchymogram (or capillary) reading. Whereas many anatomical variations have been described regarding the Willis circle, parenchymogram phase is a reproducible, reliable tool to identify abnormal brain perfusion ^[1]. The abnormalities are characterized by a lack of enhancement in the pathological areas. As this abnormal opacification can be multifactorial, radiologists usually review the arterial phase in the light of the parenchymogram phase.

The third step is the venous phase. Hemorrhagic changes and alternative diagnosis can be depicted at this phase. Since the venous drainage is inconsistent between subject, it provides less information than the two previous phases. Nevertheless, some slowed down arteries can still be seen at this late phase, highly suggestive for occlusion. The third step is the venous phase. Hemorrhagic changes and alternative diagnosis can be depicted at this phase. Since the venous drainage is inconsistent between subject, it provides less information than the two previous phases. Nevertheless, some slowed down arteries can still be seen at this late phase, highly suggestive for occlusion. Frequently, due to superposition and artifacts, the angiogram reading cannot be made this way. It relies on an analysis of the phases altogether mixed and therefore could lead to misinterpretation or discrepancies. A possibility to decompose the three phases in three set of images provide homogenized data cleaned from artifact and superposition. Therefore, it facilitates the 3 steps reading grid. Nevertheless, the utility has to be proven and mandate a clinical trial.

From a technical point of view, we will discriminate the three aforementioned phases by scale and time of occurrence. Indeed the arterial phase is characterized by relatively small scale features (blood vessels) in the image, and occurs at early stages; the parenchymogram phase is a large scale phenomenon (whole organ) and appears at intermediate times; and finally the venous phase is a small scale (blood vessel) phenomenon appearing at later stage of the perfusion angiography. In this paper we will use a redundant base for the representation of the image that is based on a multi-scale decomposition. This representation of the image is based on the Isotropic Undecimated Wavelet Transform (IUWT) which was proposed in ^[10,11] and used in astronomy where it is labeled "starlet transform”, and in medicine for then analysis of retinal images ^[2]. One of the difficulties here is to separate positive contributions to the initial image with the sparsest set of coefficients. In order to achieve this we use a filter bank designed so that synthesis from positive coefficients produces a positive image (see ^[10]), and we apply a basis pursuit algorithm (see ^[7,10]).

The article is organised as follows: first, in section 2 we give a description of the method we have used for the separation of scales in original images, then in section 3 we show how this method has allowed us to obtain a proper parenchymographic and arterial/venous components on real clinical data. In section 3 we will first give an description of the origin of the data along with a primary medical analysis of the cases, then we will describe the pre-treatment pipeline that is applied on these sequences of images before using the scale separation procedure. Finally, a later medical analysis of the images is performed by using the information provided by the parenchymogram and the arterial/venous component.

2.   Separation of space scales

In this section we focus on the separation of the different scales present in the images by the means of the isotropic undecimated wavelet transform (IUWT). The IUWT is a popular method for the analysis, denoising and the compression of images. Initially developed for astronomy (see ^[11]), it relies on the same idea of multiscale decomposition as wavelets, but lacks the decimation step (elimination of even or odd coefficients) classically applied between each step of the decomposition. One consequence of this is the redundancy of the representation of the image through its coefficients: an image of $N$ pixels will be represented with $N(3M+1)$ coefficients, where $M$ is the number of levels in the decomposition.

Let $Y_{i, j}$ for $i \in \left[\kern-0.15em\left[{0, N_1-1} \right]\kern-0.15em\right]$ and $j \in \left[\kern-0.15em\left[{0, N_2-1} \right]\kern-0.15em\right]$ be the (real) greyscale values of an image of size $N_1 \times N_2$, let $f_i$ for $i \in \left[\kern-0.15em\left[{-a, a} \right]\kern-0.15em\right]$ be the real coefficients of a filter kernel. Let us define the vertical and horizontal convolution operators:

Let us define the 2D separable kernel convolution operator:

Finally, in order to compute the Undecimated Wavelet Transform (UWT) of an image with Mallat's a trous algorithm ^[6,9], we introduce the $n$-th level filter $f^{(n)}$:

and $f^{(n)}_i = 0$ for indexes $i$ that are not a multiple of $2^n$. Given two analysis filter banks $h$ and $g$, the $M$-level UWT of an image $Y$ is defined by the coefficients $\alpha = (\alpha^{\text{a}}_{(M)}, \, \alpha^{\text{h, v, d}}_{(M)}, \, \alpha^{\text{h, v, d}}_{(M-1)}, …, \alpha^{\text{h, v, d}}_{(1)})$ which can be obtained with the following recursive definition: set $\alpha^{\text{a}}_{(0)} = Y$ and then, for $n \in \left[\kern-0.15em\left[{0, M-1} \right]\kern-0.15em\right]$,

where $\bar{h}$ and $\bar{g}$ denote the "reversed" filters defined by $\bar{h}_k = h_{-k}$ (r.p. for $g$). The superscripts a, h, v, and d in the above description refer respectively to approximation, horizontal, vertical and diagonal sub-bands (see Figure 1).

Let $\tilde{h}$ and $\tilde{g}$ be two reconstruction filters, the original image $Y$ can be recovered from $\alpha$ with:

The original image is then recovered with: $Y = Y^{(0)}$. The filter bank $(h, g, \tilde{h}, \tilde{g})$ has to satisfy the following condition to achieve perfect reconstruction with the above equations (see ^[7] for details):

where $H$, $\tilde{H}$, $G$, $\tilde{G}$ are the z-transforms of $h$, $\tilde{h}$, $g$, $\tilde{g}$, respectively.

We define $\Phi : \alpha \rightarrow Y$ as defined in (2.4)-(2.7), and $\mathcal{W} : Y \rightarrow \alpha$ as defined in (2.8). Clearly, from the reconstruction condition, $\mathcal{W}$ is a right-inverse of $\Phi$, hence the following mapping $\Pi$ defines a projection onto the set $\{ \alpha \in {\mathbb{R}}^{3NM+1} \; \; | \; \; \Phi \, \alpha = Y \}$:

Using $\Phi$ as an overcomplete basis, our goal is to find the sparsest set of coefficients $\alpha$ that recovers the original image i.e. $\Phi \alpha = Y$. One difficulty is to ensure the positivity of the image that is reconstructed from coefficients at each scale. In ^[10], J.-L. Starck et. al. have proposed a symmetric filter bank satisfying the property (2.9), for which the reconstruction filters are positive:

where the above greyed superscripts denote the array indexes. With the above synthesis filters, if the coefficients $\alpha$ are all positive, the image obtained from any subset of coefficient will be positive, which is particularly useful when separating scales. Our goal is hence to solve the basis pursuit problem:

which can be recast as:

where $\iota_{\alpha \geq 0}$ is the indicator of the cone of positive values of $\alpha$. Denoting $J(\alpha) = |\alpha|_1 + \iota_{\alpha \geq 0}$, we have that (see ref for details) :

The (projected) iterative soft thresholding algorithm 1 proposed in ^[12] aims at solving the above problem. It is based on the proximal algorithm with a projection step for the constraint $\Phi \alpha = Y$.

The projected iterative threshold algorithm proposed in (ref)

Note that in this algorithm the given projection is not an orthogonal projection since $\mathcal{W}$ is a right inverse of $\Phi$ (ensured by the condition (2.9)) but not its pseudo-inverse. This would require $\mathcal{W} \Phi$ to be Hermitian, which in terms of the filter bank writes :

which cannot be achieved with a set of symmetric FIR filters (see appendix for the proof in one dimension of space).

3.   Results

In the following section we show the result of the method described in the previous section on the separation of the different morphological scales involved in the perfusion (i.e. arteries, veins and the parenchyma).

3.1.   Data acquisition

All anonymized data were acquired from anonymized patients (Patient 1 and Patient 2) referring for at the emergency department for cerebral thrombus on the middle cerebral artery requiring mechanical thrombectomy. For the 2 anonymized patient, informed consent was obtained, and study was approved by the local ethics committee. Patient 1 is a 50 y.o male patient, referring for right hemicorporeal motor impairment. Patient 2 is a 67 y.o female referring for left hemicorporeal paresthesia.

For each patient, a cerebral angiography was performed using the same methodology: A Newton TERUMO catheter (n RFEH15010M), 5 French, was inserted through the femoral artery to the internal carotid artery and the radiological study was performed on a conventional antero-posterior view and left-right view with an Axium - Artis biplane sensor (Siemens Healthineers, Erlagen, Germany). The acquisition parameters were a voltage of 76kV and an amperage of 120mA, collimation of 4.8 cm, and a rate of 3 images per second). The contrast agent (Iomeron 300mg/ml, Iomeprol, BRACCO, Italy), was injected into the internal carotid artery at a rate of 4 cc/sec, at a dose of 8 mL, using a MEDRAD automatic injector (Mark V pro VIS, Bayer Healthcare, USA). Contrast agent are used to opacified the brain vasculature: injection is made in arterial phase to have access to the arterial tree. The, the agent process into the capillaries and reflects the parenchymal phase. Finally, the venous drainage is visible.

An X-Ray acquisition was performed before and after mechanical thrombectomy. Dynamic X-rays acquisitions were performed over time: 2 images per sec during 20 seconds after intra arterial injection resulting in a set of 40 images from arterial phase, parenchymal phase and then venous phase. The data were then extracted using Maincare PACS station in DICOM files and then processed using our algorithms.

3.2.   Primary medical analysis of the images

The primary analysis was performed on real time reading. The site of occlusion was determined using heuristic method based on the defect from expected arterial vasculature. Then the parenchymal phase was analyzed to assess abnormal hypodensity in a brain area expected to be opacified. Latter opacification was also recorded as an indirect sign of hypoperfusion. Finally, the venous phase was used to assess delayed venous drainage as another indirect sign for occlusion.

After removing the thrombus from the cerebral middle artery using mechanical thrombectomy, the same 3 step analysis was performed to assess complete filling of all of the expected vascular territory. The speed of the filling was also qualitatively evaluated in order to evaluate the vascular restitution and its consequences on brain perfusion.

In the two patients, same findings were identified. Before mechanical thrombectomy, arterial vasculature was impaired at the site of occlusion (proximal middle cerebral artery). No parenchymal phase nor venous opacification were visible in the corresponding brain territory. After mechanical thrombectomy, arterial vasculature was restituted and both parenchymal phase and venous drainage were identified. Filling speed and venous drainage were similar as expected.

The results of the mechanical thrombectomy was scored using an established reading grid: the TICI score. This qualitative and semi quantitative evaluation was performed by the neurointerventionnal radiologist.

3.3.   Pre-treatment

A pre-treatment is applied to these sequences of images in order to obtain proper opacification maps. A summary of this procedure is described in Figure 2. One crucial step is is the rigid alignment of the time sequences in order to compensate for the patient's motion during the aquisition. This rigid registration is obtained from a modified version of the images that show only the skull contour obtained from a simple thresholding of the initial images. The alignement procedure is based on the phase-correlation method, involving the Fourier-Mellin transfom for the estimation of the rotation. This phase-corellation pipeline has been developped in ^[8] and, once refined, allows for a sub-pixel alignement of the images. We follow the same procedures as in ^[8] to estimate the alignement parameters (translation, rotation and scale) on the skull contours. The initial sequence of images is then aligned and a subtraction with the initial image is applied in order to reveal the opacifiation of the valcular structures only.

While the rigid alignment compensates for the patients motion in the plane of aquisition, it cannot take into account any rotation in a plane orthogonal to the plane of aquisition. Hence some undesired artifacts remain in the final images that reveal anatomical features of the skull. After subtraction, only positive values are supposed to appear since no pixel should appear brighter during perfusion than before perfusion. Finally an affine transform is applied to the pixel values in order to obtain maps of opacification rate instead of raw subtracted values of illumination.

The Figure 3 shows a sample from a sequence of 82 images obtained by the procedure described above (patient 1, frontal view, post-thrombolysis). We can clearly see the different perfusion components: the first image essentially reveals the arteries, the opacification of the parenchyma (capillary flow) can be seen on the second image, and the veins are shown in the third image. Our goal for the rest of the paper is to use the mathematical methods introduced in section 2 in order to separate these three perfusion phases.

3.4.   Separation of perfusion components

Once the above described post-treatment is applied on the raw data, we use the algorithm detailed in section ref, with the goal of separating small scale features (blood vessels) and the large scale (parenchym perfusion) features in the images sequences. The Figures 4 and 5 show the separation of these two components on the frontal and sagittal views before and after thrombolysis for respectively Patient 1 and Patient 2 at the time of maximal capillary perfusion (between arterial perfusion and venous drainage). The small scale component (middle column) is obtained by a back-projection (through $\Phi$) of the details coefficients $(\alpha_i^h, \alpha_i^v, \alpha_i^d)_{i = 1..6}$ on 6 levels obtained through the basis pursuit algorithm described in section 2. The large scale (right column) is simply the last approximation coefficient array $\alpha^a_6$. The fact that $\alpha$ is projected on the set of constraints $\{\Phi \alpha = Y\}$ ensures that, by summing of the small scale and the large scale image, we recover is the initial image (left column).

3.5.   Posterior medical analysis of the processed images

The Figure 4 illustrates the angiogram before and after successful mechanical thrombectomy for Patient 1. The even lines correspond to the pathological condition and the odd lines to the images after successful mechanical thrombectomy. The black arrow illustrates the lack of enhancement of the impaired ischemic area that is reperfused on the final set of post therapeutic images. The parenchymogram phase can be used as a confirmation of the normal findings, nevertheless the retrograde opacification of the external carotid artery has to be considered when interpretating the parenchymogram. It is in agreement with the reading methodology that includes the three phases analysis to overcome the artifacts that can be misinterpreted using only one of the three phases.

The Figure 5 illustrates another set of images from Patient 2. On the arterial and venous phases form the pretreatment set of images it is difficult to identify pathological ischemic area. Whereas on the parenchymogram phase, a lack of enhanced area is easily depicted (grey arrow). The findings are similar on the post treatment acquisition. The perfusion seems to be better on the post treatment parenchymogram whereas still confusing at the arterial and venous phases. Indeed, it could be misinterpreted as a superposition arising from the external carotid artery. On the light of the parenchymal results, the two others phases provide supplementary medical information that were not initially identified.

4.   Conclusion

In this article we have shown that basis pursuit in the context of undecimated wavelet transforms yield interesting applications for the analysis of angiographic sequences. In particular the discrimination of the arterial/venous and the capillary (or parenchymographic) phases with a scale criterion can be very efficient. Further improvements of this method (for instance by exploiting the time dynamics of each phase) should lead to better results and are currently under investigation. Our long term goal is to introduce physiologically-driven image representations that are both computationally tractable and that capture in the sparsest way possible the essential information for the medical experts.

By seeking the sparsest representation of angiograms we also set the scene for the use of statistical learning for the detection of subtle anomalies in a patient, which would lead to improve the existing qualitative scoring methods for the prediction of the patient state after thrombolysis.

In a near future we will consider applying this type of methodology to other organs, such as heart, kidney or liver angiographies. We expect each case to require a different pre-treatment in order to take into acocunt the specificity of each case. For instance, coronarographies may require non-rigid registration before any substraction is applied in order to compensate for the large mechanical deformations of the heard during the aquisition (this may even require the use of a biomechnical model of the heart).

Acknowledgements

The authors whis to thank the LabCom I3M involving Siemens for financial support. The authors with to thank Frederic Bosio (LMA, Poitiers) for the insight he provided us with on Laurent Polynomials. All the numerical calculations performed in thi article were conducted using the following python libraries : numpy/scipy, pydicom, pywavelets ^[5], scikit-image and matplotlib.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

In this appendix we give further details on the exact reconstruction condition and pseudo-inverse condition for FIR filters bank, restricting the study to one dimensional transforms on one level (for the sake of simplicity). We define the convolution of $y$ with a Finite Impulse Response (FIR) filter $h$ by :

we assume that $y$ is extended by periodicity (i.e. $y_{i+N} = y_i$). The discrete Fourier transform of $f$ writes:

hence :

If we define the $z$-transform of $(f_l)_{l = -a, ..a}$ as :

and if we denote $z_j = e^{-2 \pi i j \over N}$, we then have :

Let us consider a 1-Level UWT with kernels $(h, g, \tilde{h}, \tilde{g})$. The analysis writes :

where $\bar{g}$ (r.p. $\bar{h}$) is the reversed filter of $g$ (r.p. $h$), i.e. $\bar{g}_l = g_{-l}$. Note that the $z$-transform of a reversed filter is the conjugate of the $z$-transform of the said filter. The synthesis writes :

Perfect reconstruction happens when $\tilde{y} = y$, in other words, when :

which, in the Fourier space writes :

hence perfect reconstruction is achieved if and only if, for $j = 1..N-1$:

which, since $\bar{G}(z_j) = \overline{G(z_j)} = G(\bar{z_j}) = G(z_j^{-1})$ (since $|z_j| = 1$), leads to the perfect reconstruction condition:

Denoting $\Phi$ the synthesis operator and $\mathcal{W}$ the analysis operator, we clearly have that $\Phi \mathcal{W} = \text{Id}$, which means that $\mathcal{W}$ is a right-inverse of $\Phi$. Because there is no decimation step in this analysis/synthesis process, there is no aliasing condition to be satisfied here. Let us seek the required conditions that need to be satisfied so that $\mathcal{W}$ is the pseudo-inverse of $\Phi$. The only additional condition we need to have is that $\mathcal{W} \Phi$ is Hermitian. We have :

Then, in the Fourier space we have :

Hence, $\mathcal{W} \Phi$ is Hermitian if and only if the matrix (again we exploit the fact that ($\bar{G}(z) = G(z^{-1})$ and $\bar{F}(z) = F(z^{-1})$ for $|z| = 1$) :

is Hermitian, which leads to the following conditions:

We now restrict our study to symmetric filters, which satisfy the property : $\overline{F(z_i)} = F(\overline{z_i}) = F(z_i^{-1}) = F(z_i)$. In this case, the two first conditions (4.15) and (4.16) are satisfied automatically. Hence, from (4.11) and (4.17) for symmetric FIR, $\mathcal{W}$ is the pseudo-inverse of $\Phi$ if and only if, for $j = 0..N-1$ :

Let $m$ be the maximum width of the stencil of the filter bank, if we consider that $N > 2m$ (which is usually the case in practice: we chose stencils much smaller than the image size), the above expression can be recast as equations on Laurent polynomials

The following theorem states that it is not possible to find a symmetric FIR filter bank that satisfy both the exact reconstruction condition, and for which the analysis operator is the pseudo-inverse of the synthesis operator.

Theorem. Let $H$, $\tilde{H}$, $G$, $\tilde{G}$ be symmetric Laurent polynomials (i.e. $Q(X) = Q(X^{-1})$ for any such polynomial). If the equations (4.20)-(4.21) are satisfied, then $H$, $\tilde{H}$, $G$, $\tilde{G}$ are constant.

Proof. Since $H$, $\tilde{H}$, $G$, $\tilde{G}$ are symetric we can define $\mathcal{H}, \, \mathcal{G}, \, \tilde{\mathcal{H}}, \, \tilde{\mathcal{H}} \in {\mathbb{R}}[X]$ by : $H(X) = \mathcal{H}(X+X^{-1})$, $G(X) = \mathcal{G}(X+X^{-1})$, $\tilde{H}(X) = \tilde{\mathcal{H}}(X+X^{-1})$, $\tilde{G}(X) = \tilde{\mathcal{G}}(X+X^{-1})$. Then equations (4.20) and (4.21) imply:

which means that the polynomials $\mathcal{G}+i \mathcal{H}$, $\tilde{\mathcal{G}}+i \tilde{\mathcal{H}} \in {\mathbb{C}}[X]$ are inverses of each other and are therefore constant. Since $\mathcal{H}, \, \mathcal{G}, \, \tilde{\mathcal{H}}, \, \tilde{\mathcal{H}} \in {\mathbb{R}}[X]$, each of these polynomials is constant.