Implementation of Region-Mura Detection Based on Recursive
Polynomial-Surface Fitting Algorithm
WangZhengyao

MaLingFast (shanghai) Imaging Technology
C-101, No.456, Bibo Rd, Pu Dong
Shanghai,China

Fast Corporation
2791-5 Shimotsuruma,Yamato
Kanagawa,242-0001,Japan

ABSTRACTdefect.This paper focus on the region-Mura defects, see
Figure1.

Mura defect appear as low contrast, non-uniform
brightness regions in the Flat Panel Display(FPD). Automated Mura inspection method is needed in FPD production. There are many types of algorithms in Mura detection, such as algorithm based on singular value decomposition and surface fitting. Robustness and efficiency isimportant in the online FPD inspection systems. This paper proposes an efficient Mura detect algorithm for the
automatic inspection of FPD. The non-uniform brightness region can be detected by the difference of one fitted
Surface and the original image at every pixel. Surface
fitting at every pixel is computing cost. This paper proposes a recursive surface fitting method to improve inspection efficiency. Gauss filter is also used in order tosegment the non-uniform region. Performance and efficiency of the proposed algorithm has been evaluated onsome TFT(Thin FilmTransistor)-LCD panel samples. It
is showed that the algorithm is very effective.

Figure 1: Region-Mura in the FPD

KEY WORDS
Several algorithms have been developed for detection of region-Mura defects[1, 2, 3, 4, 5, 6, 7]. Since
the boundary between the regional Mura and the background is indistinct, multiple resolution analysis[2] isused to delete noise and segment the defect region for
low-contrast image. For defect detection of periodical,
repetitive textured patterns image, Lu and Tsai[3] developed an defect inspection algorithm based on singularvalue decomposition and reconstruction process. Lee andYoo[5] use surface fitting and regression to remove thenon-uniform background as the first phase to detect candidate region-Mura defects, in the second phase, basedon visual perception model, they quantify the candidate
Mura defects to identify real muras.
In FPD inspection process, needs of small defects detection on large size FPD demands high speed defectdetection algorithms. This paper presents an efficientmethod for the surface fitting scheme, making it possible for real time FPD automatic inspection systems.

Mura detection,FPD inspection,polynomial-surface fit

1

Introduction

Quality control is becoming a critical task in the FPD
manufacturing process. Surface defects on TFT panel
will cause visual failures. Cost and efficient production of FPD demand automatic inspection systems. At
present, human visual inspection remains the standardpractice for evaluating display luminance and color uniformity through all stages of display development andproduction. However, human inspection needs higher labor power and usually causes inefficient production andinconsistent inspection: individuals observe and interpretthings differently and a single person may report different results over time. Moreover it is difficult for human
to inspect large size FPD.
Automated inspection methods use imaging systems
combined with scientific digital cameras and sophisticated algorithms. High speed and high resolution digital camera makes automatic inspection possible. Automated defect detection algorithms have also improved tothe point that it can replace human visual inspection.
In general, Mura defect appear as low contrast, nonuniform brightness regions. They are larger than a single
pixel when the screen is driven to a constant gray level[1].
Depending on the shapes and sizes, Mura defects may
be classified into spot-Mura, line-Mura, and region-Mura

2 Mura Detection based on Surface Fitting
2.1

Basic Defect Detect Algorithm

Lee, Yoo and Choi[5, 6] proposed a surface-fitting
based region-Mura detection algorithm. For the sake of
performance and efficiency improvement, after some experiments, here we proposed a modified version algorithm. Let I(r, c), 1 ≤ r ≤ R, 1 ≤ c ≤ C is the observed
1

graylevel image of FPD use high resolution digital camera. Then the algorithm works as follows:

Use Matrix notion, let



f1
1 r1
 f2 
 1 r2
f =
 ...  , A =  ...
