OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 18, Iss. 5 — Mar. 1, 2010
  • pp: 4148–4157
« Show journal navigation

Simulation study on the detection of size, shape and position of three different scatterers using Non-standard time domain time inverse Maxwell’s algorithm

Kisalaya Chakrabarti and James B. Cole  »View Author Affiliations


Optics Express, Vol. 18, Issue 5, pp. 4148-4157 (2010)
http://dx.doi.org/10.1364/OE.18.004148


View Full Text Article

Acrobat PDF (199 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Inverse method has wide application on medical diagnosis where non-destructive evaluation is the key factor .Back scattered waves or echoes generated from the forward moving waves has information about its geometry, size and location. In this paper we have investigated how well different geometries of the object is determined from the back scattered waves by a high accuracy Non-Standard Finite Difference Time Inverse (NSFD-TI) Maxwell’s algorithm and how the refractive index of the object plays a deterministic role on its size.

© 2010 OSA

1. Introduction

In medical diagnosis, detection of malignant tumors of various size and shape from back scattered waves [1

1. L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).

] with the application of the Inverse theorem is widespread. The challenge lies how well this method can detect the tumors when the refractive index of the malignant cells and the normal cells differs by a small quantity. In this paper we have investigated three different reconstructive spectrums of three different shapes of scatterer (circle, rectangle and triangle) using NSFD-TI algorithm based on Non-Standard Finite Difference [2

2. R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).

,4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

] and also we have observed the effect of refractive index upon the size of the circular scatterer.

2. Model and Calculation Methods

Maxwell’s equations for a non-conducting dielectric are given by:
tB(x,t)=×E(x,t)
(1)
tD(x,t)=×H(x,t)J(x,t)
(2)
·B(x,t)=0
(3)
·D(x,t)=ρ(x)
(4)
where, H is the magnetic field, B is the magnetic induction, E the electric field, D the electric displacement, J the electric current density, and ρ is the density of free charge. Materials, such as iron or glass, modify the electric and magnetic fields. B describes how the material affects the magnetic field, and D describes how it affects electric field. The vector position is x = (x,y,z), t=/t, and =(x,y,z) is a vector differential operator. The curl of a vector is defined by×V=(yVzzVy,zVxxVz,xVyyVz), and the divergence of V is·V=xVx+yVy+zVz, whereV=(Vx,Vy,Vz).

In most optical materials there is no free charge andρ=0. In a non-conducting dielectric, such as glass, no electric current flow so J = 0. Also to an excellent approximation
B(x,t)=μH(x,t)
(5)
where, μ, the magnetic permeability, is a constant. In general, the electrical permittivity, ε depends upon the frequency, but when this dependence is weak, it is a good approximation to take
D(x,t)=ε(x)E(x,t)
(6)
The permittivity is different for different materials, however so it is a function of position. For example, in glassε1.2. Assuming that Eq. (5) and Eq. (6) hold, Maxwell’s equations for a non-magnetic, non-conducting medium with no free charges or electrical currents is
μtH(x,t)=×E(x,t)
(7)
ε(x)tE(x,t)=×H(x,t)
(8)
Now, let us introduce what is called a Finite Difference (FD) approximation. The first derivative is approximately given by
f(x)f(x+h2)f(xh2)h
(9)
For convenience let us define the difference operator
dxf(x)=f(x+h2)f(xh2)
(10)
with this definition the FD approximation to f is
f(x)dxhf(x)
(11)
Defining the vector difference operator:
d=(dx,dy,dz)
(12)
and,
f(x)dhf(x)
(13)
where, Δx=Δy=Δz=h. The second derivativef can be approximated by applying Eq. (11) twice:
f(x)dx2h2f(x).
(14)
and from the definition of dx it is easy to show that dx2=dxdx is given by:
dx2f(x)=f(x+h)+f(xh)2f(x)
(15)
Defining,
d2=dx2+dy2+dz2
(16)
We now have,
2f(x)d2h2f(x)
(17)
and also replacing the derivatives of Maxwell’s Eq. (7) by FD approximations we obtain
dtH(x,t)=1μΔthd×E(x,t)
(18)
Here the electromagnetic field components are arranged on the grid in such a way those central FDs can be taken. ExpandingdtH(x,t), and solving for H(x,t+Δt2) gives
H(x,t+Δt2)=H(x,tΔt2)1μΔthd×E(x,t)
(19)
Doing the same for Eq. (8) we have
dtE(x,t+Δt2)=1ε(x)Δthd×H(x,t+Δt2)
(20)
and we computedtE(x,t+Δt2), rather than dtE(x,t) in order to use the updated value of H. ExpandingdtE(x,t+Δt2), and solving for E(x,t+Δt) gives
E(x,t+Δt)=E(x,t)1ε(x)Δthd×H(x,t+Δt2)
(21)
Equation (19) and Eq. (21) are called the Yee algorithm [3

3. K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).

]. The notation does not show it, but each component of H and E is evaluated at a different position. For example to computedtHzdxEydyEx, if Hz lies on (x,y,z) then Ey must be known at(x±h2,y,z), and Ex at(x,y±h2,z).

