# A Combined Approach on RBC Image Segmentation through Shape Feature Extraction

Authors: Ruihu Wang; Bin Fang
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The classification of erythrocyte plays an important role in clinic diagnosis. In terms of the fact that the shape deformability of red blood cell brings more difficulty in detecting and recognize for operating automatically, we believed that the recovered 3D shape surface feature would give more information than traditional 2D intensity image processing methods. This paper proposed a combined approach for complex surface segmentation of red blood cell based on shape-from-shading technique and multiscale surface fitting. By means of the image irradiance equation under SEM imaging condition, the 3D height field could be recovered from the varied shading. Afterwards the depth maps of each point on the surfaces were applied to calculate Gaussian curvature and mean curvature, which were used to produce surface-type label image. Accordingly the surface was segmented into different parts through multiscale bivariate polynomials function fitting. The experimental results showed that this approach was easily implemented and promising.


---

## Body

## 1. Introduction

The erythrocyte shape deformability is critical to the filterability of blood. It has drawn considerable attentions into the pathology research in clinical relevant blood diseases. Unfortunately, the diagnosing is usually performed by a human expert, and it shows some drawbacks such as time-cost consuming and inaccuracy. Conventionally, the experts deal with erythrocyte images segmentation issue with 2D gray scale image. However, in order to obtain a satisfied performance, the classification and recognition should be based on the real shape of RBCs. In fact, the shape feature of red blood cell provides more useful information for diagnosing accurately than intensity level image. So it is necessary to take the shape information into consideration in this real problem. For our experiments, we use a database of 100 RBC images obtained by scanned electron microscope rather than optical imaging system. In [1], Egerton elaborated the operation principle of SEM, which creates images that are particularly easy to implement because the brightness in it is a function of the slope of the specimen at that point and forms a varied shading image. It is unlike the optical and transmission electron microscope, whose brightness depends on the thickness and optical and electron density instead.In [2], Russ mentioned that many of the two-dimensional images have been sectioned through three-dimensional structures. This is especially true in the various types of microscopy, where either polished flat planes or cut thin sections are needed in order to form the images in the first place. But the specimens thus sampled are three dimensional, and the goal of microscopists is to understand the three-dimensional structure. Some works have been done on optical blood cell images with traditional 2D methods [3–5]. In order to detect and classify malaria parasites in images of Giemsa-stained blood slides, Di Ruberto et al. proposed a morphological approach to evaluate the parasitaemia of the blood. They segmented the cells (red and white) from the background firstly and then detected and classified the parasites infecting them [3, 6]. Equally with malaria parasite detection mentioned above, Kumarasmy et al. presented an automated method for the robust analysis of RBC images via Gestalt laws [4], and Mandal et al. presented a segmentation method of blood smear images using normalized cut [5]. As we well know that, noise estimation is a challenge problem for complex structures of images, Liao et al. presented a method determining neighborhoods of the image pixels automatically with adaptive denoising and estimate noise for a single-slice sonogram of low-dose CT based on the homogenous patches centered at a special pixel. In their method, the noisy image is viewed as an observation of a nonlinear time series. The true state of the NTS must be recovered from the observation to realize image denoising [7, 8]. Furthermore, Hu et al. proposed an image smoothing using nonlinear anisotropic diffusion [9]. They suggested that the diffusion should be performed both among the time variants and spatial variant.Figure1 shows a typical example of such kind of red blood cell images captured by SEM with which we are going to deal. They were obtained at 600 times magnification using a scanned electron microscope. As shown in Figure 1, there are some outstanding characteristics about this image which make our problem be significant and challenging. On the one hand, the image represents very highly good quality with varied shading illuminated by light source. On the other hand, some light gridlines are superimposed in the image. The gridlines were added for the sake of manual counting. In addition, those cells’ shape takes on lots of irregular deformation, which is the primary problem we have to solve to segment effectively. As we well know, image segmentation is the bridge to classification properly. We aimed to develop a satisfied algorithm to classify the red blood cells into different groups accurately. And also we believed that conventional segmentation methods based on gray value could be unsuitable to this case. In this paper we proposed a new strategy to segment RBC image according to surface feature extraction. At first, we have to estimate the distribution of erythrocyte shapes from scanned electron microscope. Then each cell’s three-dimension shape was reconstructed as 3D height field using shape-from-shading technique. Lastly we implement multiscale surface fitting segmentation algorithm to partition the cells based on the depth data acquired in the previous procedures.Figure 1
A typical SEM image of red blood cells.This paper is organized as follow. Firstly the preliminary work including system framework and image preprocessing are introduced in Section2. In Section 3, a guided contour tracing method was used to extract the boundary and center point information. Accordingly all pixels of each cell which is located on the top of overlapped cells can be further obtained, whose intensity tone is disposed as shading information. The 3-D reconstruction of each cell is introduced in Section 4. We deduced an image irradiance equation under SEM imaging condition with linear approximation. Shape-from-shading technique is applied in this project, such as shape-from-shading technique using linear approximation. The next section is about to divide surface into several different types after computing mean curvature and Gaussian curvature. Multiscale surface fitting segmentation algorithm is proposed in Section 6 which involves seed extraction, region growing, and so on. In the end, Section 7 draws some conclusions and expects a few future works.


## 2. Image Preprocessing and System Framework

As shown in Figure2(a), there existed some bright gridlines superimposed on the original image, which have some side effects on the subsequent work. As mentioned before, the images show perfect quality other than these lines, which is used to count and classification manually. The system we developed here is to relieve human from exhausted hand work and run automatically. Additionally, in terms of the shading information being critical in our case, we regard the lines as noise and we have to remove them before recovering the 3D shape from gray tone image.RBC image preprocessing.
(a)
Gridline superimposed image(b)
Denoising using median locally(c)
Denoising using median filter directly on the whole image