f9
1 r9I) Input: I(r, c), 1 ≤ r ≤ R, 1 ≤ c ≤ C, window size
W , polynomial order d, parameter of Gaussian filterσ, fixed threshold T ;thenII) Compute the difference of each pixel from the region surround it: for each pixel p = (r0 , c0 ) in theimage I except some marginal pixels;

f = Akk = (AT A)−1 AT f = Bfii) Fit a polynomial-surface model f (d) (r, c) in
the square window image use the least-squares
method(LSM),
X

αi,j ri cj(5)

Obviously, the least-squares solution of k isi) Get square window image centered at p, thesquare window size is W ;f (d) (r, c) =




c1
k1
c2 
 , k =  k2  (4)

k3
c9

(6)

where, A9×3 is known, 9 is the number of pixels and 3
is the number of polynomial terms. Then B3×9 is also
known and can be pre-computed and for each window
we can compute k1 , k2 , k3 using k = Bf .
Letb11 b12 b13 . . . b17 b18 b19
B =  b21 b22 b23 . . . b27 b28 b29 (7)b31 b32 b33 . . . b37 b38 b39(1)i≥0,j≥0,i+j≤diii) Compute the absolute value of difference ofdI(r0 , c0 ) = |I(r0 , c0 ) − f (d) (r0 , c0 )|

b1
b2
b3

III) Spread the difference of pixel to the region: D =dI ∗ g(r, c; σ), where ∗ means linear filtering andg(r, c; σ) is Gaussian filter; then D is the regionMura measures;= [b11 , b12 , b13 , b14 , b15 , b16 , b17 , b18 , b19 ] (8)= [b21 , b22 , b23 , b24 , b25 , b26 , b27 , b28 , b29 ] (9)= [b31 , b32 , b33 , b34 , b35 , b36 , b37 , b38 , b39 ] (10)

According to (6),we can get
IV) Construct a binary image use a single fixed threshold
T acquired from experiments.

k1 = b11 f1 + b12 f2 + ... + b19 f9 = b1 ∗ f
k2 = b21 f1 + b22 f2 + ... + b29 f9 = b2 ∗ f
k3 = b31 f1 + b32 f2 + ... + b39 f9 = b3 ∗ f

As described in §2.5, experiments showed that d = 3,
W = 101 and σ = 3.0 is good for the most images.
There are some key problems in this algorithm:
polynomial-surface model f (d) (r, c) fit, Gaussian filtering and single fixed threshold T acquired from experiments. They will be described in the following subsections.

(11)
(12)
(13)

that is, k1 , k2 , k3 can be computed from linear convolution and the convolution templates, corresponding to eachrows of B3×9 , can be pre-computed. This algorithm canbe generalized for any W and any d.
So this method transform the polynomial-surface fitting problem to template convolution. Noticed that only
k1 is necessary in computing dI if local coordinate system, like (2), is used. Obviously, the result of this method
is equivalent to the direct numerical method but the speed
is faster than the direct. But this method also very slow
for large W which is needed in region-Mura detection,
the computing cost of the method is too high to use.

2.2 Two Types Surface Fitting Algorithms
To fit the polynomial-surface model f (d) (x, y) use
LSM, the basic algorithm is direct numerical method as
described in[8]. Suppose M = N = 800 and W = 101,
there will be 700 × 700 = 490, 000 times such fits andeach fit use 101 ∗ 101 = 10, 201 pixels to compute themodel. It is terrible.
The algorithm can be improved by convolution. For
the sake of simplicity, here we give a simple example:
W = 3 and d = 1, that is try to fit f (r, c) = k1 + k2 r +k3 c use 9 pixels. Use local coordinate system to describe
the pixels in the window,(−1, −1) (−1, 0) (−1, 1)
 (0, −1)(0, 0) (0, 1) 
(2)
(1, −1)(1, 0) (1, 1)

2.3

Recursive Polynomial-Surface Fitting