To develop the Nonstandard-FDTD (NS-FDTD) model or in short NSFD we first construct an alternative Standard-FDTD (S-FDTD) model or SFD of Maxwell’s equations. DefiningH=μhH/Δt, we obtain
dtH(x,t)=d1×E(x,t),
(22)
dtE(x,t+Δt2)=1εμΔt2h2d1×H(x,t+Δt2)
(23)
We introduced in earlier publication [4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

,7

7. J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002). [CrossRef]

], the vector difference operator d0 which is also highly isotropic in nature with the property that
d1·d0=d0·d1=d02.
(24)
whereas, d1·d1=d02,d0·d0d02. In two dimensionsd0=(dx(0),dy(0)), where
dx(0)=dx+(1γ4)dxdy2
(25)
dy(0)=dy+(1γ4)dydx2
(26)
andd0=d1(1+1γ4dxdy)
(27)
Here, γ which is also a function of k=k0ε (assuming uniform μ), given by the expression γ=23190(kh)2115120(kh)4(1152) is chosen to minimize the error of the NSFD approximation [4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

].

Making the substitutions d1d0 and Δt2/(h2εiμ)u0(ε)2 in Eq. (23) we obtain a new NSFD model of Maxwell’s equations,
dtH(x,t)=d1×E(x,t)
(28)
dtE(x,t+Δt2)=u02εd0(ε)×H(x,t+Δt2)
(29)
Solving Eq. (28) and Eq. (29) for H(x,t+Δt2) andE(x,t+Δt), we obtain the NS-Yee algorithm
H(x,t+Δt2)=H(x,tΔt2)d1×E(x,t)
(30)
E(x,t+Δt)=E(x,t)+u02εd0(ε)×H(x,t+Δt2)
(31)
This NSFD algorithm is somewhat different from the one introduced in [4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

]. We find that it gives better results when μ is constant. To reverse the time direction we need to reverse the sequence of the above computations; doing so we get,
H(x,tΔt2)=H(x,t+Δt2)+d1×E(x,t+Δt)
(32)
and
E(x,t)=E(x,t+Δt)u02εd0(ε)×H(x,tΔt2)
(33)
Here, we observed that we basically get the same equations as we already derived for forward propagating wave. This depicts that even if we iterate the wave in time reverse direction the nature of the wave remain same. For this purpose we have applied the periodic boundary conditions instead of Mur boundary conditions because the later causes field explosion when we iterate in time reverse condition [5

5. R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993). [CrossRef]

]. Let, the computational domain bex=0,h,2h,Nxh. In periodic boundary we define, ψ(h,t)=ψ(Nxh,t)where, waves that exit the computational domain at x=0 re-enter it fromx=Nxh.

Algorithm Stability

When any algorithm is iterated is can be become unstable. Let A be an algorithm to solve for an unknown vector ψ. Let φ0 be the initial solution vector. Thenφ1=Aφ0, φ2=Aφ1, ⋯, and thus φn=Anφ0. If lim|φn|n when |ψ| is finite, the algorithm is said to be unstable.

It can be shown that the D-dimensional FDTD algorithm is stable ifvΔthDD. Taking finite-difference algorithm of the form:
ψ(t+Δt)=pψ(tΔt)+2qψ(t)
(34)
Defining t=nΔtand ψn=ψ(nΔt)we get,
ψn+1=pψn1+2qψn
(35)
Now we postulate a solution:Inserting ψn=λn and dividing byλn1 we get,
λ22qλp=0
(36)
Solving for λ we get,
λ±=q±q2+p
(37)
Thus the most general solution is:
ψn=c+λ+n+cλn
(38)
The condition for the stability will be|ψn| being finite as n is|λ±|1.

We want finite oscillatory solutions of the difference Eq. (36), for the oscillation the solutions must have an imaginary part. If p<0p=|p| and p+q2<0 and λ±=q±q2|p|=q±i|p|q2|λ±|2=q2|p|+q2=|p| . Thus |p|1 and q2<|p| are the conditions for stability for an oscillating solution.

Applying these results to the wave equation given below:
(t2v22)ψ(x,t)=0,
(39)
Generalized Finite Difference model is:
(dt2u02d2)ψ(x,t)=0
(40)
and the FDTD algorithm is
ψ(x,t+Δt)=ψ(x,tΔt)+(2+u02d2)ψ(x,t)
(41)
While in S-FDTD algorithm, d2=d12 andu02=v2Δt2h2, whereas in NS-FDTD d2=d02 andu02=sin2(ωΔt/2)sin2(kh/2). To direct the FDTD algorithm in the form of (34) we need to computed2ψ. Substituting ψ=φk=ei(k.xωt) we can write, d2φk=2d2φk(x,t) where the values of d2 vary case by case [4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

]:

