• No results found

Image Reconstruction for Pet Using Fuzzy Potential

N/A
N/A
Protected

Academic year: 2023

Share "Image Reconstruction for Pet Using Fuzzy Potential"

Copied!
10
0
0

Loading.... (view fulltext now)

Full text

(1)

2004 IEEE Workshop on Machine Learning for Signal Processing

IMAGE RECONSTRUCTION FOR P E T USING FUZZY POTENTIAL

Partha Pratim Mondal and K. Rajan

Department of Physics, Indian Institute of Science, Bangalore 560012, India Phone: +91 080 22933280

E-mail: partha(rajan) Qphysics.iisc.ernet.in

Abstract. Maximum likelihood (ML) and Maximum a-posteriori (MAP) are the most widely used algorithms for emission tomogra- phy (ET). MAP-approach heavily depends on Gibbs hyper-parameter for noise free reconstruction. Choosing a correct hyperparameter is dimcult and time consuming task. Recently proposed median root prior (MRP) algorithm is a good alternative, but are prone t o step like streaking effect. In this research work, a fuzzy logic based approach is proposed to overcome these shortcomings. Unlike tra- ditional potential function, a fuzzy potential is used for modeling inter pixel interaction. Two basic operations viz. edge detection and fuzzy smoothing are performed sequentially during each iter- ation. The first operation is employed for the detection of edges (if present) in all the eight directions of a 3 x 3 neighborhood win- dow. The second operation uses this edge information t o perform fuzzy smoothing. Due t o recursive nature of the reconstruction al- gorithm, both these operations are employed iteratively to reduce heavy noise produced due to dimensional instability

[lo].

Simu- lated experimental results are obtained to show the feasibility of the proposed approach. These algorithms are also compared with other approaches such as, MAP and M R P by numerical measures and visual inspection. Algorithm evaluation shows promising re- sults.

INTRODUCTION

In most of the medical imaging techniques, the data aquisition system ends up acquiring incomplete data. Incompleteness is due to the geometry of the imaging system and the nature of physical process involved. Penalized algorithms are found to be most efficient and reliable when dealing with incomplete data problem. A descriptive understanding of iterative image re- construction algorithms for emission tomography (ET) can be found in [1][2].

Ilerative algorithnis offer good reconstruction at the cost of large computa-

(2)

tional complexity can be addressed to some extent by parallel computation [3] and by using fast processors (like, digital signal processor (DSP), me- dia processor(MP) etc.) [4]. Many of the present day diagnostic imaging modalities like positron emission tomography (PET), single photon emission computed tomography (SPECT) demand good quality images. Though iter- ative algorithms like, maximum likelihood (ML) [1][2], maximum a-posteriori (MAP) [5][6][7] and median root prior (MW) [8][9] algorithms are capable of generating good quality images but at the cost of artifacts like noise, over- smoothness, streaking effect

[SI

and hence require improvement.

Both ML and MAP are iterative data fitting problem, in which t h a t esti- mate is chosen as the solution which makes the measurement data most likely.

Regularization involves controlling the solutions (reconstructed images) that are desired. Regularization can be interpreted as a feedback process which demands continuous passage from image space t o data space and vice-versa until pseudo-projections and measured projections match to the desired level of accuracy. The prior knowledge helps in producing the desired effect in the reconstrncted image. For example, smoothing priors produce smooth reconstruction [5][6][7] whereas edge-preserving priors [8][9] produce sharp reconstruction.

Fuzzy techniques have been successfully applied in image processing ap- plications such as image restoration [ll], morphology

[E]

etc.. In the present mrork we have extended fuzzy concepts t o image reconstruction in PET. Prior distribution is defined by Gibbs distribution and the potential (which defines the nature of nearest neighbor interactions) is modeled using fuzzy rules.

IMAGE RECONSTRUCTION ALGORITHMS FOR PET The measurements in PET

, aj ,

j=l:

...,

M are modeled as independent Pois- son random variables i.e, yj

-

P d s s a ( C ~ , X i p i j ) for j=l,

...,

M where X i

,

i=1,

..., N

are the mean parameters of the emission process and pjj is the probability that an annihilation in the i t h pixel is detected in jfh detector.

The likelihood function i.e, the conditional probability for observing

Y

= y given the emission parameter A = X is the joint probability of the individual Poisson process i.e,

Reconstruction algorithm proceeds by finding that estimate which maxi- mizes the objective function i.e, ML or MAP depending upon Lhe estimation technique employed. In ML, the objective function is taken as the likelihood function or equivalently log-likelihood function. Then, the task is t h e de- termination of that estimate A“‘ which maximizes the objective function

(3)

i.e,

y$

lOSP(Y/X)

1

(2)

A"L =

The iterative equation solution for the above ill-posed reconstruction problem is first. given by Shepp and Logan for P E T [l]. The iterative update for generating (IC

+

l ) t h estimate from kth estimate is,

On the other hand, MAP algorithm determines t h a t estimate

X M L

as the solution tihich maximizes the post,erior densit.y function P(A/y) or equiva- lently the log of P(A/y). Given a suitable prior P(X), MAP-reconstruction can be formulated

as,

AA'Ap =

y2y[

logP(y/X)

+

logP(X) ] (4) Image field is assumed as Markov random Geld (MRP) [6] and by Hammerseley- Clifford theorem 1131, image

X

is characterized by Gibbs distribution,

( 5 ) 1

-BC,CjrNi

W i j V ( A i : A j )

P(A) =

F e

where, Z is t.he normalizing constant. for the dist,ribut,ion,

fi

is the Gibbs hyper-parameter, wij is the weight of pixel j e N ; 151, .Vt is the nearest neighbor set of pixcl i and V(Xc>Aj) is termed as t,hc potential at site i duc t o the nearest neighbor elements j .

Solution for eqn.(l) is very difficult due to the complicated nature of prior.

Green [5] has proposed one step late (OSL) approximation for an iterative update to the h4.4P-problem and is given by,

Given thc iterative OSL-algorithm (cqn.(S)), the next step is the proper modeling of the interacting potential V(X;,X,) between the pixel at site i and its neighbors jeA'i. 4 large number of potentials (see [6][5][7]) have been suggested in the literature to produce desired image characteristics. It is a geueral nature of these potentials function to sniooth the edges irrespective 01 the density class. The M4P reconstruction also involves prior hiowledge of the hyperparameter values of the prior distribution, which are rarely avail- able.

Recently, Alenius et. al. [S][9] has suggested median root prior algorithm for P E T image reconst.ruction. The MRP formulation requires comparing of

(4)

of

More details can be found in [8][9]. The MW-algorithm is given by:

The penalty term in (7) is piecewise linear with respect to the difference between the center pixel and local median. This difference can be non-zcro only if the image is not locally monotonic. Thus, noise is PO nalized because it is non-monotonic inside a small neighborhood. The local monotonic part represents the signal [9].

FUZZY LOGIC BASED POTENTIAL FUNCTION

Fuzzy rules provide a robust estimation for edge detection. A fuzzy logic based potential is proposed for edge-preserving reconstruction. This c0nsist.s of two basic operations : fuzzy filtering and fuzzy smoothing. The proposed fuzzy technique is similar t o that used by Ville et. al. for image restoration (111. In the present work, the idea is expanded and adapted for PET.

Fuzzy Derivatives for Edge Detection

The deribative V k ( i , j ) for pixel at (i, j ) along the direction A at kth iteration is defined as,

V"(i,j)

A =

( A"(*,*)

-

P ( i , j )

)A

(8) where, A&(*,

*)A

represents the nearest pixel value along the unit directional vector fi. For example, the derivative along the west direction is defined as,

Vk(i,j) f i 7 = ( X"i,j

-

1)

-

X k ( i , j ) )I@

For identifying edge in a particular direction, three elementary derivatives are chosen (see fig.1). For example to detect an edge along Wpirection, the derivatives are : V k ( i ; j ) W , V'(i - 1 , j ) W and V k ( i

+

1 , j ) W for a 3 x 3 window. It is safe to assume that if 2 out of 3 elementary deritatives are small: the edge is absent in the neighborhood. This is termed a 2 3 rule. For A = W the fiizzy derivative is defined as follows :