In order to get much faster algorithm to fit apolynomial-surface, Haralick, Zuniga and Qiang Ji[9, 10,11] have proposed many algorithms and the best one is
the recursive algorithm.
Def. 1 Symmetric Index Set:R is a symmetric index
set, if n ∈ R implies −n ∈ R, let the number of elements
in R be N .

corresponding grayvalue are f1 , f2 , ..., f9 . Then we have,


 

1 r1 c1 f1
 f2   1 r2 c2  k1
  k2 

 
(3)

 ...  =  ...
k3
1 r9 c9
f9

Use Pn (r) express n-order polynomial with the coefficient of the highest ordered item is 1, here call it nordered polynomial for short. General speaking,Pn (r) = rn + an−1 rn−1 + ... + a1 r + a0
2

(14)

+K9 r(c2 − a) + K10 (c3 − bc)Def. 2 Discrete Orthogonal Polynomial on Symmetric Index Set:{P0 (r) = 1, P1 (r), ..., PN (r)} is a groupof polynomials, which have N + 1 elements, with different degrees, where N is the number of elements of thesymmetric index set. It is called discrete orthogonal polynomials on symmetric index set R if
X
Pm Pn = 0, 0 ≤ m < n ≤ N
(15)

=

(23)

K1 , ..., K10 and k1 , ..., k10 have following relations:
k1
k2
k3
ki

Haralick[9]proposed a recursive formula of discrete
orthogonal polynomials on symmetric index set,

where,

Ki gi (r, c)= 1= rPi (r) − βi Pi−1 (r)

=
=
=
=

K1 − aK4 − aK6
K2 − bK7 − aK9
K3 − bk10 − aK8
Ki , i = 4, ..., 10r∈R Pi (r)Pi−1 (r)P
2r∈R Pi−1 (r)

(24)
(25)
(26)
(27)

Let
G(r, c) =10
X

Ki gi (r, c)

(28)

i=1

P
βi =

(22)

i=1

r∈R

P0 (r)
Pi+1 (r)10X

so fit the polynomial-surface model f (d) (x, y) use LSM
is equivalent to the following problem: minimize e about
Ki , i = 1, ..., 10,where(16)

So we can get arbitrary degree of discrete orthogonal
polynomial on R using this recursive formulas. In our
problem, n = 3 is enough andXe=[I(r, c) − G(r, c)]2(29)(r,c)∈R×C

P0 (r) = 1(17)
P1 (r) = r(18)
2P2 (r) = r − a(19)
3P3 (r) = r − br
(20)
P
µ4
k
where, a = µµ20 , b = µ2
, µk =r∈R s .When R ={−2, −1, 0, 1, 2}, µ0 = 5, µ2 = 10, µ3 = 17, a = 2, b =3.4.
The two dimensional discrete orthogonal polynomial can be constructed from the tensor product of the two sets of one dimensional discrete
polynomials.
Let {P0 (r), P1 (r), ..., PN (r)} and{Q0 (c), Q1 (c), ..., QM (c)} be discrete orthogonal polynomials on the symmetric index setR and C respectively.
Then {P0 (r)Q0 (c), ..., Pm (r)Qn (c), ..., PN (r)QM (c)}is the discrete orthogonal polynomials on the 2Dsymmetric index set R × C.
When R = C = {−2, −1, 0, 1, 2}, the tensor productdiscrete orthogonal polynomial on R × C is {1, r, c, r2 −a, rc, c2 − a, r3 − br, (r2 − a)c, r(c2 − a), c3 − br, (r3 −bc)c, (r2 −a)(c2 −a), r(c3 −bc), (r3 −br)(c2 −a), (r2 −2)(c3 − bc), (r3 − br)(c3 − bc)}. Because we only carepolynomials with order not more than d = 3, so weshould make a selection from above set. Tensor productdiscrete orthogonal polynomials with degree not exceed
d = 3 on R×C are {1, r, c, r2 −a, rc, c2 −a, r3 −br, (r2 −a)c, r(c2 − a), c3 − bc}, we use gi (r, c), i = 1, ..., 10 toexpress them.
Obviously, any bivariate cubic polynomial can be expressed by linear combination of those discrete orthogonal polynomials,