For SFD:
d2=d12=2sin2(kh/2): 1Dimension
d2=d12=2[sin2(kxh/2)+sin2(kyh/2)]: 2Dimensions
Max(d12)=2Dwhere, D = Dimension

For NSFD:
d02=d12=dx2and Max(d02)=2: 1Dimension
also,
d22=1cos(kxh/2)cos(kxh/2): 2Dimensions
Till now we have not specified the value ofγ. For the stability we optimize the value of γ by γ0 (also termed as weight coefficient ford12andd22)Where
γ0=1sin2(kxh/2)+sin2(kyh/2)sin2(kh/2)2sin2(kxh/2)sin2(kyh/2)
(42)
and(kx,ky)=k(21/4,121/2).Thus for the NSFD algorithm we have,
d02=γ0d12+(1γ0)d22
(43)
From the expanded value of γ given above, we get Max (γ0) = 2/3. The above equation can be rewritten as:
d02(k)=γ0d12(k)+(1γ0)d22(k)
(44)
where, k=k(cosθ,sinθ)and Max(d02)=8/3 in 2 Dimensions.If we write, ψ(x,t)=ψxt then the FDTD algorithm (41) becomes:
φxt+1=φxt1+2(1d2u02)φxt
(45)
This equation is in the form of (36) and the solutions are:
λ±=1u02d2±(1u02d2)21
(46)
To accomplish the stability condition|λ±|1, we have to have(1u02d2)21.

If d2>0andu02>0, for non-vanishing monochromatic plane waves, the stability condition turns out to be [4

4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

]:
u022d2
(47)
Whiled2is a function of the magnitude and direction of k, stability is guaranteed by inserting the maximum value ofd2, Max(d2) instead ofd2.
u022Max(d2)
(48)
Considering NSFD algorithm, we have d02=d12=d2and Max(d02)=2,8/3for dimensions, D = 1, 2 respectively. Defining earlier, u02=sin2(ωΔt/2)sin2(kh/2)the stability condition (48) becomes,
sin(ωΔt/2)sin(kh/2)2Max(d02)
(49)
Replacing, ϖ=ωΔt,k¯=khandϖ=ck¯, where the value of c is to be found out so that (49) does not violate. In D dimensions expressing Max(d02)=mdnsfd the Eq. (49) turns out to be:
sin(ck¯/2)2mdnsfdsin(k¯/2)
(50)
In the range0θπ/2, sin(θ)rises monotonically with increasing θ, so it is obvious thatsinθ1sinθ2θ1θ2. The Nyquist limit can be articulated in the formk<π/D. Hence we shall never figure out a FDTD algorithm which disobeys λh>2D condition, and presuming that, c1 we conclude from (50) that
c<2dπarcsin[2mdnsfdsin(π2D)]
(51)
Putting the numerical values of mdnsfdfrom above, we have

c1nsfd=1: In One Dimension
(52)
c1nsfd22πarcsin[32sin(π22)]0.84: In Two Dimensions
(53)

3. Verifications and Practical Tests

In Fig. 2(a)
Fig. 2 Angular distribution of scattered intensity about a circular contour of radius λ0 centered on the axis of the cylinder. (a) Analytic solution (black) compared with the NSFD calculation (NSFD-8, red) and the SFD one (SFD-8, blue) at λ0/h=8. (b) Analytic solution (black) compared with the SFD (SFD-24, blue) one at λ0/h=24.
the scattered fields shown in Fig. 1 are plotted as a function of angle the (θ) from the +x-axis on a circular contour of radius λ0 centered on the axis of the cylinder. At discretization λ0/h=8 the root mean square error for the NS-FDTD and S-FDTD algorithms is εNS-8=0.04, and εS-8=0.20, respectively. In Fig. 2(b) at discretization λ0/h=24 we compare the S-FDTD calculation with the analytic solution. The root mean square error is εS-24=0.04 the same as the NS-FDTD algorithm atλ0/h=8. At λ0/h=8 the NS-FDTD algorithm delivers the same accuracy as the S-FDTD one does at λ0/h=24. Because of the lower coarser discretization, the computational cost of NSFDTD algorithm is only 1/27th that of the Standard FDTD one.

