# Patent

## Title
Variational approach on whole body SPECT/CT registration and zipping

## Abstract

A method and a system for implementing the method for simultaneously registering and zipping a multiple scan whole body SPECT/CT image. The method includes the steps of simultaneously registering and zipping multiple input images and re-sampling the registered images. The step of simultaneously registering and zipping multiple input images is accomplished by initially aligning the images to be registered with each other, aligning the images with a reference image, and adjusting the alignment of the images with each other.

## Background
BACKGROUND DESCRIPTION
1. Technical Field
This invention relates generally to SPECT and CT imaging. Specifically, it relates to registering SPECT images and CT images of the same patient regions, and “zipping” together SPECT images of different portions of a whole body scan to provide a single whole body image.
2. Background of the Invention
When taking whole body Single Photon Emission Computed Tomography (“SPECT”) and Computed Topography (“CT”) scans, in many machines the detector's field of view (“FOV”) is limited. It is therefore often necessary to take several separate scans for SPECT, which overlap in the z direction (see FIG. 1), at two or more different positions with respect to a patient. These separate scans must then be “zipped” or appended together after reconstruction.
When such zipping takes place, it is usually impossible to determine the proper zipping position in the overlapping region where zipping occurs, in the z direction (see FIG. 1). When two adjacent images are so overlapped, the proper dividing line could be anywhere in the overlapping region. Presently there are no satisfactory methods to determine this position and to zip the two images together based either on relative bed positions or on image positions.
Even in cases where the zipping position can be approximated, i.e. when the full reconstruction range is used in the overlapping region, other factors may hinder a satisfactory zipped whole body image. These factors may include: bed deflection, patient motion, and image edge handling in 3D reconstruction algorithms with CT attenuation correction.
In many current methods, image registration (i.e., between the SPECT and CT images) and zipping (i.e., of two SPECT images of overlapping adjacent patient regions) are done completely separately. Auto-zipping is done after the multiple whole body SPECT images have been individually registered with the CT image. The separate registration and zipping processes do not generate satisfactory whole body images.
To solve these problems, it is desired to merge the registration task and the zipping task into a single optimization task.
SUMMARY OF THE INVENTION
Therefore, according to the present invention a method for simultaneously registering and zipping a multiple scan whole body SPECT/CT image is provided. The method includes the steps of (a) simultaneously registering and zipping multiple input images and (b) re-sampling the registered images. The step of simultaneously registering and zipping multiple input images is accomplished by (i) initially aligning the images to be registered with each other, (ii) aligning the images with a reference CT image, and (iii) adjusting the alignment of the images with each other.
In order to determine the best registration, which is used to generate a registered output, the method uses the equation:

Mtotal(φ):=ΣjM(U,Vj∘φj)=max,

which is subject to

Vtotal(φ):=Σi≠jν(Vi∘φi,Vj∘φj)=min

where U is a CT reference image, V is a set of SPECT images, and φ is a transform.

Further provided is a system for implementing the method that includes a SPECT/CT scanning device, a processor that receives scans from the SPECT/CT scanning device, and software that inputs multiple images, simultaneously registers and zips the images and outputs a single, unified, registered image.



BRIEF DESCRIPTION OF THE DRAWINGS
The invention will now be described in greater detail in the following by way of example only and with reference to the attached drawings, in which:
FIG. 1 is an example of two SPECT scans with an overlapping region.
FIG. 2 is an example of unsatisfactory whole body zipping for two SPECT scans.
FIG. 3 is a diagram of the method for simultaneously registering and zipping a multiple scan whole body SPECT/CT image
FIG. 4 is a diagram of an algorithm for executing the method.
FIG. 5 is a system for simultaneously registering and zipping a multiple scan whole body SPECT/CT image.



DETAILED DESCRIPTION OF THE INVENTION
As required, disclosures herein provide detailed embodiments of the present invention; however, the disclosed embodiments are merely exemplary of the invention that may be embodied in various and alternative forms. Therefore, there is no intent that specific structural and functional details should be limiting, but rather the intention is that they provide a basis for the claims and as a representative basis for teaching one skilled in the art to variously employ the present invention.
FIG. 1 depicts two whole body SPECT scans 110 and 120 of overlapping adjacent regions of a patient, where the solid lines represent the image reconstruction range for each image. It is desirable to zip these two SPECT scans 110 and 120 together to create one image. The two SPECT scans 110 and 120 will be zipped together somewhere in the overlapping region 140. However, the proper dividing line 130 could be anywhere in the overlapping region 140. If the wrong dividing line 130 is chosen, an unsatisfactory final image will be produced (see FIG. 2).
FIG. 3 diagrams the method 300 according to one embodiment of the present invention, which may find the proper dividing line 130 from which to zip the SPECT images together to form a single whole body image. The two main processes in the method 300 are (1) simultaneous registration and zipping (340) and (2) re-sampling (350). The simultaneous registration and zipping 340 may include imputing a reference CT image 310 and imputing whole body SPECT scans 1 (320) through K (330). The simultaneous registration and zipping 340 may not only determine the best alignment (registration) between each individual SPECT image 320 through 330 and the reference CT image 310, but also may determine the best alignment among the K SPECT images 320 through 330 themselves. After the simultaneous registration and zipping 340 is complete, the registered multiple images 360 may be re-sampled in process 350. The re-sampling 350 may sample the multiple registered images 360 (which may have overlaps) to generate a single unified output 370.
FIG. 4 depicts an algorithm 400 for completing the method 300. The algorithm may consist of the steps of (1) scanning the SPECT images (“V1, . . . , VK”) 410, (2) finding an initial alignment among V1, . . . , VK images 420, (3) imputing a reference CT image (“U”) 480, (4) maximizing the total distance between U 480 and each SPECT image (“Vi”) 430, (5) minimizing the overall variation among overlaps of V1, . . . , VK 450 (6) repeating (440) the maximizing 430 and minimizing 450 steps until the best registration is found, (6) re-sampling 460, and (7) registering an output 470.
Finding the initial alignment of the K images V1, . . . , VK 420 may be done based on either bed positions or image positions of the K images. The maximizing 430 and minimizing 450 steps may be formulated as an optimization problem as follows.
Let R3 denote the usual three-dimensional Euclidian space. An image may be defined as a function from R3 to R which satisfies certain regularity conditions. Given two images U and V, where U is a reference image and V is the image to be registered towards U, the objective of the registration between these two images is to find a proper transformation

φ:R3→R3 

such that U and V∘φ are best matched in accordance with a certain objective measure, where V∘φ denotes the registered version of V with V∘φ(x)=V(φ(x)) for xεR3.

In a multiple input registration and zipping setting, there may be one reference image U (the CT image) and a set of K SPECT images {Vj} to be registered (the multiple whole body SPECT images). It may be necessary to find K best transformations {φj} under certain optimization criteria, where each φj represents the best registration between U and Vj. The set of functions {φj} cannot be found separately because their domains have overlaps in general, and those are the regions where transformations need to be adjusted to make the best zipping for the neighboring two images.
If φ=(φ1, . . . , φK), the maximizing (430) and minimizing (450) steps may be formulated as an optimization problem as follows:
Given one reference image U:R3→R3 and a set of K images {Vj} to be registered, where Vj: R3→R3, j=1, . . . , K, find a transformation φ such that

Mtotal(φ):=ΣjM(U,Vj∘φj)=max

subject to

Vtotal(φ):=Σi≠jν(Vi∘φi,Vj∘φj)=min  (1)

where

M(U,Vj∘φj) measures the similarity between the reference image U and the transformed image Vj∘φj; Mtotal is the sum of all M(U,Vj∘φj); ν(Vi∘φi,Vj∘φj) measures the variation between the two registered images Vi∘φi and Vj∘φj at their overlapped region; and Vtotal is the sum of all ν(Vi∘φi,Vj∘φj).
One implementation for the optimization of problem (1) is to set the objective functional as

J(φ):=−Mtotal(φ)+λVtotal(φ)  (2)

and search for φ* such that J(φ*)=min, where λ is a constant to be determined. A gradient based steepest descent method may be used to seek the minimum of the functional. First the gradient ∇J(φ) may be calculated, and then updates in the search for the optimal transformation φ may be made according to

φn+1=φn−μ∇J(φn),μ>0,n=1,2,3, . . .

where μ is a constant used to control the convergence rate.

Thus, this registration algorithm searches for the best registration φ between the set of images {Vj} and the reference image U in such a way that the individual images Vj(1≦j≦K) are optimally aligned with respect to the reference image U (in the sense of Mtotal=min).
Once the best registration φ has been found, it may be used in the final re-sampling operation to generate a registered output. Note that in the conventional image registration setting where the re-sample is based on one transformation function φ only, the multiple input re-sample algorithm in this operation must handle the multiple transformation functions {φj}. In particular, interpolation is needed in the overlapped domain of the functions.
Often some type of regularization is needed because the image registration problem is ill-posed.
Let the transformation function φ:R3→R3 be the deformation map defined by

φ(x)=x+u(x)

where u is a proper function from R3 to R3.

For the similarity measure M between two images U and V, one may use the popular mutual information defined by




  
    M
    ⁡
    
      (
      
        U
        ,
        V
      
      )
    
  
  =
  
    
      ∫
      
        
          R
          3
        
        ×
        
          R
          3
        
      
      
        
      
    
    ⁢
    
      
        
          p
          
            U
            ,
            V
          
        
        ⁡
        
          (
          
            u
            ,
            v
          
          )
        
      
      ⁢
      log
      ⁢
      
        
          
            
              p
              U
            
            ⁡
            
              (
              u
              )
            
          
          ⁢
          
            
              p
              V
            
            ⁡
            
              (
              v
              )
            
          
        
        
          
            p
            
              U
              ,
              V
            
          
          ⁡
          
            (
            
              u
              ,
              v
            
            )
          
        
      
      ⁢
      
        
      
      ⁢
      
        d
        
          (
          
            u
            ,
            v
          
          )
        
      
    
  




where pU and pV are the probability densities of the pixel values of the images U and V, respectively; pU,V is the joint probability density of the pixel values of images U and V.

For the variation measure v between two overlapped images F and G, one can use the sum of the squared difference defined by

v(F,G)=∫O(F(x)−G(x))2dx 

where O denotes the overlapped region between the two images.

Under these notation, equation (1) can be formulated as the following variational problem:

minimize J(u1, . . . , uk)+rS(u1, . . . , uk)  (3)

where J is defined as in equation (2), S is a regularization term, and r>0 is a regularization parameter. In many cases, the regularization term S can be defined as a bi-linear form of B:

S(u1, . . . , uk)=Σj=1K∫D<B(uj),B(uj)>dx,D⊂R3 

where B is a differential operator, and <,> denotes the inner product in L2(R3). L2(R3) is the completion of the continuous functions with respect to the L2-norm. For example, for elastic registration, the integral term in the above expression can be represented as





  
    
      ∫
      D
      
        
      
    
    ⁢
    
      
        〈
        
          
            B
            ⁡
            
              (
              
                u
                j
              
              )
            
          
          ,
          
            B
            ⁡
            
              (
              
                u
                j
              
              )
            
          
        
        〉
      
      ⁢
      
        
      
      ⁢
      
        d
        x
      
    
  
  =
  
    
      ∫
      D
      
        
      
    
    ⁢
    
      
        {
        
          
            
              α
              4
            
            ⁢
            
              
                ∑
                
                  i
                  ,
                  
                    j
                    =
                    1
                  
                
                3
              
              ⁢
              
                
                  (
                  
                    
                      
                        ∂
                        
                          x
                          i
                        
                      
                      ⁢
                      
                        u
                        j
                      
                    
                    +
                    
                      
                        ∂
                        
                          x
                          j
                        
                      
                      ⁢
                      
                        u
                        i
                      
                    
                  
                  )
                
                2
              
            
          
          +
          
            
              β
              2
            
            ⁢
            
              
                (
                
                  ∇
                  
                    ·
                    u
                  
                
                )
              
              2
            
          
        
        }
      
      ⁢
      
        
      
      ⁢
      
        d
        x
      
    
  




where α and β are the so-called Lamé constants, and ∇● is the divergence operator. Note that u is the function used to define the non-rigid transformation function φ(x)=x+u(x). In a multiple input registration setting, with K input images to be registered, there should be K such u's.

Using a proper discretization technique, the regularized minimization problem can be implemented as an iterative algorithm.
FIG. 5 shows an example system 500 that uses the method 300. The system may be comprised of a SPECT/CT scanning device 510 and a processing device. The processing device may obtain scanned images from the scanning device 510 and may run software that implements the algorithm 400 to output a single registered and zipped whole body image. The system may include a monitor 520 for displaying data, operating instructions, etc. from the processing device.

## Claims

1. A method for simultaneously registering and zipping a multiple scan non-whole body SPECT image with a whole body CT image of a same region, comprising:
(a) initially aligning multiple non-whole body input SPECT images with each other;
(b) separately aligning each of said multiple non-whole body input SPECT images with said whole body CT image;
(c) adjusting alignment between said multiple non-whole body input SPECT images based on alignment with said whole body CT image; and
(d) re-sampling the registered images to obtain a single output zipped image.

2. The method of claim 1, wherein step (a) includes:
(i) initially aligning the images to be registered with each other;
(ii) aligning the images with a reference image; and
(iii) adjusting the alignment of the images with each other.

3. The method of claim 2, wherein step (i) is based on bed positions.

4. The method of claim 2, wherein step (i) is based on image positions.

5. The method of claim 2, wherein the reference image is a CT.

6. The method of claim 2, wherein steps (ii) and (iii) are repeated until a final alignment has been attained.

7. The method of claim 6, wherein steps (ii) and (iii) are accomplished by finding a transformation Φ wherein:


Mtotal(φ):=ΣjM(U,Vj∘φ=max

subject to


Vtotal(φ) :=Σi≠jν(Vi∘φi,Vj∘φj)=min


wherein:
M (U,Vj∘φj) measures the similarity between the reference image U and the transformed image Vj∘φj;
Mtotal is the sum of all M(U,Vj∘φj);
ν(Vi∘φi,Vj∘φj) measures the variation between the two registered images Vi∘φi, and Vj∘φj at their overlapped region; and
Vtotal is the sum of all ν(Vi∘φi, Vj∘φj).

8. The method of claim 7, wherein the transformation φ is found by searching for a φ* such that J(φ*)=min using the equation:
J(φ):=−Mtotal(φ)+λVtotal(φ)
wherein:
λ is a constant to be determined.

9. The method of claim 8, wherein φ* is found by calculating the gradient ∇J(φ) and updating the search for the optimal transformation φ according to:


φn+1=φn−μ∇J(φn),μ>0,n=1,2,3, . . .

wherein:
μ is a constant used to control the convergence rate.

10. The method of claim 7, wherein the transformation φ(for φ(x)=x+u(x)) is found by minimizing:


J(u1, . . . , uk)+rS(u1, . . . , uk)

wherein:


J(φ):=−Mtotal(φ)+λVtotal(φ)


wherein, λ is a constant to be determined;
S is a regularization term; and
r>0 is a regularization parameter.

11. The method of claim 10, wherein:


S(u1, . . . , uk)=Σj=1K∫D<B(uj),B(uj)>dx,D⊂R3 

wherein:
B is a differential operator; and
<,> denotes the inner product in L2(R3);
wherein L2(R3) is the completion of the continuous functions with respect to the L2-norm.

12. The method of claim 1, further comprising:
(c) generating a registered output.

13. A system for simultaneously registering and zipping multiple scan whole body SPECT/CT images, comprising:
a SPECT/CT scanning device;
a processor in communication with the SPECT/CT scanning device; and
software executing on the processor, wherein the software
initially aligns multiple non-whole body input SPECT images with each other;
separately aligns each of said multiple non-whole body input SPECT images with the whole body CT image
adjusting alignment between the multiple non-whole body input SPECT images based on alignment with the whole body CT image; and
re-samples the registered images to obtain a single output zipped image.

14. The system of claim 13, wherein the software:
initially aligns the images to be registered with each other;
aligns the images with a reference image; and
adjusts the alignment of the images with each other.

15. The system of claim 14, wherein the software finds a transformation φ wherein:


Mtotal(φ):=ΣjM(U,Vj∘φj)=max

subject to


Vtotal(φ):=Σi≠jμ(Vi∘φi,Vj∘φj)=min


wherein:
M(U,Vj∘φj) measures the similarity between the reference image U and the transformed image Vj∘j;
Mtotal is the sum of all M(U,Vj∘φj);
ν(Vi∘φi, Vj∘φj) measures the variation between the two registered images Vi∘φi and Vj∘φj at their overlapped region; and
Vtotal is the sum of all ν(Vi∘φi, Vj∘φj).

16. The system of claim 15, wherein the software finds the transformation φ by searching for a φ* such that J(φ*)=min using the equation:


J(φ):=Mtotal(φ)+λVtotal(φ)

wherein:
λ is a constant to be determined.

17. The method of claim 16, wherein the software finds φ* by calculating the gradient ∇J(φ) and updating the search for the optimal transformation φ according to:


φn+1=φn−μ∇J(φn),μ>0,n=1,2,3, . . .

wherein:
μ is a constant used to control the convergence rate.

18. The system of claim 15, wherein the software finds transformation φ (for φ(x)=x+u(x)) by minimizing:


J(u1, . . . , uk)+rS(u1, . . . , uk)

wherein:


J(φ):=−Mtotal(φ)+λVtotal(φ);


wherein, λ is a constant to be determined;
S is a regularization term; and
r>0 is a regularization parameter.

19. The system of claim 18, wherein:


S(u1, . . . , uk)=Σj=1K∫D<B(uj,B(uj)>dx,D⊂R3 

wherein:
B is a differential operator; and
<,> denotes the inner product in L2(R3);
wherein L2(R3) is the completion of the continuous functions with respect to the L2-norm.