Because gi (r, c) and gj (r, c), i 6= j is orthogonal,PKi =(r,c)∈R×C gi (r, c)I(r, c)P2(r,c)∈R×C gi (r, c)(30)(30) tell us that each coefficient Ki can be computed using I(r, c) by template convolution independently and thecorresponding template is Wi ,gi (r, c)2(r,c)∈R×C gi (r, c)

Wi = P(31)

Compute Ki using Wi is possible but not the best. Because gi (r, c) is separable,
g1 (r, c)g2 (r, c)g3 (r, c)g4 (r, c)g5 (r, c)g6 (r, c)g7 (r, c)g8 (r, c)g9 (r, c)g10 (r, c)

=
=
=
=
=
=
=
=
=
=

1=1×1
r =r×1c=1×c
r2 − a = [r2 − a] × 1
rc = r × cc2 − a = 1 × [c2 − a]r3 − br = [r3 − br] × 1(r2 − a)c = [r2 − a] × cr(c2 − a) = r × [c2 − a]c3 − bc = 1 × [c3 − bc]

so 2D template convolution can be decomposed intosome simple 1D template convolution. For example, template convolution using g1 (r, c) can be decomposed into two 1D template convolution: first template convolution using row template[1, 1, 1, ..., 1, 1, 1]
and then using column template [1, 1, 1, ..., 1, 1, 1], here
template[1, 1, 1, ..., 1, 1, 1] is generated by polynomial 1.
Now we only focus on some simple 1D template convolution. Suppose the input is I1 , I2 , ..., IN , we are only

f (r, c)= k1 + k2 r + k3 c + k4 r2 + k5 rc + k6 c2+k7 r3 + k8 r2 c + k9 rc2 + k10 c3(21)2= K1 + K2 r + K3 c + K4 (r − a) + K5 rc+K6 (c2 − a) + K7 (r3 − b) + K8 (r2 − a)c3interested in following 4 types template convolution:
xn

=

K
X

In+r

IV) Filtering the image column by column:(32)i) Filter x using the same method as that in III)to get K1 , K3 , K6 , K10 ;(33)ii) Filter y using the same method as that in III) to
get K2 , K5 , K6 , K9 ;

r=−K

yn

=

K
X

rIn+r

iii) Filter z using the same method as that in III) to
get K4 , K8 ;

r=−K

zn

=

K
X

(r2 − a)In+r(34)iv) Filter v using the same method as that in III) to
get K7 ;(r3 − br)In+r(35)v) Compute ki from Ki using (24), (25), (26) and(27) respectively.

r=−K

vn

=

K
X
r=−K

All filtering above use polynomial as filter,
f ilter1 (s)f ilter2 (s)f ilter3 (s)f ilter4 (s)

=
=
=
=

2.4

1
s
s2 − as3 − bs

Now we will construct the recursive formulas of (32),(33), (34), (35).
xn+1 =K
X

In+1+r = xn + In+K+1 − In−K

(36)

r=−K

yn+1 = yn − xn+1 + (K + 1)In+K+1 + KIn−K (37)zn+1

= zn − 2yn+1 − xn+1+[(K + 1)2 − a)]In+K+1−[K 2 − a]In−K

(38)

2.5
vn+1

Gaussian filtering and Mura Region Label

After difference surface fitting image with original image, we can segment Mura defects from the image. Before labeling mura area, Gauss filter can be used to remove the noise.
Gaussian filter is one of the important filters in the image processing and computer vision and then many researchers have studied their properties and implementation. This paper have implemented some of them, and theconclusion is: Deriche’method[12] is a good implementation for parallel computer, Demnigy’s method[13] is agood implementation for real-time hardware and Young’s
method[14] is a good implementation for serial computer.
We can choose appropriate method in different system.
The best way to label Mura in the region-Mura measures D is the method based on human visual perception
like the ones in Lee, Yoo and Choi[5, 6] and Chen[15].
Even then, this paper try to choose a fixed threshold T ,based on human visual perception experiments, for all
images.