Understanding the fact that NSFD algorithm is better than S-FDTD one, we extended our study for the Time Inverse algorithm on the detection of size, shape of different scatterers. for this purpose we have generated a weighted source array continuous wave Gaussian beam that turns on at the time t = 1. The beam width of the Gaussian pulse taken for the simulation purpose is 1.25λ0, where λ0 is the wavelength of the central frequency of the pulse. For the simulation we have chosen three types of the scatterer; Circular, Rectangular and triangular. Now we have allowed the light to be passed through the scatterer and it further moves to a distance which is very close to the periodic boundary of the computational space (Nx=140h), then we remove the scatterer and iterate the forward moving wave in the time reverse in order to trace the location and shape of the scatterer and examine how well it will match with the actual location and size. For this experiment, we have confined our observation on TE mode only, i.e. we will compute x-intensity and y-intensity as depicted in Fig. 3
Fig. 3 Time Inverse Computed Intensities using NS-FDTD for: (a)Circular scatterer (diameter = 2λ0)(b) Rectangular scatterer (side = 2λ0) (c)Triangular scatterer(base = 2λ0,height = 1λ0),where, λ0/h=10
.

4. Results and discussions

From the Fig. 3 we noticed that the Time Inverse algorithm (TE mode) has two intensities, x-intensity and y intensity in which left one is the computed x-intensity and right one is the computed y-intensity. It is clearly seen from the figure that x-intensities gives the clear indication of the shape, which means the first one is circular, the middle one is the rectangular and the last one is the triangular scatterer. The y- intensity is perpendicularly oriented to the x-intensity and it does not differentiate the edges well as it does for the x-intensity. But considering both the intensities we can distinguish the exact shape, size and location of the scatterer by Non standard Time Domain Time Inverse algorithm. It is very helpful in medical imaging because a small cancerous tumor of arbitrary shape is perfectly detected and chance of misdetection or false detection is eliminated.

In Fig. 4, shown below applying the TM Mode of polarization of incoming light we have noticed that as the refractive index of the circular scatterer gradually increases the distinguishable size of the scatterer gradually decreases. One of the probable reasons is that, the scatterer having the refractive index close to vacuum, it size must be large enough in order to get some considerable contrast. We have also noticed that our Non-Standard Time Domain Inverse algorithm is good enough to detect the shape, size and location of the object in comparison to Standard Time Domain Inverse algorithm even if the refractive index is very close to unity. In the later case more time steps are required and the picture contrast is not good as that of NSFD algorithm.

Acknowledgements

We would like to acknowledge the Ministry of Education, Government of Japan for granting one of the authors (KC) with a scholarship during the period of this work.

References and links

1.

L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).

2.

R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).

3.

K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).

4.

J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.

5.

R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993). [CrossRef]

6.

P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).

7.

J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002). [CrossRef]

8.

H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART.

OCIS Codes
(200.0200) Optics in computing : Optics in computing
(070.7345) Fourier optics and signal processing : Wave propagation

ToC Category:
Scattering

History
Original Manuscript: October 7, 2009
Revised Manuscript: December 14, 2009
Manuscript Accepted: December 16, 2009
Published: February 17, 2010

Virtual Issues
Vol. 5, Iss. 6 Virtual Journal for Biomedical Optics

Citation
Kisalaya Chakrabarti and James B. Cole, "Simulation study on the detection of size, shape and position of three different scatterers using Non-standard time domain time inverse Maxwell’s algorithm," Opt. Express 18, 4148-4157 (2010)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-5-4148


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. L. W. Schmerr, Jr., in Fundamentals of Ultrasonic Nondestructive Evaluation: A Modeling Approach, (Plenum Press, New York, 1998).
  2. R. Mickens, E, in Nonstandard finite difference models of differential equations, (World Scientific Publishing Co., Inc., River Edge, NJ, 1994).
  3. K. S. Kunz, and R. J. Luebbers, in The Finite Difference Time Domain Method for Electromagnetics, (CRC Press, New York, 1993).
  4. J. B. Cole, S. Banerjee, and M. Haftel, in Advances in the Applications of Nonstandard Finite Difference Schemes, ed. R. E. Mickens (World Scientific Singapore, 2007), Chap.4.
  5. R. Sorrentino, L. Roselli, and P. Mezzanotte, “Time reversal in finite difference time domain method,” IEEE Microw. Guid. Wave Lett. 3(11), 402–404 (1993). [CrossRef]
  6. P. W. Barber, and S. C. Hill, in Light Scattering by Particles: Computational Methods, (World Scientific, Singapore, 1990).
  7. J. B. Cole, “High-Accuracy Yee algorithm Based on Nonstandard Finite Difference: New Developments and Verifications,” IEEE Trans. Antenn. Propag. 50(9), 1185–1191 (2002). [CrossRef]
  8. H. Kudo, et al., “Numerical Dispersion and Stability Condition of the Nonstandard FDTD Method” Electronics and Communications in Japan, 85, 22–30(2002), http://www3.interscience.wiley.com/cgi-bin/fulltext/93514073/PDFSTART .

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1 Fig. 2 Fig. 3
 
Fig. 4
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited