Method of calculating shading correction coefficients of imaging systems from non-uniform and unknown calibration standards

[0001] This application claims the benefit of U.S. Provisional Patent application Ser. No. 60/236,945, filed Sep. 28, 2000 which is incorporated herein by this reference.

[0002] The present invention relates to a process for analyzing shading errors in an imaging system and determining shading correction coefficients to compensate for non-uniformity in illumination, optics and imaging sensors and electronics.

[0003] Images acquired by an imaging system are not perfect. One common type of distortion is shading error (systematic intensity variation across acquired images), which is mainly caused by:

• [0004] i. Photo-electronic conversion (sensors) inhomogenity;
• [0005] ii. Illumination inhomogenity;
• [0006] iii. Optical inhomogenity;
• [0007] iv. Electronic inhomogenity.

[0008] If not corrected, shading error will distort quantification performed from image intensities. It will also distort the visual appearance of the acquired images.

[0009] In the prior art, shading errors and shading correction coefficients of imaging systems are calculated from the image of a calibration standard (or reference) with uniform reflectance, transmittance, or emission across the entire field of view. Typically, the materials in a calibration standard are the same or have the same characteristics of the materials to be imaged. For example, to correct a fluorescent imaging system, a slide covered by a thin layer of fluorescent material with uniform thickness and density may be used as a calibration standard. The fluorescent substance in the calibration slide should have the same excitation/emission properties as those to be imaged and/or analyzed subsequently.

[0010] Typically, an imaging system is used to image and/or to analyze objects in a continuous fashion. That is, every point in the field of view is of interest. In this case, pixel-by-pixel shading corrections should be performed. FIG. 1 illustrates the mapping between object plane and image sensor plane. The rectangular area 200 in the object plane of an imaging system 250 is projected to the image sensor 202 of the imaging system 250. The area 200 is also referred to as the field of view of the system 250. In the example shown in FIG. 1, the image sensor 202 consists of 48 sensor elements 206 (or pixels) in an 8×6 array, denoted as pxl(x, y), x=1, . . . , 8 and y=1, . . . , 6. Correspondingly, there are 8×6 divisions 205, denoted as Rp(x, y), x=1, . . . , 8 and y=1, . . . , 6, in the field of view 200 in the object plane. The region Rp(x, y) is projected and focused to sensor element pxl(x, y) through the lens 201 of the imaging system 250. In this example, there will be 8×6 shading correction coefficients, each of which is an integrated (average) coefficient of the corresponding region.

[0011] If the image sensor of an imaging system has N_P by M_P elements, and IC_P(x, y) is the image pixel intensity of pixel pxl(x, y) of an uniform calibration standard, the pixel-by-pixel shading error SE_P(x, y) is defined in EQ. 1. Typically, shading errors are normalized so that the average of all errors equals to one. $\begin{array}{cc}{\mathrm{SE}}_P\left(x,y\right)=\frac{{\mathrm{IC}}_P\left(x,y\right)}{\frac{1}{N_P\times M_P}\sum _{i=1}^{N_P}\sum _{j=1}^{M_P}{\mathrm{IC}}_P\left(i,j\right)}\mathrm{For}x=1,\dots ,N_P\mathrm{and}y=1,\dots ,M_P& \mathrm{EQ}.1\end{array}$

[0012] Typically, an imaging system is a linear device. That is, the raw pixel gray levels of the imaging system have a linear relationship with the input signal. For example, when imaging an illuminant object, the input signal is the self-emitted light intensity from the object. In the linear case, image intensity IC_P(x, y) is the raw pixel gray-level of the pixel pxl(x, y) minus a stored background value of the corresponding pixel. In the non-linear case, the non-linear relationship first is modeled, and then this non-linear model is used in calculating the imaging intensity IC_P(x, y) from the raw gray-levels so that the image intensity IC_P(x, y) is linear with the input signal. If the image intensity is not linear to the input signal or the background values are not properly removed, the calculated shading error will not be accurate. In such a system, a skilled person will be capable of adjusting the image intensity values so that they have a linear relationship with the input signal.

[0013] FIG. 2 generalizes the pixel-by-pixel shading error concept. Instead of establishing shading correction coefficients of a region corresponding to each single pixel for the entire field of view, shading correction coefficients can be established for regions corresponding to more than one pixel. The pixel-by-pixel example of FIG. 1 is a special case of this generalized model.

[0014] In the example shown in FIG. 2, there are 3×2 regions 215, denoted as R(x, y): x=1, 2, 3 and y=1, 2. Each of these regions projects to its corresponding pixel set 216, denoted as P(x, y): x=1, 2, 3 and y=1, 2, in the imaging sensor 202. Each pixel set in this example contains four pixels. For example, P(1, 1), corresponding to R(1, 1), consists of: {pxl(1, 1), pxl(2, 1), pxl(1, 2), pxl(2, 2)}. The generalized shading correction has 3×2 coefficients, and each coefficient corresponds to an area of 4 pixels. Of course, the regions do not have to be in a regular 2-dimensional matrix form. They can be in any form. For the convenience of discussion, the regular 2-dimensional matrix form is used in this document. The calibration standard may be 1-dimensional or 3-dimensional or may have a non-rectangular 2-dimensional form.

[0015] Given a calibration standard with uniform regions R(x, y): x=1, . . . , N and y=1, . . . , M, the shading error SE(x, y) for the region R(x, y) (or pixel set P(x, y)) is defined in EQ. 2. $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{\mathrm{IC}\left(x,y\right)}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M\mathrm{IC}\left(i,j\right)}\mathrm{For}x=1,\dots ,N\mathrm{and}y=1,\dots ,M& \mathrm{EQ}.2\end{array}$
where: IC(x, y) is the image intensity of region R(x, y). It is calculated by averaging the image pixel intensities of all the pixels corresponding to region R(x, y).

[0016] Once the shading error is defined, the shading correction coefficient SC(x, y) for the region R(x, y) is defined in EQ. 3. $\begin{array}{cc}\mathrm{SC}\left(x,y\right)=\frac{1}{\mathrm{SE}\left(x,y\right)}\mathrm{For}x=1,\dots ,N\mathrm{and}y=1,\dots ,M& \mathrm{EQ}.3\end{array}$

[0017] The shading correction coefficients SC(x, y) will then be stored and applied to subsequent images to compensate for the system's shading error.

[0018] The key requirement of the prior art is uniformity. However, uniform calibration standards are often difficult or impossible to prepare, and, therefore, the resulting estimated shading correction coefficients are often inaccurate.

[0019] Another method for calculating the shading correction coefficients is using calibration standards with non-uniform but known local reflectance, transmittance or emission values. Again, it is often very difficult to accurately determine the required reflectance, transmittance or emission values at all the required locations, especially in microscopic imaging applications.

[0020] The present invention provides a method for calculating the shading correction coefficients for an imaging system from non-uniform and unknown calibration standards.

[0021] A calibration standard which corresponds to the item or material to be imaged is first selected. The calibration standard may have one, two or three dimensions. Potentially, the calibration standard may have more dimensions. Preferably, the calibration standard covers all or most of the field of view of the imaging system. The calibration standard is divided into regions and imaged to calculate an initial image intensity. The calibration standard is re-arranged and imaged again to determine a set of re-arranged image intensities. The process of re-arranging and image is repeated once for each dimension of the calibration standard.

[0022] The initial image intensity and re-arranged image intensities are combined mathematically to form a set of equations in which a shading error for each region may be calculated based on the shading error of another region or regions. By selecting a value for the shading error of one region, an intermediate shading error of each region may be calculated. The shading error may then be calculated as the normalized value of the intermediate shading error. The shading correction coefficient for each region may be calculated by inverting the shading error for that region.

[0023] The shading correction coefficients may then be used to correct the image intensity values for an actual object imaged using the imaging system.

[0024] The calibration standard need not have uniform regions, as long as each region has the same characteristics as the object to be imaged. Differences between the regions are abstracted out of the determination of the shading correction coefficient by the present invention.

[0025] In one aspect, the present invention provides a method for calculating shading correction coefficients for an imaging system, comprising the steps of defining a set of calibration regions on a calibration standard; defining a first part and a second part of said calibration standard, wherein each of said first and second parts contains at least one calibration region; taking an image of said calibration standard and calculating an initial image intensity for each of said calibration regions; re-arranging said first and second parts to form a re-arrangement of said calibration standard; taking an image of said re-arrangement and calculating a re-arranged image intensity for each of said calibration regions; calculating a shading error for each of said calibration regions; and calculating a shading correction coefficient for each of said calibration regions.

[0026] In another aspect, the present invention provides a method for calculating shading correction coefficients for an imaging system, comprising the steps of: defining a set of calibration regions on a calibration standard; defining a first part, a second part, a third part and a fourth part of said calibration standard, wherein each of said first, second, third and fourth parts contains at least one calibration region; taking an image of said calibration standard and calculating an initial image intensity for each of said calibration regions; re-arranging said first, second, third and fourth parts in a first direction to form a first re-arrangement of said calibration standard; taking an image of said first re-arrangement and calculating a first re-arranged image intensity for each of said calibration regions; re-arranging said first, second, third and fourth parts in a second direction to form a second re-arrangement of said calibration standard; taking an image of said second re-arrangement and calculating a second re-arranged image intensity for each of said calibration regions; calculating a shading error for each of said calibration regions; and calculating a shading correction coefficient for each of said calibration regions.

[0047] The present invention provides a method for calculating shading correction coefficients for an imaging system, using a calibration standard with N×M regions. Reference is first made to FIG. 3, which illustrates a first calibration standard 10 according to the present invention. The described calibration standard 10 consists of N×M (N columns and M rows) calibration regions 5, which are referred to as R(x, y). Calibration standard 10 is divided into four parts. The variables x and y are the array index. x takes values form 1 to N, denoted as x=1, . . . , N, and y takes values from 1 to M, denoted as y=1, . . . , M. The top-left corner 4 consists of one calibration region 5 at its center; the top-row part 2 consists of N−1 calibration regions 5; the left-column part 3 consists of M−1 calibration regions 5, and bottom-right part 1 consists of (N−1)×(M−1) calibration regions 5. In this exemplary embodiment, N>1 and M>1. This is referred to as the 2D (2-dimensional) case. In the 2D case, both N and M have to be greater than one. If either N or M is equal to one, it becomes a one dimensional case, which is described later later.

[0048] The location of part 4 (i.e., the part with one calibration region R(1, 1)) of calibration standard 10 is used to name the configuration of the calibration standard. The calibration standard shown in FIG. 3 has a top-left configuration, and the arrangement of the four parts shown in FIG. 3 is the initial setup of the top-left configuration.

[0049] Calibration standard 10 is imaged by an imaging system and a shading correction coefficient for each region 5 is calculated according to the present invention. The shading correction coefficients for regions 5 will depend on the object being imaged and the particular characteristic being imaged. For example, a fluorescent object may be imaged to determine the amount of excited light from the object. Subsequently, the same imaging system may be used to image an illuminant object to determine the amount of self-emitted light. In these two different cases, the shading correction coefficients may be different, even for comparable regions of the imaging field. In the first case, the shading error is due to the inhomogenity of the lens, the excitation light source, sensor and other electronic sources. In the second case, the shading error is due to the inhomogenity of the lens, sensor and other electronic sources (there is no excitation light source).

[0050] It is preferable that the material in each calibration region 5 is the same or has characteristics similar to the object which will ultimately be imaged. For example, the same fluorescent material is preferably deposited on the calibration regions 5 if the object to be imaged is fluorescent. However, unlike the prior art, it is not necessary that the densities or concentrations of the deposited material in regions 5 of calibration standard 10 are the same and the differences between regions do not need to be known.

[0051] The physical shape of the calibration regions 5 can be in any form. The square shape shown in the Figures is for illustration purpose only. Also, the gaps between calibration regions 5 are for the purpose of illustration only. Their physical presence is not required. Also, the central placement of calibration regions 5 in parts 1, 2, 3 and 4 are for the illustration purpose only. Their placement could be off the center as long as the region intervals remain uniform.

[0052] The signal intensity of region R(x, y) is denoted as VC(x, y). In the case of illuminant calibration standards, the signal intensity is the emitted light intensity. In the case of fluorescent calibration standards, the signal intensity is the excited light intensity under uniform illumination since the illumination inhomogenity is already incorporated into the shading error. The signal intensity VC(x, y) of a particular region R(x, y) should remain constant during the calibration process described below.

[0053] The dimensions of the calibration standard 10 are shown in FIG. 4, in which “a” is the horizontal interval of regions 5 and “b” is the vertical interval of regions 5. Preferably, all the calibration regions 5 should be within the field of view of the imaging system. Also, preferably, the spatial dimension of the calibration standard 10 should be equal to that of the field of view of the imaging system to be calibrated. That is:
W_FOV=a×N and H_FOV=b×M
where: W_FOV and H_FOV are the width and the height of the field of view of the imaging system to be calibrated.

[0054] The selections of a, b, N, M and the shape of the calibration regions 5 are determined by requirements of a particular application of the present invention. Where a large variation in shading errors is expected, a and b may be selected to be small and consequently N and M (assuming a fixed field of view) will be large. That is, the regions 5 are spaced more closely together.

[0055] When the calibration standard shown in FIG. 3 is imaged in an imaging system, the imaging intensity IC^O(x, y) of the region R(x, y) is defined by the following equation:
IC^O(x, y)=G×SE(x, y)×VC(x, y)  EQ. 4

• • [0056] for x=1, . . . , N and y=1, . . . , M
    where: G is a constant representing the imaging system gain.

[0057] Taking a natural log transform on both side of EQ. 4, we obtained:
IC_L^O(x, y)=G_L+SE_L(x, y)+VC_L(x, y)  EQ. 5

• • [0058] for x=1, . . . , N and y=1, . . . , M
    where:
    IC_L^O(x, y)=In (IC^O(x, y))
    SE_L(x, y)=In (SE(x, y))
    VC_L(x, y)=In (VC(x, y))
    G_L=In (G)
  • [0059] ln(d) is the natural log of number d.

[0060] FIG. 5 shows a horizontally shifted-and-rotated arrangement 11 of the top-left configuration of calibration standard 10 (FIGS. 3 and 4). The horizontal shift-and-rotation is done by shifting the entire calibration standard 10 to the left by one column and then moving parts 3 and 4 to the right of parts 1 and 2. The re-arranged calibration regions 5 are referred to as R^H(x, y), and the relationship between the contents (calibration materials) of R^H(x, y) and R(x, y) is defined as:
Content(R^H(x, y))=Content(R(Mod(x+1, N), y))  EQ. 6

• • [0061] for x=1, . . ., N and y=1, . . ., M

[0062] The function Mod (q, d) in EQ. 6 and subsequent equations is defined as: $\mathrm{Mod}\left(q,d\right)=\{\begin{array}{cc}\left(q+m\times d-1\right)\%d+1& \mathrm{if}q\mathrm{and}d\mathrm{are}\mathrm{integers}\mathrm{and}q<1\mathrm{and}d\ge 1\\ \left(q-1\right)\%d+1& \mathrm{if}q\mathrm{and}d\mathrm{are}\mathrm{integers}\mathrm{and}q\ge 1\mathrm{and}d\ge 1\\ \mathrm{Invalid}& \mathrm{Otherwise}\end{array}$
where: % is a remainder operator (e.g., 1%3=10%3=1).

[0063] Given integer q and d, m is the smallest positive integer such that q+m×d≧1.

[0064] FIG. 6 shows a vertically shifted-and-rotated arrangement 12 of the top-left configuration of calibration standard 10 (FIGS. 3 and 4). The vertical shift-and-rotation is done by shifting the entire calibration standard 10 up by one row and then moving parts 2 and 4 to the bottom of parts 1 and 3. The re-arranged calibration regions 5 are referred to as R^V(x, y), and the relationship between the contents of R^V(x, y) and R(x, y) is defined as:
Content(R^V(x, y)) Content(R(x, Mod(y+1, M)))  EQ. 7

• • [0065] for x=1, . . . , N and y=1, . . . , M

[0066] According to the present invention, the re-arranged calibration regions 5 of the shifted-and-rotated arrangements 11 and 12 are physically aligned to the calibration regions 5 in the initial setup of calibration standard 10. That is:
P(x, y)=P^H(x, y)=P^V(x, y)  RQ. 1

• • [0067] for x=1, . . . , N and y=1, . . . , M
    where: P(x, y), P^H(x, y), and P^V(x, y) are the mapped pixel sets of calibration regions R(x, y) (FIG. 3), R^H(x, y) (FIG. 4) and R^V(x, y) (FIG. 5) respectively.

[0068] Let IC_L^H(x, y) be the natural log transformed image intensity of the region R^H(x, y) for the horizontally shifted-and-rotated arrangement 11 shown in FIG. 5. From RQ. 1, we know that the region R(x, y) and the region R^H(x, y) map to the same pixels. Therefore, they share the same shading error SE(x, y). From EQ. 6, we further know that the content of region R^H(x, y) is the content of region R(Mod(x+1, N), y), i.e., the signal intensity of calibration region R^H(x, y) is VC(Mod(x+1, N), y). Therefore, we have:
IC_L^H(x, y)=G_L+SE_L(x, y)+VC_L(Mod(x+1, N), y)  EQ. 8

• • [0069] for x=1, . . . , N and y=1, . . . , M

[0070] Let IC_L^V(x, y) be the natural log transformed image intensity of the region R^V(x, y) for the vertically shifted-and-rotated arrangement 12 shown in FIG. 6. Similar to the deduction of EQ. 8 from EQ. 7 and RQ. 1, we have:
IC_L^V(x, y)=G_L+SE_L(x, y)+VC_L(x, Mod(y+1, M))   EQ. 9

• • [0071] for x=1, . . . , N and y=1, . . . , M

[0072] Multiplying 2 to both sides of EQ. 5, multiplying−1 to both sides of EQ. 8 and EQ. 9, and then adding the three new equations together, i.e., 2×(EQ. 5)−(EQ. 8)−(EQ. 9), we have:
ICV(x, y)=2×VC_L(x, y)−VC_L(Mod(x+1, N), y)−VC_L(x, Mod(y+1, M))  EQ. 10

• • [0073] for x=1, . . . , N and y=1, . . . , M
    where:
    ICV(x, y)=2×IC_L^O(x, y)−IC_L^H(x, y)−IC_L^V(x, y)

[0074] Using EQ. 5, we can obtain the following relationship:
VC_L(x, y)=IC_L^O(x, y)−G_L−SE_L(x, y)
VC_L(Mod(x+1, N), y)=IC_L^O(Mod(x+1, N), y)−G_L−SE_L(Mod(x+1, N), y)
VC_L(x, Mod(y+1, M))=IC_L(x, Mod(y+1, M))−G_L−SE_L(x, Mod(y+1, M))

[0075] Substituting VC_L(x, y), VC_L(Mod(x+1, N), y) and VC_L(x, Mod(y+1, M)) in EQ. 10 with the above relationship provides EQ. 11:
ICS(x, y)=2×SE_L(x, y)−SE_L(Mod(x+1, N), y)−SE_L(x, Mod(y+1, M))  EQ. 11

• • [0076] for x=1, . . . , N and y=1, . . . , M
    where:
    ICS(x, y)=IC_L^H(x, y)+IC_L^V(x, y)−IC_L^O(Mod(x+1, N), y)−IC_L^O(x, Mod(y++1, M))

[0077] If we order the two dimensional index (x, y) in EQ. 11 by the column first method, i.e., (1, 1), (2, 1), . . . , (N, 1), (1, 2), (2, 2), . . . , (N, 2), . . . , (1, M), (2, M), . . ., (N, M), we can represent EQ. 11 in a matrix representation, as defined in EQ. 12.
ICS=A SE_L  EQ. 12
where: A is a (N×M)×(N×M) matrix, i.e., it has N×M rows and N×M columns. It can be represented by a partitioned matrix, as defined below. The matrix E is a N×N matrix as defined below, the matrix I is a N×N identity matrix and matrix 0 is a N×N zero matrix (all coefficients are 0). $A={\left[\begin{array}{ccccccc}E& -I& 0& .& 0& 0& 0\\ 0& E& -I& .& 0& 0& 0\\ 0& 0& E& .& 0& 0& 0\\ .& .& .& .& .& .& .\\ 0& 0& 0& .& E& -I& 0\\ 0& 0& 0& .& 0& E& -I\\ -I& 0& 0& .& 0& 0& E\end{array}\right]}_{\mathrm{NxM}x\mathrm{NxM}}$ $E={\left[\begin{array}{ccccc}2& -1& .& 0& 0\\ 0& 2& .& 0& 0\\ .& .& .& .& .\\ 0& 0& .& 2& -1\\ -1& 0& .& 0& 2\end{array}\right]}_{\mathrm{NxN}}$ $I={\left[\begin{array}{ccccc}1& 0& .& 0& 0\\ 0& 1& .& 0& 0\\ .& .& .& .& .\\ 0& 0& .& 1& 0\\ 0& 0& .& 0& 1\end{array}\right]}_{\mathrm{NxN}}$ $\mathrm{ICS}={\lfloor \begin{array}{c}\mathrm{ICS}\left(1,1\right)\\ \mathrm{ICS}\left(2,1\right)\\ .\\ \mathrm{ICS}\left(N,1\right)\\ .\\ .\\ \mathrm{ICS}\left(1,M\right)\\ \mathrm{ICS}\left(2,M\right)\\ .\\ \mathrm{ICS}\left(N,M\right)\end{array}\rfloor }_{1\times \left(N\times M\right)}$ ${\mathrm{SE}}_L={\lfloor \begin{array}{c}{\mathrm{SE}}_L\left(1,1\right)\\ {\mathrm{SE}}_L\left(2,1\right)\\ .\\ {\mathrm{SE}}_L\left(N,1\right)\\ .\\ .\\ {\mathrm{SE}}_L\left(1,M\right)\\ {\mathrm{SE}}_L\left(2,M\right)\\ .\\ {\mathrm{SE}}_L\left(N,M\right)\end{array}\rfloor }_{1\times \left(N\times M\right)}$

[0078] Using standard mathematical techniques, it can be proved that the rank of matrix A is N×M−1. Therefore the solution of SE_L(x, y) is not unique. However, since the rank of matrix A is N×M−1, SE_L(x, y) can be resolved by arbitrarily specifying only one of the values (e.g., SE_L(1, 1) is specified).

[0079] Let SEO be a solution of EQ. 11 or EQ. 12, by giving one of the shading errors a known value. SEO can be defined as
{SEO(x, y, ax, ay, K): x=1, . . . , N, y=1, . . . , M, 1≦ax≦N, 1≦ay≦M, K is a given value and SE_L(ax, ay)=K}  DF. 1
Taking the exponential transform (inverse of the log transform) and normalization (as in EQ. 2), the corresponding shading error SE(x, y) is defined as: $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{e^{\mathrm{SEO}\left(x,y,\mathrm{ax},\mathrm{ay},K\right)}}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{e}^{\mathrm{SEO}\left(i,j,\mathrm{ax},\mathrm{ay},K\right)}}\mathrm{for}x=1,\dots ,N\mathrm{and}y=1,\dots ,M& \mathrm{EQ}.13\end{array}$

[0080] Using standard mathematical techniques, it can be proved that the difference between any given two solutions SEO(x, y, ax1, ay1, K1) and SEO(x, y, ax2, ay2, K2) of EQ. 12 is constant, for a given set of measured image intensities (IC_L^O(x, y), IC_L^H(x, y) and IC_L^V(x, y)). That is:
SEO(x, y, ax2, ay2, K2)−SEO(x, y, ax1, ay1, K1)=F

• • [0081] for x=1, . . . , N and y=1, . . . , M
    where: F is a constant.

[0082] This give the following equation: $\begin{array}{c}\frac{e^{\mathrm{SEO}\left(x,y,\mathrm{ax2},\mathrm{ay2},\mathrm{K2}\right)}}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{e}^{\mathrm{SEO}\left(i,j,\mathrm{ax2},\mathrm{ay2},\mathrm{K2}\right)}}=\frac{e^{\mathrm{SEO}\left(x,y,\mathrm{ax1},\mathrm{ay1},\mathrm{K1}\right)+F}}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{e}^{\mathrm{SEO}\left(i,j,\mathrm{ax1},\mathrm{ay1},\mathrm{K1}\right)+F}}\\ =\frac{e^{\mathrm{SEO}\left(x,y,\mathrm{ax1},\mathrm{ay1},\mathrm{K1}\right)}}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{e}^{\mathrm{SEO}\left(i,j,\mathrm{ax1},\mathrm{ay1},\mathrm{K1}\right)}}\end{array}$

[0083] That is, the shading error SE(x, y) defined in EQ. 13 is unique, regardless of the values of ax, ay and K.

[0084] As an example, select ax=N, ay=M and K=1 as the solution for calculating the shading error. This solution is referred to as intermediate shading error and denoted as S1(x, y), i.e., S1(x, y)=SEO(x, y, N, M, 1) , whose solution is defined by EQ. 14A and EQ. 14B. $\begin{array}{cc}\mathrm{S1}\left(N,M\right)=1& \text{EQ. 14A}\\ {\lfloor \begin{array}{c}\mathrm{ICS}\left(1,1\right)\\ \mathrm{ICS}\left(2,1\right)\\ .\\ \mathrm{ICS}\left(N,1\right)\\ .\\ .\\ \mathrm{ICS}\left(N,M-1\right)+1\\ \mathrm{ICS}\left(1,M\right)\\ \mathrm{ICS}\left(2,M\right)\\ .\\ \mathrm{ICS}\left(N-1,M\right)+1\end{array}\rfloor }_{1\times \left(N\times M-1\right)}=\mathrm{A1}{\lfloor \begin{array}{c}\mathrm{S1}\left(1,1\right)\\ \mathrm{S1}\left(2,1\right)\\ .\\ \mathrm{S1}\left(N,1\right)\\ .\\ .\\ \mathrm{S1}\left(1,M\right)\\ \mathrm{S1}\left(2,M\right)\\ .\\ \mathrm{S1}\left(N-1,M\right)\end{array}\rfloor }_{1\times \left(N\times M-1\right)}& \text{EQ. 14B}\end{array}$
where: Matrix A1 is the top-left (N×M−1)×(N×M−1) sub-matrix of A defined in EQ. 12. That is:
a1_i,j=a_i,j for i=1, . . . , N−1 and j=1, . . . , M−1.

• • [0085] a1_i,j and a_i,j are the coefficient of matrix A1 and A at column i and row j respectively.
    The constant “1” in the terms “ICS(N, M−1)+1” and “ICS(N, M−1)+1” results from moving the variable SE_L(N, M) to the left side of the equation and SE_L(N, M)=1.

[0086] EQ. 14B can be solved by directly calculating the inverse matrix A1⁻¹, as defined in EQ. 15. It can also be solved by standard iterative methods. $\begin{array}{cc}{\lfloor \begin{array}{c}\mathrm{S1}\left(1,1\right)\\ \mathrm{S1}\left(2,1\right)\\ .\\ \mathrm{S1}\left(N,1\right)\\ .\\ .\\ \mathrm{S1}\left(1,M\right)\\ \mathrm{S1}\left(2,M\right)\\ .\\ \mathrm{S1}\left(N-1,M\right)\end{array}\rfloor }_{1\times \left(N\times M-1\right)}={\mathrm{A1}}^{-1}{\lfloor \begin{array}{c}\mathrm{ICS}\left(1,1\right)\\ \mathrm{ICS}\left(2,1\right)\\ .\\ \mathrm{ICS}\left(N,1\right)\\ .\\ .\\ \mathrm{ICS}\left(N,M-1\right)+1\\ \mathrm{ICS}\left(1,M\right)\\ \mathrm{ICS}\left(2,M\right)\\ .\\ \mathrm{ICS}\left(N-1,M\right)+1\end{array}\rfloor }_{1\times \left(N\times M-1\right)}& \mathrm{EQ}.15\end{array}$

[0087] The shading error SE(x, y) can then be calculated from the intermediate shading error by EQ. 16: $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{e^{\mathrm{S1}\left(x,y\right)}}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{e}^{\mathrm{S1}\left(i,j\right)}}\mathrm{for}x=1,\dots ,N\mathrm{and}y=1,\dots ,M& \mathrm{EQ}.16\end{array}$

[0088] Once shading error SE(x, y) is determined, the shading correction coefficients SC(x, y) can be determined from EQ. 3.

[0089] FIG. 7 summarizes the shading correction coefficient calculation. A calibration standard 10 with a top-left configuration is present to an imaging system to be calibrated. Step 101 utilizes the calibration standard 10 in its initial setup as illustrated in FIG. 3. Step 102 manipulates calibration standard 10 to form a horizontally shifted-and-rotated arrangement 11 (FIG. 5). Step 103 manipulates the standard to form the vertically shifted-and-rotated arrangement 12 (FIG. 6). In step 104, images are taken of initial setup of calibration standard 10, shifted-and-rotated arrangement 11 and shifted-and-rotated arrangement 12, and the pixel sets P(x, y), P^H(x, y) and P^V(x, y) of the corresponding calibration regions R(x, y), R^H(x, y) and R^V(x, y) are determined using standard image processing algorithms. Step 105 calculates the image intensity IC^O(x, y), IC^H(x, y) and IC^V(x, y) from the pixel sets P(x, y), P^H(x, y) and P^V(x, y) and the acquired images. Step 106 calculates the intermediate shading error S1(x, y) from EQ. 14A and EQ. 14B. Step 107 calculates the shading error SE(x, y) from the intermediate shading error S1(x, y), by EQ. 16. Finally, step 108 calculates the shading correction coefficient SC(x, y) from SE(x, y) by EQ. 3.

[0090] The region signal intensity VC(x, y) can also be determined from EQ. 4, if the gain of the imaging system is known. That is: $\begin{array}{cc}\mathrm{VC}\left(x,y\right)=\frac{{\mathrm{IC}}^0\left(x,y\right)}{G\times \mathrm{SE}\left(x,y\right)}& \mathrm{EQ}.17\end{array}$

[0091] If only the relative region signal intensity (VC_R(x, y)) is required, it can be calculated without knowledge of the system gain. That is: $\begin{array}{cc}{\mathrm{VC}}_R\left(x,y\right)=\frac{{\mathrm{IC}}^0\left(x,y\right)/\mathrm{SE}\left(x,y\right)}{\frac{1}{N\times M}\sum _{i=1}^N\sum _{j=1}^M{\mathrm{IC}}^0\left(i,j\right)/\mathrm{SE}\left(i,j\right)}& \mathrm{EQ}.18\end{array}$

[0092] FIG. 3 shows calibration standard 10 which has a top-left configuration. There are many other configurations that can be used to calculate shading errors and the region signal intensities. For example, FIG. 8 shows a calibration standard 60 with a bottom-right configuration. As with the top-left calibration standard 10 described in relation to FIG. 3, this bottom-right calibration standard 60 has four detachable parts with N×M uniformly arranged calibration regions 5. The bottom-right corner 24 consists of one calibration region 5 at the center; the bottom-row part 22 consists of N−1 calibration regions 5; the right-column part 23 consists of M−1 calibration regions 5, and the part 21 consists of (N−1)×(M−1) calibration regions 5. The calibration regions 5 are referred to as R^BR(x, y).

[0093] FIG. 9 shows a horizontally shifted-and-rotated arrangement 61 of the bottom-right calibration standard 60. The shift-and-rotation is done by shifting the entire calibration standard 60 (in FIG. 8) to right by one column and then moving the part 23 and 24 to the left of part 21 and 22. The re-arranged calibration regions 5 are referred to as P^HR(x, y), and:
Content(R^HR(x, y))=Content (R^BR(Mod(x−1, N), y))  EQ. 19

• • [0094] for x=1, . . . , N and y=1, . . . , M

[0095] FIG. 10 shows a vertically shifted-and-rotated arrangement 62 of the bottom-right calibration standard 60. The shift-and-rotation is done by shifting the entire calibration standard 60 (in FIG. 8) down by one row and then moving the part 22 and 24 to the top of part 21 and 23. The re-arranged calibration regions 5 are referred to as P^VD(x, y), and:
Content(R^VD(x, y))=Content(R^BR(x, Mod(y−1, M)))  EQ. 20

• • [0096] for x=1, . . . , N and y=1, . . . , M

[0097] P^BR(x, y), P^HR(x, y) and P^VD(x, y) are the pixel sets corresponding to regions R^BR(x, y), R^HR(x, y) and R^VD(x, y) respectively. Similar to the top-left configuration, it is required that the re-arranged calibration regions of the horizontally and vertically shifted-and-rotated arrangements 61 and 62 are aligned to the regions in the initial calibration standard 60. That is:
P^BR(x, y)=P^HR(x, y)=P^VD(x, y)  RQ. 2

• • [0098] for x=1, . . . , N and y=1, . . . , M

[0099] Similar to the deduction of EQ. 11, the equation for the shading error SE_L(x, y) of the bottom-right calibration standard 60 is:
ICS_BR(x, y)=2×SE_L(x, y)−SE_L(Mod(x−1, N), y)−SE_L(x, Mod(y−1, M))  EQ. 21

• • [0100] for x=1, . . . , N and y=1, . . . , M
    where:
    ICS_BR(x, y)=IC_L^HR(x, y)+IC_L^VD(x, y)−IC_L^BR(Mod(x−1, N), y)−IC_L^BR(x, Mod(y−1, M))
  • [0101] IC_L^BR(x, y) is the log transformed image intensity of the region
  • [0102] R^BR(x, y) shown in FIG. 8;
  • [0103] IC_L^HR(x, y) is the log transformed image intensity of the region R^HR(x, y) shown in FIG. 9; and
  • [0104] IC_L^VD(x, y) is the log transformed image intensity of the region R^VD(x, y) shown in FIG. 10.

[0105] Similar to the deduction of EQ. 12 (the matrix form of EQ. 11), the matrix form of EQ. 21 can be easily deduced. EQ. 21 has the same properties as those of EQ. 11, and the procedure for solving EQ. 21 is exactly the same as that for solving EQ. 11, except the coefficients of the matrix corresponding to EQ. 21 is different than those in EQ. 12.

[0106] In the previous discussion, the solution of shading error SE(x, y) from a top-left or a bottom-right calibration standard is derived from the image intensities of the initial calibration standard (10 or 60), a horizontally shifted-and-rotated arrangement (11 or 61) and the vertically shifted-and-rotated arrangement (12 or 62). There are other configurations and arrangements, which can be used to derive the shading error SE(x, y) for an imaging system.

[0107] Reference is next made to FIG. 11. In the previous discussion, the calibration standards (10 and 60) presented are rectangular. Calibration standards may be in other forms as well. For example, calibration standard 310 is a parallelogram. Calibration standard 310 has four parts 301, 302, 303 and 304. Each part has a number of calibration regions 305. In fact, a calibration standard can be of any form as long as the calibration regions of shifted-and-rotated arrangements can be aligned with the regions of the initial calibration standard.

[0108] The corner part (e.g., part 4 in FIG. 3) of calibration standard 10 (FIG. 3) consists of only one calibration region 5. Referring to FIG. 12, a calibration standard may alternatively have more than one region in the corner part. FIG. 12 shows a calibration standard 30 with a top-left configuration and with 2 calibration regions in the corner part 34. The top-row part 32 consists of N-2 calibration regions; the left-column part 33 consists of 2×(M−1) calibration regions, and bottom-right part 31 consists of (N-2)×(M−1) calibration regions. The horizontally shifted-and-rotated arrangement corresponding to calibration standard 30 will require shifting TWO columns. Similar to the deduction of EQ. 11, we can also derive an equation to solve the shading error SE(x, y) from the image intensities of the initial calibration standard, a horizontally shifted-and-rotated arrangement, and a vertically shifted-and-rotated arrangement. However, unlike EQ. 11, solving this equation requires specifying TWO values of SE(x, y), e.g., SE(1, 1) and SE(2, 1) are given. Therefore, more than one initial value of the shading error SE(x, y) needs to be defined if the corner parts of a calibration standard consists of more than one calibration region. In order to obtain the correct shading correction coefficients for an imaging system using a calibration standard 30, the ratio of the shading errors for the two specified regions must be known a priori.

[0109] In certain applications, fabrication of a physically detachable calibration standard may not be convenient or possible. In this case, the shift-and-rotation can be achieved by performing two shifts of a calibration standard. FIG. 13 shows calibration standard 70 with N×M calibration regions RS(x, y) with an imaginary top-left configuration. The dotted lines 81 and 82 in FIG. 13 represent the imaginary divisions to the standard 70 into a top-left configuration, and imaginary parts 71, 72, 73 and 74 correspond to part 1, 2, 3 and 4 in FIG. 3. It is to be noted that the use of the word “imaginary” here indicates only that calibration standard 70 cannot be conveniently detached to separate parts 71, 72, 73 and 74.

[0110] FIG. 14 shows the process to obtain the image of a horizontally shifted-and-rotated arrangement 42 of calibration standard 70 (FIG. 13).

[0111] First, calibration standard 70 is shifted to the left by 1 column so that line 81 is aligned with the left edge of the image of calibration standard 70 in its original position. Image 40 is acquired with calibration standard 70 in this position. The shaded area of image 40 contains regions in the imaginary parts 71 and 72 of calibration standard 70.

[0112] Second, calibration standard 70 is shifted to the right by N−1 column (relative to its original position) so that line 81 is aligned with the right edge of the image of calibration standard 70 in its original position. Image 41 is acquired with calibration standard 70 in this position. The shaded area of image 41 contains regions in the imaginary parts 73 and 74 of calibration standard 70.

[0113] Finally, images 40 and 41 are combined to form horizontally shifted-and-rotated arrangement 42. Areas 71 and 72 in horizontally shifted-and-rotated arrangement 42 are copied from the corresponding areas in image 40, and areas 73 and 74 in horizontally shifted-and-rotated arrangement 42 are copied from the corresponding areas in image 41. The combined horizontally shifted-and-rotated arrangement 42 is essentially the same as if calibration standard 70 had been detached along line 81 and then re-assembled in a horizontally shifted-and-rotated arrangement and then imaged.

[0114] FIG. 15 shows the process to obtain the image of a vertically shifted-and-rotated arrangement 52 of calibration standard 70.

[0115] First, calibration standard 70 is shifted up by 1 row so that line 82 is aligned with the top edge of the image of calibration standard 70 in its original position. Image 50 is acquired with calibration standard 70 in this position. The shaded area of image 50 contains regions in the imaginary parts 71 and 73 of calibration standard 70.

[0116] Secondly, calibration standard 70 down by M−1 rows (relative to its original position) so that line 82 is aligned with the bottom edge of the image of calibration standard 70 in its original position. Image 51 is then acquired with calibration standard in this position. The shaded area of image 51 contains regions in the imaginary parts 72 and 74 of calibration standard 70.

[0117] Finally, the image 50 and 51 are combined to form vertically shifted-and-rotated arrangement 52. Areas 71 and 73 in vertically shifted-and-rotated arrangement 52 are copied from the corresponding areas in image 50 and areas 72 and 74 in vertically shifted-and-rotated arrangement 52 are copied from the corresponding areas in image 51. The combined shifted-and-rotated arrangement 52 is essentially the same as if calibration standard 70 had been detached along line 82 and then re-assembled in a shifted-and-rotated calibration array and then imaged.

[0118] The shift operations in FIG. 14 and FIG. 15 can be accomplished by using either manual or computer-controlled XY positioning devices. It is preferable to have an accurate and reliable positioning device.

[0119] It is desirable that the regions RS(x, y) of horizontally and vertically shifted-and-rotated arrangements 42 and 52 align precisely with the regions RS(x, y) of the initial calibration standard. Due to various mechanical errors, this precise alignment may not be possible, especially in a microscopy application. Using a calibration standard with a top-left configurations as an example, a mis-alignment is indicated by: P(x, y)≠P^H(x, y) or P(x, y)≠P^V(x, y) for some or all (x, y). For example, P(1, 1)={pxl(1, 1), pxl(2, 1), pxl(3, 1)} and p^H(1, 1)={pxl(2, 1), pxl(3, 1), pxl(4, 1)}. If P(x, y) is used as reference, i.e., the calibration regions of the initial setup of the calibration standard, the mis-alignment of P^H(x, y) and/or P^V(x, y) with P(x, y) introduces measurement error in IC^H(x, y) and/or IC^V(x, y). If the mis-alignment is significant, local correction functions (one for each calibration region) can be established and used to correct the measurement error in IC^H(x, y) and/or IC^V(x, y) due to the mis-alignment. Given image intensity IC(x, y) of a region with corresponding pixel set P′(x, y), this local correction function estimates the imaging intensity IC′(x, y) of the same region at pixel set P(x, y). This location correction function can be established by shifting and imaging the initial setup of a calibration standard around its default location. The image intensities of the calibration regions at different locations can then be used to define the local correction function parameters, using standard mathematical procedures. The established correction function should then be incorporated in the procedure 105 as shown in FIG. 7 to correct the measurement error if mis-alignment happens.

[0120] The measured image intensities, such as IC^O(x, y), IC^H(x, y) and IC^V(x, y) in a top-left calibration standard (i.e. calibration standard 10 (FIG. 3)), are not error free. In additional to the mechanical positioning error mentioned above, the other sources of error are electronic, quantum noises and etc. These measurement errors may cause errors in calculation of shading errors, shading correction coefficients and region signal intensities.

[0121] For a top-left calibration standard, {overscore (IC)}^O(x, y), {overscore (IC)}^H(x, y) and {overscore (IC)}^V(x, y) are the true image intensities (without measurement errors). The actual measurements can be expressed as:
IC^O(x, y)={overscore (IC)}^O(x, y)×(1+e_P^O(x, y))+e_A^O(x, y)
IC^H(x, y)={overscore (IC)}^H(x, y)×(1+e_P^H(x, y))+e_A^H(x, y)
IC^V(x, y)={overscore (IC)}^V(x, y)×(1+e_P^V(x, y))+e_A^V(x, y)  EQ. 22

• • [0122] for x=1, . . . , N and y=1, . . . , M
    where: e_P^O(x, y), ep_P^H(x, y) and e_P^V(x, y) are the proportional measurement errors, and

[0123] e_A^O(x, y), e_A^H(x, y) and e_p^V(x, y) are the additive measurement errors.

[0124] The corresponding solutions of shading error can be expressed as:
SE(x, y)={overscore (SE)}(x, y)×(1+e_SE(x, y))  EQ. 23

• • [0125] for x=1, . . . , N and y=1, . . . , M
    where: {overscore (SE)}(x, y) is the true shading error without error, and
  • [0126] e_SE(x, y) is the calculation error due to the measurement error
  • [0127] e_P^O(x, y), e_P^H(x, y), e_P^V(x, y), e_A^O(x, y), e_A^H(x, y) and e_A^V(x, y).

[0128] The calculation error e_SE(x, y) is summarized as a proportional error because of the relative nature of the shading error. Error e_SE(x, y) depends on the measurement errors and size of the calibration standard (i.e., N and M). The larger the measurement errors, the larger the calculation error e_SE(x, y) , and the larger the size of the calibration array, the larger the calculation error e_SE(x, y) .

[0129] Similarly, we can define the calculation errors for other calibration standards (bottom right, etc.).

[0130] The magnitude of errors can be reduced by averaging multiple acquisitions and calculations from the same calibration standard, by averaging the results of different calibration standards, or by averaging the combination of both.

[0131] As an example of averaging multiple calculation from the same calibration standard, SE⁽¹⁾(x, y) and SE⁽²⁾(x, y) may be defined as the results of two acquisitions and calculations of shading errors of an imaging system from a calibration standard with a top-left configuration. The averaged shading error can be defined as: $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{{\mathrm{SE}}^{\left(1\right)}\left(x,y\right)+{\mathrm{SE}}^{\left(2\right)}\left(x,y\right)}{2}& \mathrm{EQ}.24\end{array}$

[0132] As an example of averaging different calibration standards, SE^TLL(x, y) may be defined as the results from a top-left calibration standard and SE^BR(x, y) may be defined as the results from a bottom-right calibration standard. A simple averaged shading error can be defined as: $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{{\mathrm{SE}}^{\mathrm{TL}}\left(x,y\right)+{\mathrm{SE}}^{\mathrm{BR}}\left(x,y\right)}{2}& \mathrm{EQ}.25\end{array}$

[0133] The statistics (such as mean and variance) of calculation error e_SE(x, y) is position-dependent, even if such statistics of measurement error is position-independent (e.g., all e_P^O(x, y), e_P^H(x, y), e_P^V(x, y), e_A^O(x, y), e_a^H(x, y) and e_A^V(x, y) are μ=0 and σ=0.1). By using standard statistical methods and knowledge of measurement error, the error mean μ_ESE(x, y) and error standard deviation σ_ESE(x, y) of e_SE(x, y) can be estimated for given measurement errors (e.g., Gaussian noise with μ=0 and σ=1 for the additive measurement errors and Gaussian noise with μ=0 and σ=0.05 for the proportional measurement errors). Using the error statistics μ_ESE(x, y) and σ_ESE(x, y), weighting factors can be derived. Instead of using a simple averaging, a weighted averaging can be performed. $\begin{array}{cc}\mathrm{SE}\left(x,y\right)=\frac{w^{\mathrm{TL}}\left(x,y\right)\times {\mathrm{SE}}^{\mathrm{TL}}\left(x,y\right)+w^{\mathrm{BR}}\left(x,y\right)\times {\mathrm{SE}}^{\mathrm{BR}}\left(x,y\right)}{w^{\mathrm{TL}}\left(x,y\right)+w^{\mathrm{BR}}\left(x,y\right)}& \mathrm{EQ}.26\end{array}$
where: w^TL(x, y) and w^BR(x, y) are the weighting factors for the top-left configuration and bottom-right configuration respectively.

[0134] The weighted average will reduce the overall error and will also reduce the variance of error statistics. The weighting factor can exist in many forms. For example: $\begin{array}{cc}{w}^{\mathrm{TL}}\left(x,y\right)=\frac{1}{{\sigma }_{\mathrm{ESE}}^{\mathrm{TL}}\left(x,y\right)}\mathrm{and}{w}^{\mathrm{BR}}\left(x,y\right)=\frac{1}{{\sigma }_{\mathrm{ESE}}^{\mathrm{BR}}\left(x,y\right)}& \mathrm{EQ}.27\end{array}$

[0135] Calculation error reduction of shading error will reduce the calculation error of shading correction coefficients. EQ. 23 to EQ. 27 describes error analysis and averaging methods for shading error. Similarly, we could also directly analyze and average shading correction coefficients as well.

[0136] In calibration standard N and M were defined to be greater than 1. If one of N or M is equal to one, the system becomes a 1D (1-dimensional) case. For example, the image sensor array consists of one row or one column of sensor elements, or only one row of regions is of interest. Without loss of generality, the following discussions will be based on a horizontal 1D array. A vertical 1D array will be simply rotated 90 degree. The naming convention for the calibration standards and shifted-and-rotated arrangements is similar to the 2D case.

[0137] FIG. 16 illustrates a one dimensional calibration standard 100 with a left configuration. The calibration standard 100 contains N uniformly spaced regions. In the 1D case, N has to be greater than one. The left part 102 consists of one region 5, and right part consists of N−1 regions 5. FIG. 17 shows the spatial dimensions of standard 100. “a” is the interval of regions 5. The regions 5 are referred to as R1(x).

[0138] FIG. 18 illustrates a shifted-and-rotated arrangement 110 of calibration standard 100. The shift-and-rotation is done by shifting the entire calibration standard 100 (in FIG. 16) to right by one column and then moving part 102 to the right of part 101. The re-arranged calibration regions 5 are referred to as R1^SR(x), and:
Content(R1^SR(x))=Content(R1(Mod(x+1, N))) for x=1, . . . , N  EQ. 28

[0139] As with the 2D case, calibration regions of the shifted-and-rotated arrangement 110 are required to be in alignment with the regions 5 of the initial calibration standard 100. That is:
P1^SR(x)=P1(x) for x=1, . . . , N  RQ. 3
where: P1(x) is the pixel set corresponding to region R1(x) , and

[0140] P1^SR(x) is the pixel set corresponding to region R1^SR(x).

Similar to the 2D case, the following relationship can be developed:
IC1^O(x)=G1×SE1(x)×VC1(x) for x=1, . . . , N  EQ. 29
where: IC1^O(x) is the image intensity of the region R1(x, y) shown in FIG. 16.

[0141] SE1(x) is the shading error of the region R1(x, y).

[0142] VC1(x) is the signal intensity of the region R1(x, y)

[0143] G1 is the imaging system gain.

And:
IC1_L^O(x)=G1_L+SE1_L(x)+VC1_L(x) for x=1, . . . , N  EQ. 30
IC1_L^SR(x)=G1_L+SE1_L(x)+VC1_L(Mod(x+1, N)) for x=1, . . . , N  EQ. 31
where:
IC1_L^O(x)=ln(IC1^O(x))

• • [0144] IC1_L^SR(x)=ln(IC^SR(x)), and IC1^SR(x) is the image intensity of region 5
  • [0145] R1^SR(x, y) shown in FIG. 18.
    SE1_L(x)=ln(SE1_L(x))
    VC1_L(x)=ln(VC1(x))
    G1_L=ln(G1)
    And:
    ICS1(x)=SE1_L(x)−SE1_L(Mod(x+1, N)) for x=1, . . . , N  EQ. 32
    where:
    ICS1(x)=IC1_L^SR(x)−IC1_L^O(Mod(x+1, N))

[0146] As with the 2D case, SE1_L(X) can be resolved from EQ. 32 if one of the values is specified.

[0147] Similar to the 2D case, SE1(x) is defined as the solution (intermediate shading error) of EQ. 32 when SE1_L(N)=1, and the shading errors, defined in EQ. 33, are uniquely defined. $\begin{array}{cc}\mathrm{SE1}\left(x\right)=\frac{e^{\mathrm{S11}\left(x\right)}}{\frac{1}{N}\sum _{i=1}^N{e}^{\mathrm{S11}\left(i\right)}}\mathrm{for}x=1,\dots ,N& \mathrm{EQ}.33\end{array}$

[0148] Once the shading error is calculated, the shading correction coefficients SC1(x), signal intensity VC1(x), and relative signal intensity VC1_R(x) can be calculated. $\begin{array}{cc}\mathrm{SC1}\left(x\right)=\frac{1}{\mathrm{SE1}\left(x\right)}\mathrm{for}x=1,\dots ,N& \mathrm{EQ}.34\\ \mathrm{VC1}\left(x\right)=\frac{{\mathrm{IC1}}^0\left(x\right)}{\mathrm{G1}\times \mathrm{SE1}\left(x\right)}\mathrm{for}x=1,\dots ,N& \mathrm{EQ}.35\\ {\mathrm{VC1}}_R\left(x\right)=\frac{{\mathrm{IC1}}^0\left(x\right)/\mathrm{SE1}\left(x\right)}{\frac{1}{N}\sum _{i=1}^N{\mathrm{IC1}}^0\left(i\right)/\mathrm{SE1}\left(i\right)}\mathrm{for}x=1,\dots ,N& \mathrm{EQ}.36\end{array}$

[0149] As with the 2D case (FIG. 14), the shading errors can also be calculated from a right configuration as well. FIG. 16 and FIG. 18 show mat the right configuration is symmetrical to the left configuration. That is, the initial setup of a right calibration standard is the shifted-and-rotated arrangement 110 of a left calibration standard and vice versa. Therefore, there is no difference between a left configuration and a right configuration, except in the order they are given to an imaging system.

[0150] As with the 2D case (FIG. 14), the shift-and-rotation can be achieved by two shifts in the case where the detachable configuration is not convenient or impossible to fabricate.

[0151] We have presented both 1D and 2D cases. The same shift-and-rotation principle can be extended to the 3D (3-dimensional) domain as well. In the 3D domain, the field of view is 3 dimensional, and a calibration region becomes a calibration cube. The extension of a 2D calibration standard with a top-left configuration is a 3D calibration standard with a top-left-front configuration as shown in FIG. 19. The calibration standard 400 has N×M×L calibration cubes and has 8 detachable parts. Part 401 consists of (N−1)×(M−1)×(L−1) calibration cubes; part 402 has (N−1)×(L−1) calibration cubes; part 403, which is not visible (blocked by parts 406 and 408) in this figure, consists of (M−1)×(L−1) calibration cubes; part 404 consists of 1 calibration cube; part 405 consists of (N−1)×(M−1) calibration cubes; part 406 consists of (L−1) calibration cubes; part 407 consists of (N−1) calibration cubes; and part 408 consists of (M−1) calibration cubes. The configuration is named as top-left-front because of the position of part 404 is at top-left-front. A horizontally shifted-and-rotated arrangement of the calibration standard 400 will be done by shifting the entire standard 400 to the left by one column and then moving parts 403, 404, 406 and 408 to the right side of the parts 401, 407, 402 and 405 respectively. Similarly, the vertically shifted-and-rotated arrangement and the front-to-back shifted-and-rotated arrangement will be form in a similar manner.

[0152] The shading errors SE3(x, y, z) of the calibration cubes can be resolved from the 3D image intensities of the initial setup, the horizontally shifted-and-rotated arrangement, vertically shifted-and-rotated arrangement and the front-to-back shifted-and-rotated arrangement. All the mathematical deductions are similar to the 2D case.

[0153] Reference is again made to FIG. 3. The method of the present invention allows the log-transformed shading error (SE_L in EQ. 11 or EQ. 12) of each region 5 to be represented as a factor of the log-transformed shading errors of other regions 5. An intermediate shading error for each region may be calculated using standard matrix techniques, if the shading error for one region is specified (assuming that corner part 4 has a single region so that matrix A in EQ. 12 has a rank of N×M−1). The shading error may then be calculated as the normalized value of the intermediate shading error (EQ. 13 or EQ. 16). The effect of the selection in the intermediate shading errors is eliminated by the normalization.

[0154] Reference is made to EQ. 13, EQ. 16 and EQ. 33. Due to relative nature of shading correction, the normalization to derive a unique shading error is not completely necessary. In other words, the intermediate shading errors could be used to derive the shading correction coefficients. In this case, there may be an embedded system gain to the subsequently corrected images.

[0155] Referring to FIGS. 5 and 6, calibration standard 10 is re-arranged into horizontally shifted-and-rotated arrangements 11 and 12. Any other re-arrangements of calibration standard 10 which (i) allow the shading error of each region 5 in calibration standard 10 to be defined in terms of the shading errors of other regions and (ii) allow the shading error of each region 5 to be solved by specifying an shading error for one (or more) regions, may be used in the method of the present invention. A person skilled in the art will be capable of developing equations corresponding to EQ. 12 which allow the relative shading error for each region 5 to be solved.

[0156] Several embodiments have been described in details as examples of the present invention. These embodiments may be modified and varied without departing from the spirit and scope of the invention, which is limited only by the appended claims.