Experiment and Parameters Selection

Degree of polynomials, d, is important in surface fitting. To get more accurate approximation, higher degree
is wanted, but Approximation Theory told that polynomial with degree more than 4 or 5 is improper because
of inflexibility. So d = 3, 4, 5 is used in most applications. In image processing, d = 3 is enough, so this paperchoose d = 3.
W is the most important parameter in this program.
This paper did many experiments, see Figure2, and foundthat W = 101 is the best.
Gaussian smooth is necessary. Figure3 showed that
serious texture exists in image and dI is hard to process.
In Figure3, (c) is the smoothed image of (a), it is clearlyshown that smooth can spread the difference of pixel tothe region and eliminate isolated noise which is createdby texture. Experiments showed that σ = 3.0 is best.
The last parameter is T . Some volunteers were asked
to observe the images and the results. This paper first use
T = 50.0 to detect Mura, there is noting in detected image so volunteers dissatisfied the result, then this paper
decrease T slowly and detect Mura again. This process
was repeated again and again until they satisfied the results. T is different for different volunteers, and the average is T = 0.35. So this paper choose T = 0.35 as thethreshold.

= vn − 3(zn − yn ) − (3a − b + 1)xn+[K 3 − bK]In+K+1+[(K + 1)3 − b(K + 1)]In−K(39)

The formula we proposed here, (39), is different fromcorresponding formula in Qiang Ji and Haralick’s[11].
After carefully examination and numerical experiments,
we see that there are some misprints in their paper.
Using (36), (37), (38), (39), we can compute ki , i =1, ..., 10 of recursively:I) Input: I(r, c), 1 ≤ r ≤ R, 1 ≤ c ≤ C, K = (W −1)/2;II) Open float buffer for x, y, z, v, each have the same
size of image;III) Filtering the image row by row:i) for the Kth pixel in each row, using (32), (33),(34), (35)to compute value of corresponding
position in x, y, z, v respectively;ii) for other pixels in each row, using (36), (37),(38), (39) to compute value of corresponding
position in x, y, z, v respectively;
4

There are some examples of region-Mura detection,with d = 3, W = 101, σ = 3.0 and T = 0.35, in Figure
4.

(a) W = 11(b) W = 41(c) W = 61(d) W = 81(e) W = 101(f) W = 121

Figure 4: Example of Region-Mura Detection.

Figure 2: W = 101 is good. dI of the image in Figure
1 with various W , d = 3. To show the dI clearly, we
multiply dI by 10.

Region-Mura may appears in two types of images:
with or without periodical, repetitive textured patterns.
Experiments showed that this algorithm can be used todetect region-Mura in both of them.
Everyone can download the executable file of the program from my homepage[16].

2.6Comparison

For the same parameters, the performance of recursive method is same as that of template convolution based
method, see Figure 5.
Table 1 show that the speed of recursive based implementation is much faster than that of template convolution based.
(a) dI of Fig.1(b) Binarized image of dI

Table 1: Comparison of elapsed time of proposed regionMura detect algorithm
Implementation
Template based
Recursive based(c) D = Smoothed dI

Average Elapsed time
(Unit:MSEL)
12200
65

Machine:Intel P4 CPU 3.0GHz;512M Memory
OS:Windows XP sp2;
Program Language of algorithm:Standard C
Compiler: Visual C++ 6.0
Image Size:832 × 621(d) Binarized image of D

Figure 3: Smooth is necessary. To show the dI clearly,
we multiply dI by 10.