V k ( ( i , j ) and O L ( i

-

1 ; j ) are small or V * ( i , j ) and V k ( i

+

1 ; j ) are small or

{

V k ( i - 1,j) and V k ( i

+

1,j) a r e small

If

Then, V;(i,j)lii is small

k . . ^

Else,VF(a,3)W is large (9)

(5)

64 (b)

Figure 1: (a) 3 5 3 neighborhood of a central pixel ( i j ) , showing the directional derivative along W , (b) All the directions along which the the derivatives are taken.

Similarly, the values of the fuzzy derivatives V$(i,j)A for all the directions i.e,

@,

J??, s,h%i', S f V , SE, AfE} are calculated.

Membership Function

Fuzzy sets are best represented by membership function. Membership func- tion gives the degree of belongingness within the set. The membership func- tion /L of a fuzzy set X, maps the dements of S t o the interval [0,1] i.e, A'

4 [O;

11. To compute the value that expresses the degree t o which the fuzzy derivative is small, we make use of fumy set small. Membership func- tion for the property small is defined as,

k . .

-

[small, if V ( i : j )

5 o

'

( t ' 3 ) n = large, otherwise.

where, k represents t,he iteration number and A represents the direction.

In

iterative image reconstruct,ion procedure, t,he estimate improves aft.er each iteration and hcnce requires an updated mcmbcrship function. Duc t o the iterative nature of the MAP algorithm, proposed menibership function automatically updates the set small after each iteration. Similarly mem- bership function is definedAfor all fuzzy 4erivatives V ~ ( i : j ) f i along all the directions viz. {@,

fi,

S, AW, S W ,

SE,

N E } .

Fumy P e n a l i z a t i o n

The next step is the penalization of pixels for ahich edges axe not detected in the considered nearest neighborhood window, keeping others unaltered. The following rule is used for penalization :

If V $ ( i : j ) A is small, Then a k ( i , j ) A = V k ( i , j ) A .

Else, A k ( i , j ) A = 0 (11)

where, Ak(i,j)6 is the feedback a t site ( i : j ) due to the adjacent pixel in the direction A. Contributions are taken from all the eight directions. Hence, the total correction term AF(i,j) for pixel a t ( i , j ) considering all the directions after kt" iteration is given by,

(6)

A $(i ),

x<=xt

the OSL-algorithm modifies to,

where, coordinates (i,j) is denoted by a single coordinate {i' = (i - 1)

* v%+j]. In

the iterative image reconstruction procedure, the final correction term is fed back t o update the pixel after each keration. The iterations are continued until acceptable convergence is obtained.

SIMULATED E X P E R I M E N T A L R E S U L T S Simulated PET System

The algorithm was tested on a simulated P E T system. The P E T system con- sists of a ring detector with 64 detectors and the object space is decomposed into 64 x64 square pixels. The object space is a square region inscribed within the circle of detectors. Each element of the probability matrix p i j defines the probability of a photon getting detected in the detector j after emanating from t h e object pixel i . For simplicity, we assumed that p;j depends only on the gcomctry of thc measuremcnt system. This is takcn as the angle B i j seen by the center of the pixel i into the detector tube j 111 i.e, p i j =

%.

Before

the reconstruction begins, the probability matrix P

=

[ P i j ] ,i = 1

,

...~ N and j = 1 , ..., M is compnted and stored.

For simulating measurement data, a Monte Carlo procedure is used in which each emission is simulated as follow [1][15]. First, a random pixel is chosen in the test phantom. Thc conccntration of the radionuclei at the given pixel is assumed t.o he proportional t o the emission density of the pixel. For each of the accepted emission point, a randomly oriented line (between 0 to ?i

radians) is selected and the pair of detectors this line intersects is found. The random line corresponds to the direction in which the pair of annihilated phobons travel and the pixel corresponds to the annihilation point in the phantom. The detectors with which this line irkersects are assumed to detect this annihilation event and the count in the corresponding detector tube is inclemented. In this way all the emissions are simulated and counted in the respective tubes. This is used as the measurement data for reconstruction.

The mathematical phantom (also known as Shepp-Logan phantom) used in the prcsent study is made up of nine elliptic objccts having different sizes, orientations and density values. In fig.3, the image of the phantom is shown.

We have used a source image with 100,000 counts.

(7)

Figure 2: (a) Log-likelihood vs. iteration plot of the proposed algorithm; (b) Resid- ual error vs. iteration plot of the proposed algorithm.

A l g o r i t h m E v a l u a t i o n

All the evaluation test,s defined in this sect.ion are carried out on a simulated P E T system. The proposed algorithm is used on a 3 x 3 neighborhood win- don,. The rcsults arc compared with MAP and MFLP algorithms. MAP with potential V(Xi - X j ) =

CjrNi(Xi

= Xj)? with

fl

= 2.5 x lo4 is used in this study. That rhoice of

fl

is considered which gives the best estimate. The per- formances of the proposed new algorithm are evaluated using two different inragchased quarititative criteria as given below :

Log-likelihood Test. The algorithms described in section

II

and III com- pute t,he estimate of the mission densities iteratively, hence log-likelihood function is an appropriate qualihtive measure. For an estimate A', the log- likelihood function l(X') at IC"' iteration is defined as,

M

l ( P )

=

E[-?+$ +

g j log$;

-

log(y,!)] (14)

j=1

where,

~;

= N

The log-likelihood d u e s of the reconstructed images obtained using MAP and proposed algorithm

are

plotted against the iterations in fig.Z(a). It is clearly evident t h a t log-likelihood for the proposed algorithm converges faster compared to MAP-algorithm.

R e s i d u a l Error. In iterative image reconstruction techniques, residual error is the most prefered evaluation test. This measures the deviation of the generated pseudc+projections $$ of the reconstructed image from the observed projection d a t a yj. Residual error

[(Ak)

at kth-iteration is given by:

A:Pjj is the pseudo-projection in the tube j .

where, @$ =

E,"=,

is the pseudo-projection in the tube j. In fig.Z(b), t h e residual errors of the reconstructed images using t.he proposed

(8)

( c ) proposed algorithm. Second row shows reconstruction afta 100 iterations for (d) MAP, (e) MFW, (f) proposed algorithm. For comparison, original test phantom is also shown.

algorithm is shown. This plots reveals the fact that residual error decreases with increase in iterat,ion.

Visual Inspection. Final and the most important test is the visual in- spection of the reconstructed images. In fig.3, (a,b,c) and (d>e,f) show the reconstructed images using MAP,

MIW,

proposed algorithm after 50 and 100 iterations respectively. For quality assessment, the original test image (phantom) is also shown. The images reconstructed using the proposed algo- rithm (fig.3(c),(f)) are more appealing and artifact free as compared t o those reconstruct.ed using MAP and MRP algorithms.

A small edge region of the reconstructed images using MAP: MRP and the proposed algorithm is shown in fig.4 for further insight. MAP-reconstructed image shows over-smoothness (see fig.4(b)), while step like streaking effects are evident in MRP reconstructed image (see fig.l(c)) [9]. It can be seen that the images generated using the proposed algorithm are more appealing and free from streaking artifacts which is prominent in MRP-reconstructed images.

Line Plot. L.ine plots (see fig.5) show- the image intensities along a prede-

(a)

6) i;j (ij

Figure 4: Zoomed version of the reconstruction after 100 iterations are show: (a) Original, (b) MAP, ( c ) MRF’, (d) Proposed Algorithm.

(9)

m

Algorithm

1 s p-->

Iru Y

-

Figure 5: Line plots for the reconstructed images are shown after 80 iterations: (n) Original, (b) MAP, (c) MRP, (d) Proposed Algorithm.

fined line in t,he respective reconstructed images. The intensities are plotted along the central horizontal line through the reconstructed images. Line plots are shown for MAP,

MRP

and original test phantom for comparison. It is evident from these plots t h a t MAP over smooths the reconstruction while MRP is capable of preserving the small details a t the cost of streaking effect (see fig.S(c)) 191. But the proposed algorithms generate reconstruction that are not too smooth and st.ill preserves the edges.

CONCLUSIONS

A new approach is presented towards a better edge preserving reconstruc- tion in PET. This is based on the application of fuzzy rule based techniques t o model the pixel-pixel interacting potential in image lattice. Two basic fuzzy operatiom are performed fuzzy rdlering and luzzy smoothing. These operations arc performed iteratively to reduce h e a y noise. Iterations are stopped when acceptable Convergence is obtained. Simulated experimental P E T studies reveal t h a t the proposed algorithm is capable of producing es- timates that are slightly better than the existing algorithms like MAP and MFlP. There is no fear of convergence as the log-likelihood function is shown to converge with increasing iteration. Decrease in residual error confirms that the proposcd algorithm is capablc of producing estimates for which thc gcnw- ated pseudoprojections match closely to the measured projections and hence is the most likely estimate t h a t has given rise t o the observed projections.

Visual representation of the recoustructed images using proposed algorithm is more appealing compared to M.4P and MRP reconstructed images. The results are very encouraging.

(10)

[l] Y. Vardi, L. A. Shepp, and L. Kaufmann, “A statktical model for positron emission tomography”, J. Amer. Stat. h o c . , vol. 80 ~ pp. 8-37, 1985.

[2] L.A. Shepp and Y. Vardi., “Maximum likelihood estimation for emission to- mography”, IEEE. Trans. on Med. Img. MI-1: pp. 113-121, 1982.

[3] C. 11. Chen and S. Y. Lee, ”Pardlelization of the EM algorithm for 3-D PET image reconstruction”, IEEE Ttans. h4ed. Img., Vol.10, No.4, pp.513-5Z2, Dec.

1991.

(41 K. Rajau; L. M. Patnaik and J. Fkmakrishna, Wigh speed computation of the EM algorithm for PET image reconstruction”, IEEE Ttans. Nucl. Sci., [5] P.J. Green, ‘Bayesian reconstruction from emission tomography data using

a modified EM algorithm”, IEEE ’hans. on Med. Img., vol.9, No.1, March, 1990.

[SI 2. Zhou, R. 14. Leaby and J. Qi, “Approximate maximum likelihood hyperpa- rameter estimation for Gibbs prior”, IEEE Trans. on Img. proc., Vo1.6, No.6, pp.841-861, June, 1997.

[7] J. Nuyts, D. Bequ, P. Dupont, and L. Mortehans, “A Concave Prior Penaliz- ing Relative Differences for Maximum-a-Posteriori R.econstruction in Emission Tomography”, IEEE Trans. on Nuc. Sci., Vol. 49, No. 1, pp.56-60, Feb. 2002.

[8] S. Alenius and U. Ruotsalainen, “Using Local Median as the Location of Prior Distribution in Iterative Emission Tomography Reconstruction”, IEEE Ttan.

Nucl. Sci., Vol.45, No.6, Dec. 1998.

[9] S. Alenius and U. Ruotsalainen, “Generalization of Median Root Prior Recon- struction”, IEEE ’hans. bled. Img., Vo1.21, No.11, Nov., 2002.

(101 E. Veclerov and .J. Llacer

,

“Stopping nile for MLE algorithm baed on sta- tistical hypothesis testing”, IEEE ’hans. an Med. Img., MI-6, pp. 313-319, 1987.

1111 D. V. Villc, M. Nachtegael, D. V. Wckcn, E. E. Kerre, W. Philips and I.

Lemahieu, “Noise Reduction by Fuzzy Image Filteringa ,IEEE Trans. Fuzzy Sys., Vol. 11, NO. 4, Aug. 2003.

(121 11. Nachtegael and E. E. Kerre, “Decomposing and constructing fuzzy mor- phological operations mer a-cuts: Continuous and discrete case”, IEEE Trans.

Fuzzy Syst., vol. 8, pp. 615626, Oct. 2000.

[13] J. Besag

,

“Spatial interaction and tlre statistical analysis of lattice systems”, JI. of Royal Stat. Soc. B, Vo1.36, pp. 192-236, 1974.

(141 S. Geman and D.Geman, “Stocastic relaxation, Gibbs distribution and the Bayesian restoration of images ”, IEEE Trans. Pat. Anal. Mac. Intell., vol.

PAMI-6, pp. 721-741, Nov., 1984.

[I51 h‘. Rajeevan, K. Rajgopal and G. Krishna. “Vector-Extrapolat,ed fast maui- nium likelihood estimation algorithms for emission tomography*, IEEE Trans.

on Med. Img., Vol. 11, h-0.1, March 1992.

(161 L. Kaufmann, “Implementing and accelerating the EM-algorithm for positron emission tomography”, E E E Trans. Med. Img., Vol. MI-6, pp.37-51, Mar.

1987.

VOI. 41, No.5, pp.o-5, 1994.

References

Related documents

Many MPPT algorithms are there to track maximum power point such as perturb and observe (P&amp;O), incremental conductance, constant voltage, constant current and

In this area researchers have developed what are now called as the Maximum Power Point Tracking (MPPT) algorithms. Mummadi Veerachary has given a detailed report on the use of

Bolcskei, “Blind Estimation of Symbol Timing and Carrier Frequency Offset in Wireless OFDM Systems,” IEEE Transactions on Communications, vol.. Swami, “Blind Frequency Offset

Voltage generated by cow dung using double chamber MFC was recorded at an interval of 1&amp;1/2 hr per day for the entire time period of 6 days as shown in Fig 4-9.. The

DISCRETE-TIME SLIP CONTROL ALGORITHMS FOR A HYBRID ELECTRIC VEHICLE 4 such as fuzzy logic control [6], neural network [7], hybrid control [8], adaptive control [9], and

4.6 Input image with Gaussian noise 0.05 corresponding image after VISU Shrink thresholding 4.7 Input image with high Gaussian noise and corresponding image after SURE Shrink..

Detailed JH NMR analysis (Table 2) revealed that the interaction between the --OH of calix(n)arene and &gt; C = O of benzoquinone is stronger at low

Chapter 3: Unsupervised segmentation of coloured textured images using Gaussian Markov random field model and Genetic algorithm This Chapter studies colour texture image