### 2.1. System Framework

We describe the system framework in Algorithm1.Algorithm 1:
1. Image Denoising1.1 Find the rough position of each grid line;1.2 Remove the noise using median filtering locally.2. Cell Extraction Individually2.1 Guided boundary contour tracing;2.2 Cell Extraction through 4-neighborhood region growing.3. Shape from Shading Using Linear Approximation3.1 Image Irradiance Equation Deductio;3.2 3-D Shape Reconstructing.4. Computing RBC shape surface feature4.1 Compute mean curvature and Gaussian curvature using image convolution;4.2 Set a threshold for this two curvature.5: Segmentation through multiscale surface fittingMore details is in Section6: Experiment.


### 2.2. Image Processing Using Median Filtering Locally

As time consumption is sensitive during RBCs classification, we make use of median filter to get rid of the gridlines. The median filter is a smoothing technique that causes minimal edge blurring, which involves replacing the pixel value at each point in an image by the median of the pixel values in a neighborhood about the point:(2.1)g(x,y)=med{f(x-k,y-l),(k,l∈W)},
where f(x,y) is the original image and g(x,y) is median filtered image, respectively. W is a 2-dimensional 7*7 template.Figure2(a) shows a scaled original RBC image with four white gridlines superimposed on it. The denoised image after median filtering directly on whole image is presented in Figure 2(c), in which the lines are removed successfully. However, the edges of cell image have been blurred and brightness changed at the same time. Consequently the issue of inaccuracy would arise from the change, because the recovered shape is relied on the irradiance mostly. Fortunately we can detect the exact positions where those lines are by horizontal and vertical projection using(2.2)∑i⊕jQ(i,j),
respectively, where Q(i,j) only represents those pixels whose gray value is approximately equal to the pixels in those white gridlines. i⊕j means that only horizontal pixels are computed when sum on i and vertical pixels are projected when sum on j as well. Experimental results show that, while determining the exact positions of vertical lines, the points located in the range of i∈[253,259] or [509,517] have to be considered only. And those points whose j coordinate is in the range of j∈[253,259], [509,515], [765,771] have to be dealt with when impose median filtering locally.As a result, this local median filtering method leads to a handily approach, namely filling locally combined with median filtering, 7*7 structure element defined. The improved result is shown in Figure2(b), in which all the cells keep the same shading information as the original image and the gridlines have been removed successfully as well.


## 2.1. System Framework

We describe the system framework in Algorithm1.Algorithm 1:
1. Image Denoising1.1 Find the rough position of each grid line;1.2 Remove the noise using median filtering locally.2. Cell Extraction Individually2.1 Guided boundary contour tracing;2.2 Cell Extraction through 4-neighborhood region growing.3. Shape from Shading Using Linear Approximation3.1 Image Irradiance Equation Deductio;3.2 3-D Shape Reconstructing.4. Computing RBC shape surface feature4.1 Compute mean curvature and Gaussian curvature using image convolution;4.2 Set a threshold for this two curvature.5: Segmentation through multiscale surface fittingMore details is in Section6: Experiment.


## 2.2. Image Processing Using Median Filtering Locally

As time consumption is sensitive during RBCs classification, we make use of median filter to get rid of the gridlines. The median filter is a smoothing technique that causes minimal edge blurring, which involves replacing the pixel value at each point in an image by the median of the pixel values in a neighborhood about the point:(2.1)g(x,y)=med{f(x-k,y-l),(k,l∈W)},
where f(x,y) is the original image and g(x,y) is median filtered image, respectively. W is a 2-dimensional 7*7 template.Figure2(a) shows a scaled original RBC image with four white gridlines superimposed on it. The denoised image after median filtering directly on whole image is presented in Figure 2(c), in which the lines are removed successfully. However, the edges of cell image have been blurred and brightness changed at the same time. Consequently the issue of inaccuracy would arise from the change, because the recovered shape is relied on the irradiance mostly. Fortunately we can detect the exact positions where those lines are by horizontal and vertical projection using(2.2)∑i⊕jQ(i,j),
respectively, where Q(i,j) only represents those pixels whose gray value is approximately equal to the pixels in those white gridlines. i⊕j means that only horizontal pixels are computed when sum on i and vertical pixels are projected when sum on j as well. Experimental results show that, while determining the exact positions of vertical lines, the points located in the range of i∈[253,259] or [509,517] have to be considered only. And those points whose j coordinate is in the range of j∈[253,259], [509,515], [765,771] have to be dealt with when impose median filtering locally.As a result, this local median filtering method leads to a handily approach, namely filling locally combined with median filtering, 7*7 structure element defined. The improved result is shown in Figure2(b), in which all the cells keep the same shading information as the original image and the gridlines have been removed successfully as well.


## 3. Tracing Contour and Cell Extraction Individually

### 3.1. Guided Contour Tracing

Vromen and McCane [10] proposed a method named contour-tracing-based approach to the problem of finding the boundary of red blood cells in a scanned electron microscope image automatically. As shown in the above Figure 1, there are considerable overlapped cells. We are just only interested in estimating the distribution of different erythrocyte shapes from SEM image rather than the accurate counting number. So it makes sense to assume that the distribution of overlapped cells is identical to the overall distribution. Consequently only those top-level cells are needed to be detected and recognized. At the very beginning, the most possible direction is chosen by taking the prior information of tracing into account so far. Since the cell contours are likely elliptical, it would be reasonable to fit a conic shape to the path. A parameterized second-degree polynomial over the last n points was modeled the local curvature:(3.1)x(t)=at2+bt+c,y(t)=et2+ft+g.The best fitting polynomial was calculated through a number of data points using least squares. In order to increase accuracy and decrease computation, consider only the directions with angles in a certain window around the predicted direction. This was represented with a set of unit vectorui∈U:(3.2)U={(cos⁡(αi+β),sin⁡(αi+β))∣βmin⁡≤β≤βmax⁡},
where αi is the angle of the predicted tangent. In our application β is uniformly sampled between βmin=-(π/2) and βmax=(π/2).Figure3(a) shows the scaled 256*256 RBC image, and the contours after guided tracing are presented in Figure 3(b). There are 7 traced contours in Figure 3(b) which are located on top level of overlapped cell image. All of the contour information is stored as attributes in an XML file. It also contains the locus of each point on the boundary.Cell image extracted using contour tracing.
(a)
Scaled 256*256 image(b)
Contour after guided tracing(c)
Extracted cells


### 3.2. Cell Extraction Individually

As shown in Figure4, there is a break point marked by a white circle in each traced boundary, which would result in a wrong region when growing we should fill them as a complete boundary point before growing. If the number of neighbors around a predicted boundary point is less than three under 8-connectivity, we consider it as a break point. The relationships between the two break points and how they break with each other are shown in Figure 5.Figure 4
Break point in cells boundary.Five different types of break points.
(a)
Horizontal break(b)
Vertical break(c)
45° break(d)
135° break(e)
OthersAccording to the extracted cell contours’ information, we can grow each cell starting at center point regionally to get the entire cell image. The algorithms can be described in Algorithm2.Algorithm 2:Region growing algorithm.
1. Draw some key points for region growingCollect the center points and contour points from the extracted contour which is stored ina XML data file.2. Region Growing2.1 Set the center point as a seed;2.2 Grow the cell from the seed in the way of 4 connectivity neighborhoods till touch certainpoint on contour;2.3 Store the pixels’ information into a XML file attributed by position and gray level as treenodes.3. Continue to calculate cell pixels included in next contour4. If All cells have been extracted successfully, then exit5. Else return 1The resulting image of region growing subject to contour boundary is shown in Figure3(c). After growing regionally, we got the number of pixels which make the whole image and their gray level value. In Table 1, there are 7 cell contours that have been extracted altogether, where CENTER_X and CENTER_Y are the cell’s center point coordinates and PIXELS denotes the pixels number involved in each cell.Table 1
Contours’ information of traced cells.
No.Center_XCenter_YPixels137.884175.93126972114.131179.52629273191.34665.334827774187.606138.71825165154.363209.66720276211.253202.63423597233.308134.118953


## 3.1. Guided Contour Tracing

Vromen and McCane [10] proposed a method named contour-tracing-based approach to the problem of finding the boundary of red blood cells in a scanned electron microscope image automatically. As shown in the above Figure 1, there are considerable overlapped cells. We are just only interested in estimating the distribution of different erythrocyte shapes from SEM image rather than the accurate counting number. So it makes sense to assume that the distribution of overlapped cells is identical to the overall distribution. Consequently only those top-level cells are needed to be detected and recognized. At the very beginning, the most possible direction is chosen by taking the prior information of tracing into account so far. Since the cell contours are likely elliptical, it would be reasonable to fit a conic shape to the path. A parameterized second-degree polynomial over the last n points was modeled the local curvature:(3.1)x(t)=at2+bt+c,y(t)=et2+ft+g.The best fitting polynomial was calculated through a number of data points using least squares. In order to increase accuracy and decrease computation, consider only the directions with angles in a certain window around the predicted direction. This was represented with a set of unit vectorui∈U:(3.2)U={(cos⁡(αi+β),sin⁡(αi+β))∣βmin⁡≤β≤βmax⁡},
where αi is the angle of the predicted tangent. In our application β is uniformly sampled between βmin=-(π/2) and βmax=(π/2).Figure3(a) shows the scaled 256*256 RBC image, and the contours after guided tracing are presented in Figure 3(b). There are 7 traced contours in Figure 3(b) which are located on top level of overlapped cell image. All of the contour information is stored as attributes in an XML file. It also contains the locus of each point on the boundary.Cell image extracted using contour tracing.
(a)
Scaled 256*256 image(b)
Contour after guided tracing(c)
Extracted cells


## 3.2. Cell Extraction Individually

As shown in Figure4, there is a break point marked by a white circle in each traced boundary, which would result in a wrong region when growing we should fill them as a complete boundary point before growing. If the number of neighbors around a predicted boundary point is less than three under 8-connectivity, we consider it as a break point. The relationships between the two break points and how they break with each other are shown in Figure 5.Figure 4
Break point in cells boundary.Five different types of break points.
(a)
Horizontal break(b)
Vertical break(c)
45° break(d)
135° break(e)
OthersAccording to the extracted cell contours’ information, we can grow each cell starting at center point regionally to get the entire cell image. The algorithms can be described in Algorithm2.Algorithm 2:Region growing algorithm.
1. Draw some key points for region growingCollect the center points and contour points from the extracted contour which is stored ina XML data file.2. Region Growing2.1 Set the center point as a seed;2.2 Grow the cell from the seed in the way of 4 connectivity neighborhoods till touch certainpoint on contour;2.3 Store the pixels’ information into a XML file attributed by position and gray level as treenodes.3. Continue to calculate cell pixels included in next contour4. If All cells have been extracted successfully, then exit5. Else return 1The resulting image of region growing subject to contour boundary is shown in Figure3(c). After growing regionally, we got the number of pixels which make the whole image and their gray level value. In Table 1, there are 7 cell contours that have been extracted altogether, where CENTER_X and CENTER_Y are the cell’s center point coordinates and PIXELS denotes the pixels number involved in each cell.Table 1
Contours’ information of traced cells.
No.Center_XCenter_YPixels137.884175.93126972114.131179.52629273191.34665.334827774187.606138.71825165154.363209.66720276211.253202.63423597233.308134.118953


## 4. Shape from Shading Using Linear Approximation

### 4.1. Shape from Shading and PDE

The “shape from shading” problem, namely, SFS, is to recover the 3-D shape of a surface from a gray-level monochrome image. In the 1970s Horn firstly proposed the approach to reconstruct the original shape from a varied shading image, which associated with obtaining a solution of a nonlinear first-order partial differential equation (PDE), that is, brightness equation. From then on, a number of articles have emerged which come up with various kinds of methods to strive to implement this technique into real or artificial synthetic images.This PDE equation arises from the(4.1)I(x1,x2)=R(n(x1,x2)).(x1,x2) is the coordinates of a point x in the image. The brightness equation connects the reflectance map (R) to the brightness image (I). Almost all the shape-from-shading methods at the exception of an extremely small number of papers [11–13] assume that the scene model is Lambertian. The reflectance map is the cosine of the angle between the light vector L(x) and the normal vector n(x) to the surface:(4.2)R=cos⁡⁡(L,n)=L|L|⋅n|n|,
where R, L, and n depend on (x1,x2) [14].Shape from shading is a fundamental issue in computer vision, and considerable research has been performed [15, 16] in trying to solve this problem including methods of medical image processing. In [17], the authors applied their method to an endoscopic image of a normal stomach and showed the result obtained by generic algorithm in the perspective case with the light source at the optical center, which is not suitable for SEM. Tankus et al. in their papers [18–20] suggested the reconstruction algorithm under an assumption of perspective shape from shading. Deguchi and Okatani [21] accomplished shape reconstruction from an endoscope image by shape-from-shading technique for a point light source at the projection center.


### 4.2. Reflectance Map under SEM Imaging Condition

Jones and Taylor [22] proposed that SEM imaging process is particularly appropriate for SFS, since it allows us to make the simplifying assumptions that the projection is orthographic and the "light source" is at infinity. The Lambertian reflectance function is given by(4.3)R(p,q)=ηn̂⋅Î,
where n̂=(-p,-q,1)/p2+q2+1 is the unit normal vector, Î is a unit vector in the direction of the light source, and η is the surface albedo.The SEM reflectance function is based upon the theoretical prediction that the number of electron emitted from a surface in the SEM is proportional to the secant of the angle between the illumination direction and the surface normal. The reflectance function is denoted by [22](4.4)R(p,q)=ηn̂⋅Î.In this paper, we assume that the image is formed by orthographic projection, because the object specimen examined by SEM is very small in comparison to distance from light source.


### 4.3. Linear Approximation

In [23], the authors believed that the linearity of the reflectance map in the depth Z, instead of p and q, is more appropriate in some cases. They presented a method for computing depth from a single-shaded image by employing the discrete approximations for p and q using finite differences and linearly approximating the reflectance in Z(x,y). It gave good results for the spherical surface and can be applied to any reflectance function.In this paper we aim to recover the red blood cell’s shape of an image captured by SEM based on linear approximation. It is extended to solve such a problem assumed orthographic projection and derive the implementation equations with the reflectance function inverse to Lambertian reflectance function.Image irradiance equation (IRE) indicates the relationship between reflection function and image irradiance. The recovered shape can be represented by depth map Z, normal (nx,ny,nz), or surface gradient (p,q). The radiance of surface patch depends on gradient, light source location, and reflectance property. The gray level of a pixel in the image is determined by light direction and normal vector, assumed Lambertian model, which can be denoted by IRE:(4.5)E(x,y)=R(p,q)=1+pps+qqs1+p2+q21+ps2+qs2,
where E(x,y) is a gray level at pixel (x,y),p=∂z/∂x and q=∂z/∂y, and (ps,qs,1) is the illumination direction.By approximating thep and q discretely, we get(4.6)p=∂z∂x=Z(x,y)-Z(x-1,y),q=∂z∂y=Z(x,y)-Z(x,y-1).
According to (4.6), the reflection function can be rewritten as(4.7)0=f(E(x,y),Z(x,y),Z(x-1,y),Z(x,y-1))=E(x,y)-R(Z(x,y)-Z(x-1,y),Z(x,y)-Z(x,y-1)).Under the assumption of the point(x,y) and image E is given, the linear approximation of function f with respect to Zn-1 after Jacobi  iteration method is:(4.8)0=f(Z(x,y))≈f(Zn-1(x,y))+(Z(x,y)-Zn-1(x,y))ddZ(x,y)f(Zn-1(x,y)).
The nth iterative result Zn(x,y) can be denoted by(4.9)Zn(x,y)=Zn-1(x,y)+-f(Zn-1(x,y))(d/dZ(x,y))f(Zn-1(x,y)).As mentioned previously, the reflection function is inversal to Lambertian model under the condition of SEM imaging. Equation (4.5) is transformed into(4.10)R(p,q)=1+p2+q21+ps2+qs21+pps+qqs.
Now we compute the partial derivatives of p and q with R:(4.11)∂R(p,q)∂p=1+ps2+qs2p(1+pps+qqs)-(1+p2+q2)⋅ps(1+pps+qqs)21+p2+q2,∂R(p,q)∂q=1+ps2+qs2q(1+pps+qqs)-(1+p2+q2)⋅qs(1+pps+qqs)21+p2+q2.
So,(4.12)∂R∂p+∂R∂q=1+ps2+qs2((p+q)(1+pps+qqs)-(1+p2+q2)⋅(ps+qs))(1+pps+qqs)21+p2+q2.The right part of (4.9) is rewritten as(4.13)ddZ(x,y)f(Zn-1(x,y))=∂R∂p+∂R∂q=1+ps2+qs2((p+q)(1+pps+qqs)-(1+p2+q2)⋅(ps+qs))(1+pps+qqs)21+p2+q2.We use the shape-from-shading method with linear approximation to reconstruct the red blood cell’s 3-D shape as in Figure6.Reconstructed RBC 3D shape.
(a)
Original shape(b)
Filtered smooth shape


## 4.1. Shape from Shading and PDE

The “shape from shading” problem, namely, SFS, is to recover the 3-D shape of a surface from a gray-level monochrome image. In the 1970s Horn firstly proposed the approach to reconstruct the original shape from a varied shading image, which associated with obtaining a solution of a nonlinear first-order partial differential equation (PDE), that is, brightness equation. From then on, a number of articles have emerged which come up with various kinds of methods to strive to implement this technique into real or artificial synthetic images.This PDE equation arises from the(4.1)I(x1,x2)=R(n(x1,x2)).(x1,x2) is the coordinates of a point x in the image. The brightness equation connects the reflectance map (R) to the brightness image (I). Almost all the shape-from-shading methods at the exception of an extremely small number of papers [11–13] assume that the scene model is Lambertian. The reflectance map is the cosine of the angle between the light vector L(x) and the normal vector n(x) to the surface:(4.2)R=cos⁡⁡(L,n)=L|L|⋅n|n|,
where R, L, and n depend on (x1,x2) [14].Shape from shading is a fundamental issue in computer vision, and considerable research has been performed [15, 16] in trying to solve this problem including methods of medical image processing. In [17], the authors applied their method to an endoscopic image of a normal stomach and showed the result obtained by generic algorithm in the perspective case with the light source at the optical center, which is not suitable for SEM. Tankus et al. in their papers [18–20] suggested the reconstruction algorithm under an assumption of perspective shape from shading. Deguchi and Okatani [21] accomplished shape reconstruction from an endoscope image by shape-from-shading technique for a point light source at the projection center.


## 4.2. Reflectance Map under SEM Imaging Condition

Jones and Taylor [22] proposed that SEM imaging process is particularly appropriate for SFS, since it allows us to make the simplifying assumptions that the projection is orthographic and the "light source" is at infinity. The Lambertian reflectance function is given by(4.3)R(p,q)=ηn̂⋅Î,
where n̂=(-p,-q,1)/p2+q2+1 is the unit normal vector, Î is a unit vector in the direction of the light source, and η is the surface albedo.The SEM reflectance function is based upon the theoretical prediction that the number of electron emitted from a surface in the SEM is proportional to the secant of the angle between the illumination direction and the surface normal. The reflectance function is denoted by [22](4.4)R(p,q)=ηn̂⋅Î.In this paper, we assume that the image is formed by orthographic projection, because the object specimen examined by SEM is very small in comparison to distance from light source.


## 4.3. Linear Approximation

In [23], the authors believed that the linearity of the reflectance map in the depth Z, instead of p and q, is more appropriate in some cases. They presented a method for computing depth from a single-shaded image by employing the discrete approximations for p and q using finite differences and linearly approximating the reflectance in Z(x,y). It gave good results for the spherical surface and can be applied to any reflectance function.In this paper we aim to recover the red blood cell’s shape of an image captured by SEM based on linear approximation. It is extended to solve such a problem assumed orthographic projection and derive the implementation equations with the reflectance function inverse to Lambertian reflectance function.Image irradiance equation (IRE) indicates the relationship between reflection function and image irradiance. The recovered shape can be represented by depth map Z, normal (nx,ny,nz), or surface gradient (p,q). The radiance of surface patch depends on gradient, light source location, and reflectance property. The gray level of a pixel in the image is determined by light direction and normal vector, assumed Lambertian model, which can be denoted by IRE:(4.5)E(x,y)=R(p,q)=1+pps+qqs1+p2+q21+ps2+qs2,
where E(x,y) is a gray level at pixel (x,y),p=∂z/∂x and q=∂z/∂y, and (ps,qs,1) is the illumination direction.By approximating thep and q discretely, we get(4.6)p=∂z∂x=Z(x,y)-Z(x-1,y),q=∂z∂y=Z(x,y)-Z(x,y-1).
According to (4.6), the reflection function can be rewritten as(4.7)0=f(E(x,y),Z(x,y),Z(x-1,y),Z(x,y-1))=E(x,y)-R(Z(x,y)-Z(x-1,y),Z(x,y)-Z(x,y-1)).Under the assumption of the point(x,y) and image E is given, the linear approximation of function f with respect to Zn-1 after Jacobi  iteration method is:(4.8)0=f(Z(x,y))≈f(Zn-1(x,y))+(Z(x,y)-Zn-1(x,y))ddZ(x,y)f(Zn-1(x,y)).
The nth iterative result Zn(x,y) can be denoted by(4.9)Zn(x,y)=Zn-1(x,y)+-f(Zn-1(x,y))(d/dZ(x,y))f(Zn-1(x,y)).As mentioned previously, the reflection function is inversal to Lambertian model under the condition of SEM imaging. Equation (4.5) is transformed into(4.10)R(p,q)=1+p2+q21+ps2+qs21+pps+qqs.
Now we compute the partial derivatives of p and q with R:(4.11)∂R(p,q)∂p=1+ps2+qs2p(1+pps+qqs)-(1+p2+q2)⋅ps(1+pps+qqs)21+p2+q2,∂R(p,q)∂q=1+ps2+qs2q(1+pps+qqs)-(1+p2+q2)⋅qs(1+pps+qqs)21+p2+q2.
So,(4.12)∂R∂p+∂R∂q=1+ps2+qs2((p+q)(1+pps+qqs)-(1+p2+q2)⋅(ps+qs))(1+pps+qqs)21+p2+q2.The right part of (4.9) is rewritten as(4.13)ddZ(x,y)f(Zn-1(x,y))=∂R∂p+∂R∂q=1+ps2+qs2((p+q)(1+pps+qqs)-(1+p2+q2)⋅(ps+qs))(1+pps+qqs)21+p2+q2.We use the shape-from-shading method with linear approximation to reconstruct the red blood cell’s 3-D shape as in Figure6.Reconstructed RBC 3D shape.
(a)
Original shape(b)
Filtered smooth shape


## 5. Curvature Calculation

There are 8 different types of surface altogether, namely peak, pit, ridge, valley, flat, minimal surface, saddle ridge, and saddle valley. The surface type of each data point on a scene object can be designated by the signs of mean curvature and Gaussian curvature uniquely. Both of these two curvatures can be calculated by local convolution [24, 25]. Each data point in a given window N*N is associated with a 2-dimensional position (u,v) from the set U×U, where(5.1)U={-N-12,…,-1,0,1,…,N-12}
and N is odd.The following discrete orthogonal polynomials provide local biquadratic surface fitting capability:(5.2)ϕ0(u)=1,ϕ1(u)=u,ϕ2(u)=u2-M(M+1)3,
where M=(N-1)/2. The biquadratic is the minimal degree polynomial surface type needed to estimate the first and second partial derivatives. A corresponding set of bi(u) functions is the normalized versions of the orthogonal polynomials ϕi(u) given by bi(u)=ϕi(u)/Pi(M), where the Pi(M) are normalizing constants. The three normalization constants are given by(5.3)P0(M)=N,P1(M)=23M3+M2+13M,P2(M)=845M5+49M4+29M3-19M2-115M.Defining a set of surfaceG in R3 space, G⊂R3, which can be parameterized by(5.4)G={g(u,v)=[g1(u,v)g2(u,v)g3(u,v)]:vmin⁡<v<vmax⁡umin⁡<u<umax⁡}.Before implementing the surface-fitting segmentation algorithm, the surface types have to be divided at first, which is based on mean curvature and Gaussian curvature. The computation of mean curvature and Gaussian curvature of digital surface is approximated by partial derivative estimation, which is calculated via the appropriate 2D image convolution (denoted by *):(5.5)gu=DuS*g*,gv=DvS*g*,guu=DuuS*g*,guv=DuvS*g*,gvv=DvvS*g*,
where S=s⃗s⃗T is a 7×7 binomial smoothing window:(5.6)s⃗=164[1615201561],[Du]=d⃗0d⃗1T,[Dv]=d⃗1d⃗0T,[Duu]=d⃗0d⃗2T,[Duv]=d⃗1d⃗1T,[Dvv]=d⃗2d⃗0T,
where(5.7)d⃗0=17[1111111]T,d⃗1=128[-3-2-10123]T,d⃗2=184[50-3-4-305]T.The mean curvature and Gaussian curvature can be calculated by partial derivative as follows:(5.8)H=(1+gv2)guu+(1+gu2)gvv-2gugvguv2(1+gu2+gv2)3,K=guugvv-guv2(1+gu2+gv2)2.
Table 2 shows the average and variance of mean curvature and Gaussian curvature of the 7 extracted cells in Figure 3.Table 2
Average and variance of mean curvature and gaussian curvature.
No.Average of MCAverage of GCVariance of MCVariance of GC1−0.00358235−0.01188670.01444670.011089320.0019092−0.00319610.00074670.001353113−0.00801067−0.0001365840.007765778.05916e−0640.00362456−0.004921720.001742060.0028007454.11563e-05−0.005027110.002811720.0067241360.00774488−0.002439660.00715360.00037371270.00235231−0.01329020.001920860.0311508


## 6. Experiment

### 6.1. Surface-Type-Based Image Segmentation

The fundamental formulation of region-based image segmentation is defined as(1)
Ui=1nRi=R;(2)
Ri is a connected region, i=1,2,…,n;(3)
Ri∩Rj=Φ,for∀(i,j),i≠j;(4)
P(Ri)=TRUE,i=1,2,…,n;(5)
P(Ri∪Rj)=FALSE, Ri is adjacent to Rj,whereP(Ri) is a uniformity predicate defined on groups of connected pixels. Ri was grown regionally via 8-connected neighborhood. All the points in region Ri satisfy the same surface function. Different regions meet different surface fitting function.The segmentation procedure is divided into two mainly different parts. Firstly we compute the surface-type label image by(6.1)T=1+3(1+3sgnεH(H))+(1-3sgnεK(K)),
where T denotes the surface type ranging from 1 to 9 as shown in Table 3.Table 3
Surface type defined by mean and Gaussian curvature.
K>0K=0K<0H<0Peak(T=1)Ridge(T=2)Saddle ridge(T=3)H=0None(T=4)Flat(T=5)Minimal surface(T=6)H>0Pit(T=7)Valley(T=8)Saddle valley(T=9)We the define root mean square error (RMSE) as(6.2)σSN2=1N∑(x,y)∈SN(zij-Φa⃗(x,y))2,
to measure whether the difference between fitted value and original depth is confined to the range preset by a threshold as in Algorithm 3.Algorithm 3:Algorithm of multiscale segmentation through mapped depth.
1. Initialization1.1 Get the whole depth map image from the output of shape from shading procedure;1.2 Obtaining the range image by filling the value of each pixel, which is associated with thecurrent processing cell by depth data;1.3 Compute an estimate of the noise variance at each pixel;1.4 Computing mean curvature and Gaussian curvature through separable convolution;1.5 Computing the surface type label image and find all connected components of each surfacetype label image, sort it to get histogram distribution;1.6 Extracted seed region through erosion (contraction) operation.2. Iterative variable order surface fitting2.1 Perform surface fit from the lowest order, if it is OK using RMS error and region test;2.2 Then goto 3;2.3 Else increase the order and fit again;2.4 if order >4, then return.3. Region Growing3.1 Find the new region consisting of compatible connected neighboring pixels.


### 6.2. Algorithm

#### 6.2.1. Experimental Result

In our experiment, we define the RMS fit error asε=ωσimg, where σimg means noise variance, and(6.3)ẑ(p)=f̂(mk,a⃗l,x(p),y(p))
is compared with(6.4)z(p)=g̃(x(p),y(p)),
to see if the pixel p is compatible with the approximating surface function. If the magnitude of the difference between the function value and the digital surface value is less than the allowed tolerance value, denoted by ω0·εlk, then the pixel p is added to the set of compatible pixels, denoted by C(mk,a⃗lk,εlk), which are compatible with the surface fit to the region R̂lk. Otherwise, the pixel is incompatible and discarded. The result of this process is the compatible pixel list: (6.5)C(mk,a⃗lk,εlk)=p∈I:|ẑ(p)-z(p)|≤ω0εlk.
We choose ωl=4.5 and ω0=8 experimentally.If we threshold the mean curvature and Gaussian curvature with an appropriate value, only three different surface type are remained among all cells. In our experiment, we choose(6.6)εH=0.0015,εK=0.03
as the mean curvature and Gaussian curvature separately. The function of sgn in (6.1) is defined as(6.7)sgn⁡(H,K)={1(ε,+∞)0(-ε,+ε)-1(-∞,-ε).The three kinds of surface type are flat, pit, and valley, respectively, as shown in Figure7.Three surface-type components.
(a)
Flat(b)
Pit(c)
ValleyFigure8 represents the process of extracting seed region through erosion (contraction) operation. In order to implement the surface fitting segmentation algorithm, the seed region has to be obtained using erosion operation firstly. After several iteration, the pixel numbers of remained regions for growing are 32, 34, 45, 6 for these four region, respectively.Figure 8
Contraction process to extract seed region of red blood cell.In Figure9, the cell is segmented into three isolated parts perfectly, which are obtained through fitting based on surface type.The segmentation result using surface fitting method.
(a)
Original cell image(b)
Flat region grew regionally(c)
Pit region 1 grew regionally(d)
Pit region 2 grew regionally


### 6.3. Evaluation

In order to evaluate our proposed combined algorithm, we used a dataset containing 800 SEM images. These images are with resolution of 1024 by 768 pixels. For our evaluation, we ran the algorithm on 100 randomly selected images.We divided the cells into different categories according to their distribution of surface type. Table4 elaborated the segmentation accuracy of each category.Table 4
Accuracy of segmentation.
CategoryNumber of flat surfacesNumber of pit surfacesNumber of valley SurfacesAccuracy111088.9%212093.3%321090.9%4010100%5330100%6221100%7310100%822094.6%9111100%10211100%11121100%


## 6.1. Surface-Type-Based Image Segmentation

The fundamental formulation of region-based image segmentation is defined as(1)
Ui=1nRi=R;(2)
Ri is a connected region, i=1,2,…,n;(3)
Ri∩Rj=Φ,for∀(i,j),i≠j;(4)
P(Ri)=TRUE,i=1,2,…,n;(5)
P(Ri∪Rj)=FALSE, Ri is adjacent to Rj,whereP(Ri) is a uniformity predicate defined on groups of connected pixels. Ri was grown regionally via 8-connected neighborhood. All the points in region Ri satisfy the same surface function. Different regions meet different surface fitting function.The segmentation procedure is divided into two mainly different parts. Firstly we compute the surface-type label image by(6.1)T=1+3(1+3sgnεH(H))+(1-3sgnεK(K)),
where T denotes the surface type ranging from 1 to 9 as shown in Table 3.Table 3
Surface type defined by mean and Gaussian curvature.
K>0K=0K<0H<0Peak(T=1)Ridge(T=2)Saddle ridge(T=3)H=0None(T=4)Flat(T=5)Minimal surface(T=6)H>0Pit(T=7)Valley(T=8)Saddle valley(T=9)We the define root mean square error (RMSE) as(6.2)σSN2=1N∑(x,y)∈SN(zij-Φa⃗(x,y))2,
to measure whether the difference between fitted value and original depth is confined to the range preset by a threshold as in Algorithm 3.Algorithm 3:Algorithm of multiscale segmentation through mapped depth.
1. Initialization1.1 Get the whole depth map image from the output of shape from shading procedure;1.2 Obtaining the range image by filling the value of each pixel, which is associated with thecurrent processing cell by depth data;1.3 Compute an estimate of the noise variance at each pixel;1.4 Computing mean curvature and Gaussian curvature through separable convolution;1.5 Computing the surface type label image and find all connected components of each surfacetype label image, sort it to get histogram distribution;1.6 Extracted seed region through erosion (contraction) operation.2. Iterative variable order surface fitting2.1 Perform surface fit from the lowest order, if it is OK using RMS error and region test;2.2 Then goto 3;2.3 Else increase the order and fit again;2.4 if order >4, then return.3. Region Growing3.1 Find the new region consisting of compatible connected neighboring pixels.


## 6.2. Algorithm

### 6.2.1. Experimental Result

In our experiment, we define the RMS fit error asε=ωσimg, where σimg means noise variance, and(6.3)ẑ(p)=f̂(mk,a⃗l,x(p),y(p))
is compared with(6.4)z(p)=g̃(x(p),y(p)),
to see if the pixel p is compatible with the approximating surface function. If the magnitude of the difference between the function value and the digital surface value is less than the allowed tolerance value, denoted by ω0·εlk, then the pixel p is added to the set of compatible pixels, denoted by C(mk,a⃗lk,εlk), which are compatible with the surface fit to the region R̂lk. Otherwise, the pixel is incompatible and discarded. The result of this process is the compatible pixel list: (6.5)C(mk,a⃗lk,εlk)=p∈I:|ẑ(p)-z(p)|≤ω0εlk.
We choose ωl=4.5 and ω0=8 experimentally.If we threshold the mean curvature and Gaussian curvature with an appropriate value, only three different surface type are remained among all cells. In our experiment, we choose(6.6)εH=0.0015,εK=0.03
as the mean curvature and Gaussian curvature separately. The function of sgn in (6.1) is defined as(6.7)sgn⁡(H,K)={1(ε,+∞)0(-ε,+ε)-1(-∞,-ε).The three kinds of surface type are flat, pit, and valley, respectively, as shown in Figure7.Three surface-type components.
(a)
Flat(b)
Pit(c)
ValleyFigure8 represents the process of extracting seed region through erosion (contraction) operation. In order to implement the surface fitting segmentation algorithm, the seed region has to be obtained using erosion operation firstly. After several iteration, the pixel numbers of remained regions for growing are 32, 34, 45, 6 for these four region, respectively.Figure 8
Contraction process to extract seed region of red blood cell.In Figure9, the cell is segmented into three isolated parts perfectly, which are obtained through fitting based on surface type.The segmentation result using surface fitting method.
(a)
Original cell image(b)
Flat region grew regionally(c)
Pit region 1 grew regionally(d)
Pit region 2 grew regionally


## 6.2.1. Experimental Result

In our experiment, we define the RMS fit error asε=ωσimg, where σimg means noise variance, and(6.3)ẑ(p)=f̂(mk,a⃗l,x(p),y(p))
is compared with(6.4)z(p)=g̃(x(p),y(p)),
to see if the pixel p is compatible with the approximating surface function. If the magnitude of the difference between the function value and the digital surface value is less than the allowed tolerance value, denoted by ω0·εlk, then the pixel p is added to the set of compatible pixels, denoted by C(mk,a⃗lk,εlk), which are compatible with the surface fit to the region R̂lk. Otherwise, the pixel is incompatible and discarded. The result of this process is the compatible pixel list: (6.5)C(mk,a⃗lk,εlk)=p∈I:|ẑ(p)-z(p)|≤ω0εlk.
We choose ωl=4.5 and ω0=8 experimentally.If we threshold the mean curvature and Gaussian curvature with an appropriate value, only three different surface type are remained among all cells. In our experiment, we choose(6.6)εH=0.0015,εK=0.03
as the mean curvature and Gaussian curvature separately. The function of sgn in (6.1) is defined as(6.7)sgn⁡(H,K)={1(ε,+∞)0(-ε,+ε)-1(-∞,-ε).The three kinds of surface type are flat, pit, and valley, respectively, as shown in Figure7.Three surface-type components.
(a)
Flat(b)
Pit(c)
ValleyFigure8 represents the process of extracting seed region through erosion (contraction) operation. In order to implement the surface fitting segmentation algorithm, the seed region has to be obtained using erosion operation firstly. After several iteration, the pixel numbers of remained regions for growing are 32, 34, 45, 6 for these four region, respectively.Figure 8
Contraction process to extract seed region of red blood cell.In Figure9, the cell is segmented into three isolated parts perfectly, which are obtained through fitting based on surface type.The segmentation result using surface fitting method.
(a)
Original cell image(b)
Flat region grew regionally(c)
Pit region 1 grew regionally(d)
Pit region 2 grew regionally


## 6.3. Evaluation

In order to evaluate our proposed combined algorithm, we used a dataset containing 800 SEM images. These images are with resolution of 1024 by 768 pixels. For our evaluation, we ran the algorithm on 100 randomly selected images.We divided the cells into different categories according to their distribution of surface type. Table4 elaborated the segmentation accuracy of each category.Table 4
Accuracy of segmentation.
CategoryNumber of flat surfacesNumber of pit surfacesNumber of valley SurfacesAccuracy111088.9%212093.3%321090.9%4010100%5330100%6221100%7310100%822094.6%9111100%10211100%11121100%


## 7. Conclusion and Future Work

This paper is about how to reconstruct the 3D shape of red blood cell from gray tone images using scanned electronic microscope based on shape-from-shading technique, as well combined with linear approximation. The result of cell surface shape is given by height field. Our algorithm can be trivially transformed to various different kinds of reflection models. In Figure7, the surface-type label image is given with cell number added manually. There are mainly three types of surfaces left after threshold. The distribution of count number of each surface type in every cell can present some useful information for classifying correctly, which will be trained as input data. In the end, we aim to construct a classifier by means of cascaded SVMs architecture to recognize whether the red blood cell is normal or not.



---

*Source: 194953-2012-03-29.xml*