5[5] Jae-Young Lee and Suk-in Yoo,”Automatic detectionof region-mura defect in TFT-LCD”, IEICE Transactions on Information and Systems, vol.E87-D, no.10,
pp.2371-2378, Oct. 2004.
[6] Kyu-Nam Choi and Suk-in Yoo, ”Area-Mura detection in TFT-LCD panel”, Electronice Imaging,
pp.151-158, Jan. 2004.

(a) Example image[7] HC Chen, LT Fang, L. Lee, CH Wen, SY Cheng,
and SJ Wang,”LOG Filter Based Inspection of Cluster Mura and Vertical Band Mura on Liquid Crystal Displays”, Proceedings of SPIE, Vol. 5679, Jan.
2005.

(b) dI based on template based[8] William H. Press, Saul A. Teukolsky, William T. Vettering, and Brian P Flannery, ”Numerical Recipes inC++: The Art of Scientific Computing(Chapter 2)”,
second edition, Cambridge University Press, 2002.

(c) Final result of (b)[9] Robert M.Haralick,”Digital Step Edges from Zero
Crossing of Second Directional Derivatives,” IEEE
Trans. Pattern Analysis and Machine Intelligence,
Vol.6, No.1, pp.58-68, JANUARY 1984(d) dI based on recursive based[10] O.A. Zuniga and Robert M. Haralick , ”Integrated Directional Derivative Gradient Operator,”
IEEE Trans. Systems, Man, and Cybernetics, vol.
17,no. 3,; pp.508-517,1987.

(e) Final result of (d)Figure 5: Performance comparison between recursive
based method and convolution based method

3[11] Qiang Ji and Robert M. Haralick,”Effcient facetedge detection and quantitative performance evaluation”,Pattern Recognition 35, pp.689-700.2002.

Conclusions[12] R. Deriche, Fast Algorithms for Low-Level Vision,IEEE Transactions on PAMI,Vol.12, No.1,
pp.:78-87, January, 1990.

This paper proposed a modified version fitting based
region-Mura detect algorithm and use recursive polynomial fitting and recursive Gaussian filtering techniques toimplement it. The result is good for many images and theprogram is fast enough in online Mura detection system.

4[13] Didier Demigny, L. Kessal, J. Pons: Fast Recursive
Implementation of the Gaussian Filter. VLSI-SOC
2001: 39-49.

Acknowledgments

Thanks to our colleagues, especially LiuWei and JiaBo, who help and encourage us to complete this paper.

[14] IT Young and LJ van Vliet, ”Recursive implementation of the Gaussian filter,” Signal Process., vol. 44,
pp. 139-151, 1995.

References[1] W.K.Pratt,S.S.Sawkar and K.O’Reilly, ”Automatic blemish detection in liquid crystal flat panel
displays”, IS&T/SPIE Symposium on Electronic
Imaging:
Science and Technilogy,3306, San
Jose,California,January,1998.

[15] HC Chen, LT Fang, L. Lee, CH Wen, SY Cheng,
and SJ Wang, LOG Filter Based Inspection of Cluster Mura and Vertical Band Mura on Liquid Crystal
Displays Proceedings of SPIE, Vol. 5679, Jan. 2005.
[16] Download the executable file of the program.
http://www.quzhi.net/chinese/vision/mura/rmura.htm.

[2] Hiroki Nakano, Yumi Mori, ”Measurement method
for low-contrast nonuniformity in liquid crystal displays by using multi-wavelet analysis”, Proc.SPIE
Vol. 5880, p. 313-318, Optical Diagnostics; Leonard
M. Hanssen, Patrick V. Farrell; Eds 2005.
[3] Chi-Jie Lu and Du-Ming Tsai, ”Automatic Defect Inspection for LCDs Using Singular Value Decomposition”, Int. J. Adv Manuf Techol(2004).
[4] Chih-Chung Chen, ”A Vision System for Dot Defect and Uniformity Automatic Inspection on LCD”,
Master’s Thesis,Mechanical Engineering, 2001.
6