GUANGCHUN

Understanding Harris corner detector

February 16, 2013 gccheng 5 comments

Let’s first go over Harris detector a little bit. For a basic idea about Harris detector, check textbooks or opencv or blogs. For its intuition, check its precursor Movarec operator, which explains why we want to maximize the variation within a window to find a corner. The key to Harris detector is the variation of intensity within a sliding window $w(x,y)$ (with displacement $u$ in the x direction and $v$ in the y direction)

$E(u,v) = \sum_{x,y} w(x,y)[I(x+u, y+v) - I(x,y)]^2$

where:

$w(x,y)$ is the window weighting function at position $(x,y)$
$I(x,y)$ is the intensity at (x,y)
$I(x+u, y+v)$ is the intensity at the moved window $(x+u, y+v)$

The aim is to maximize $E(u,v)$ . Applying Taylor expansion and some arithmetic operations, we can get

$E(u,v) \approx \sum_{x,y} w(x,y) (u^2I_x^2+2uvI_xI_y+v^2I_y^2) \\ = (u,v)\left( \sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2\end{bmatrix}\right) \begin{pmatrix}u\\v \end{pmatrix} \\ =(u,v)M\begin{pmatrix}u\\v \end{pmatrix}$

where $(I_x, I_y)$ is the gradient at $(x,y)$ . Matrix $M$ is known as structure tensor of a pixel. It actually a characterization of information of all pixels within the window. Then we calculate the eigenvalues $\lambda_1$ and $\lambda_2$ of matrix $M$ to determine if window corresponds to a corner. (I)

If $\lambda_1 \approx 0$ and $\lambda_2 \approx 0$ then this pixel $(x,y)$ has no features of interest.
If $\lambda_1 \approx 0$ and $\lambda_2$ has some large positive value, then an edge is found.
If $\lambda_1$ and $\lambda_2$ have large positive values, then a corner is found.

[Note: to avoid computing the eigenvalues which is computationally expensive, in original Harris detector, the authors suggested using the following score:

$R = det(M) - k(trace(M))^{2}$

to determine if a point is a corner.]

The question is WHY? Why does the matrix $M$ have the capacity to determine the cornerness? Here, we look around to the covariance matrix $\Sigma$ of gradient $(I_x, I_y)$ at a pixel. Let $\mu_{x}=E[ I_{x} ],\mu_{y}=E[ I_{y}]$ . By definition of covariance matrix, we have

$\begin{array}{rl} \Sigma & = \begin{bmatrix} E[ (I_{x}-\mu_{x})(I_{x}-\mu_{x}) ] & E[ (I_{x}-\mu_{x})(I_{y}-\mu_{y}) ] \\ E[ (I_{y}-\mu_{y})(I_{x}-\mu_{x}) ] & E[ (I_{y}-\mu_{y})(I_{y}-\mu_{y}) ] \\ \end{bmatrix} \\ & = \begin{bmatrix} E[ I_{x}I_{x} ] - \mu_{x}^{2} & E[ I_{x}I_{y} ] - \mu_{x}\mu_{y}\\ E[ I_{y}I_{x} ] - \mu_{y}\mu_{x} & E[ I_{y}I_{y} ] - \mu_{y}^{2}\\ \end{bmatrix} \\ & \triangleq \begin{bmatrix} E[ I_{x}I_{x} ] & E[ I_{x}I_{y} ]\\ E[ I_{y}I_{x} ] & E[ I_{y}I_{y} ]\\ \end{bmatrix}\\ \end{array}$

The last step is derived based on the the assumption that the mean of derivatives $I_x$ and $I_y$ is zeros, i.e. $(\mu_{x},\mu_{y})=0$ .

On the other hand, the structure tensor $M$ is as follows if we apply an averaging filter over the time dimension (= we give uniform weight over the pixels within the pixels):

$\begin{array}{rl} M & = \frac{1}{f}\sum_{x,y} \begin{bmatrix} I_{x}^2 & I_{x}I_{y}\\ I_{x}I_{y} & I_{y}^2 \\ \end{bmatrix} \\ & = \begin{bmatrix} \widehat{I_{x}^2} & \widehat{I_{x}I_{y}}\\ \widehat{I_{x}I_{y}} & \widehat{I_{y}^2} \\ \end{bmatrix} \\ \end{array}$

where $f$ is total number of pixels within the window, and $\widehat{f}$ is the average of $f$ .

It can be seen that structure tensor $M$ is actually an unbiased estimate of the covariance matrix of the gradients for the pixels within the window. Therefore, any findings/conclusions about covariance matrix could be used to analyze structure tensors.

Let’s look at an example showing how the derivatives are distributed for different types of image patches (i.e. the window). The following figures are adapted from Robert Collins‘s slides. In this example, the image patches are from three different types: linear edge, flat area and corner. You can image that we are testing if the pixel at the center is a corner (feature point).

The first step is calculating the derivatives $I_x$ and $I_y$ for each pixel in the image, then the distributions of ( $I_x$ , $I_y$ ) are shown in the figure below. Obviously, the distributions are discriminative fro these three different type of pixels (on the edge, on flat area and a corner).(II)

For corners, the ranges in both $I_x$ and $I_y$ are large.
For points on an edge, one derivative has wide distribution while the other is almost all near zeros.
For points of a flat area, both derivatives are around zeros.

This observation is indeed the reason behind (I) we discussed above.

In Relationship between eigenvectors and directions of data distribution, we have shown that the eigenvectors of a covariance matrix are the same as the directions of principle component (i.e. the axes of a ellipsoid enclosing the data points), and the eigenvalues measure the variances of date points along eigenvectors. (Click Read more to see the relationship between eigenvalues and the semi-axes) Therefore, if the smallest eigenvalue is big enough, then it means the data have large variance along all different directions–it’s a corner!

In a word, the structure tensor is actually the covariance matrix of gradients for pixels around the pixel investigated, the the reason for current eigenvalue-based analysis is based on the property of covariance matrices.

[OPENCV] Feature detectors and descriptors

February 16, 2013 gccheng 2 comments

For a simple example of how to setup and use them, see [OPENCV] OpenCV+Eclipse|Qt|VS2010

For many computer vision tasks, the first step is usually extracting features from images. It consists of two parts: feature detection and description. By feature detection, we mean to find the pixels (or regions, objects) of interest in an image, which are represented by cv::KeyPoint. They are corresponding to the corners of objects in the image. OpenCV implements many popular feature points detection methods, such as SIFT, SURF, ORB, GoodFeatureToTrack and so on.

After these feature points are located, the next step is to describe them. An intuitive description of a feature point could be its intensity/color value. But a single intensity/color value does not provide sufficiently discriminative information so that this feature point can be uniquely matched to feature points in another image. Usually we use the information around a feature point to represent it, i.e. we use information of many pixels to represent one feature point. Based on Hubel & Wiesel’s research on information processing, orientation is the most sensitive factor to visual system. That’s why most feature techniques are based on orientation to describe feature points. For example, SIFT descriptor uses the normalized orientation histogram to describe a feature point. See figure below. Since this post is not about how the descriptors are designed, interested readers should read their papers.

To provide consistent programming interface, OpenCV designs the feature detectors and descriptor extractors following OO principle. Below is the class hierarchy of most detectors and descriptors for version 2.4.1. FeatureDetector and DescriptorExtractor are the base class for detectors and descriptors, respectively. Feature2d simplifies detection and extraction by deriving from FeatureDetector and DescriptorExtractor. But Feature2d is still an abstract class because it doesn’t implement the virtual function detectImpl() and computeImpl().

Classes SURF, SIFT, BRISK and ORB (and others) are concrete classes which are both detectors and descriptor extractors. In order to explicitly use them as detectors or descriptor extractors, aliases are defined. For example,

typedef SURF SurfFeatureDetector;
typedef SURF SurfDescriptorExtractor;

typedef SIFT SiftFeatureDetector;
typedef SIFT SiftDescriptorExtractor;

typedef ORB OrbFeatureDetector;
typedef ORB OrbDescriptorExtractor;

From the perspective of functionality, there’s no distinction between SurfFeatureDetector and SurfDescriptorExtractor. They are exactly the same thing except they have different names. The same is true for SiftFeatureDetecor and SiftDescriptorExtractor, and others.

From this class hierarchy diagram, it can also be seen that FAST, GFTTDetector and DenseFeatureDetector are only feature point detectors. After obtaining the feature points (std::vector<cv::KeyPoint>), we can use SURF, SIFT or any other descriptor extractors to get a description for these feature points. std::vector<cv::KeyPoint> is the data structure for data passing.

Categories: 计算机视觉, 杂记 Tags: opencv, software

Relationship between eigenvectors and directions of data distribution

February 7, 2013 gccheng 1 comment

For approaches based on dimension reduction, such as Principle Component Analysis, one critical question is “To which direction should we project the data?” Given $N$ observations of variable $X$ , we might compute the covariance matrix, and then use its eigenvectors as the direction of projection. Lots of approaches use this technique (e.g. Harris corner detector), but WHY? What’s the reason behind? This posts explains.

In the figure above, we illustrate a distribution of a 2D dataset. It’s already demeaned. The discussion below applies to any dimension size. Suppose each data point is $N$ -dimensional $X$ . Intuitively, we want to find linear functions of $X$ , which project $X$ into a new space. Since the data is already demeaned, $\bar{X}=0$ which means the origin of the axes is now at the center of the data, and the linear function has no intercept part.

We start by defining the first linear function (or projection): $\alpha_1^TX$ , where $\alpha_1$ is the coefficients of the linear function (see the figure above for 2D projection where $\alpha_1=(a,b)^T)$ . Looking at the data plot above, the line of $\alpha_1$ should pass through the data points, and the projected points should have largest span. Mathematically speaking, we want to maximize the data variance after projecting them to a line $\alpha$ , i.e. to find $\alpha_1$ so that

$\alpha_1 = \arg\max_\alpha \{Var(\alpha^TX)\}$

Example: For 2D data, $X=(X_1, X_2)$ , $\alpha=(a,b)$ . $\alpha^TX=aX_1+bX_2$ . It’s easy to show that the expectation $E(\alpha^TX)=aE(X_1)+bE(X_2)=a\mu_1+b\mu_2$ . The variance

$\begin{array} {rl} Var(\alpha^TX) & =E[ (aX_1+bX_2)-(a\mu_1+b\mu_2) ]^2\\ & =E\left[a(X_1-\mu_1)+b(X_2-\mu_2)\right]^2 \\ & =E\left[a^2(X_1-\mu_1)^2+b^2(X_2-\mu_2)^2+2ab(X_1-\mu_1)(X_2-\mu_2)\right] \\ & =a^2Var(X_1)+b^2Var(X_2)+2abCov(X_1, X_2) \\ & =a^2\sigma_{11}+b^2\sigma_{22}+2ab\sigma_{12} \\ & =(a,b)\begin{pmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{12} & \sigma_{22} \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix}\\ & = \alpha^T\Sigma\alpha \end{array}$

where $\Sigma$ is the covariance matrix of the original data. It is worth noticing that the finding $Var(\alpha^TX)=\alpha^T\Sigma\alpha$ can be extended to any higher dimensional variable $X$ . So right now our aim becomes

$\alpha_1 = \arg\max_\alpha \alpha^T\Sigma\alpha$

Without any constraint, we could always pick a larger $\alpha$ so there’s no solution to $\alpha_1$ . We choose normalized $\alpha$ , which means $\|\alpha\|^2=\alpha^T\alpha=1$ .

Solution: Now our problem becomes a optimization problem subject to a constraint. We use Lagrange multipliers to maximize the function

$\alpha^T\Sigma\alpha-\lambda(\alpha^T\alpha-1)$

with respect to $\alpha$ . This results in

$\begin{array}{rl} \frac{d}{d\alpha}\left(\alpha\Sigma\alpha-\lambda(\alpha^T\alpha-1)\right) & =0\\ \Sigma\alpha-\lambda\alpha & = 0 \\ \Sigma\alpha & = \lambda\alpha\\ \end{array}$

Suppose the solution to this equation is $\alpha=\alpha^*$ and $\lambda=\lambda_1$ .

It can be easily seen that the $\alpha_1 = \alpha^*$ is the eigenvector of $\Sigma$ . When the line (represented by $\alpha$ ) is in the same direction as the eigenvector associated with largest eigenvalue, the variance will be the largest when we project the original data onto this line, and the variable is:

$\alpha_1^T\Sigma\alpha_1 = \alpha_1^T\lambda\alpha_1 = \lambda_1\alpha_1^T\alpha_1 = \lambda_1$

Therefore, we should choose the eigenvector of the covariance matrix with largest eigenvalue as the basis of first dimension (first principle component), and the corresponding eigenvalue is indeed the variance of that projected data in that dimension.

__________________

Following the same deduction, we can conclude that the $k$ -th eigenvector $\alpha_k$ is the $k$-th principle component, and $Var(\alpha^TX)=\lambda_k$ . For example, for $k=1$ ,

$\begin{array}{rl} &\arg\max_{\alpha_2} \alpha_2^T\Sigma\alpha_2 \\ s.t. &\alpha_2^T\alpha_2 = 1 \\ &\alpha_2\alpha_1=0 \end{array}$

(The second due to $cov(\alpha_1^TX,\alpha_2^TX)=\alpha_1^T\Sigma\alpha_2=\alpha_2^T\Sigma\alpha_1=\alpha_2^T\lambda_1\alpha_1=\lambda_1\alpha_2\alpha_1=0$ .)

Use Lagrange multiplier again, you will get

$\Sigma\alpha_2 = \lambda_2\alpha_2$

___________________

If you read here, I hopeit is clear now the reason why the eigenvector with largest eigenvalue coincide with the direction of biggest variance and why the eigenvectors of a covariance matrix are used as the directions of principle components.

Note: One more question in my mind is why we are so furtunate to have such a simple yet powerful solution, but I want to stop here today.

Categories: 统计, 学习&研究 Tags: statistics, math, PCA

Statistical Data Mining Basics

November 14, 2012 gccheng Leave a comment

Statistical Data Mining Tutorials
by Andrew Moore

The following links point to a set of tutorials on many aspects of statistical data mining, including the foundations of probability, the foundations of statistical data analysis, and most of the classic machine learning and data mining algorithms.

Articles on Doctoral Education

November 8, 2012 gccheng Leave a comment

A deeper look into the “80% of PhDs who do not become professors”

Doctorate recipients holding tenure and tenure-track appointments at academic institutions

Improving Graduate Education to Support a Branching Career Pipeline

Say NO to the Second Post Doc!

Categories: 杂记

[OPENGL]Some OpenGL Examples

August 30, 2012 gccheng Leave a comment

Here are some OpenGL programs written for course assignments, from the very basic one to some advanced application. The codes are hosted on github：https://github.com/gccheng/openGL-tutorials

1. Rotating pendulum

2. Hand & fingers

3. Solar System

基于线框模型的简单物体绘制与操作。

4. Array operation (Rotation)
绘制一个函数的图像，基于Array进行绘制。下面的图形是函数 $z=f(x,y)=\frac{1}{2}e^{-0.4\sqrt{(80x-40)^2+(90y-45)^2}}cos(0.15\sqrt{(80x-40)^2+(90y-45)^2})$ 实时绘制的结果。

function from gc-cheng@hotmail.com on Vimeo.

5. Quaternion algebra
使用quaternion表示旋转、缩放等操作。由于使用了Quaternion，极大地简化了操作和代码：

planet-quaternion algebra from gc-cheng@hotmail.com on Vimeo.

6. 分段贝塞尔(Bezier)曲线

一段Bezier曲线是由四个点确定的，两段Bezier曲线是将7个点组成的一条连续可导的Bezier曲线，的原理可以查询Wikipedia中的文章。

实现原理：Piecewise Bezier Curve

实现代码：prog6_cheng.cpp

实现效果：

Categories: 计算机图形学, 杂记 Tags: graphics, OpenGL

Last week I attended the commencement ceremony, and I really found most of the graduates were females.

August 16, 2012 gccheng Leave a comment

http://usat.ly/MBiBha Women’s colleges struggle to keep enrollment As enrollment numbers drop, women’s colleges are forced to make decisions about their admission policies. To view this story, click the link or paste it into your browser.

Categories: 杂记

温故知新：CXX/JAVA

July 21, 2012 gccheng Leave a comment

上周花了点时间把《The C++ Standard Library: A Totorial and Reference》浏览了一下，有段时间没有写长的C++程序了，温习知新。天天捣腾那么点Matlab程序，要不是CV领域最强大的程序库是C++写的，不然编程能力真的要倒退好几年了。会继续在闲暇之余温习温习基本的数据结构和算法。

Project STL: http://compilr.com/gccheng/stl

1. 使用迭代器适配器（Iterator Adapter)进行容器拷贝: container_copy

2. 使用流迭代器: stream_iterator

3. 函数对象(Function Object)如何实现和使用: my_foreach

4. C++11中新增的array容器，提供了cbegin(), cend()，没深研究如何以前的容器没提供: it_traverse

Project JavaAlgorithm: http://compilr.com/gccheng/javaalgorithm

1. 树的四种遍历: TreeTraversal

2. 统计词频: WordCount

Categories: Programming Languages Tags: C++, java, STL

Office VBA二次开发

July 18, 2012 gccheng Leave a comment

看到两篇文章：《为什么我们不重视Office开发——与郭安定谈话想到的》和《也谈为什么我们不重视Office开发－－看孟岩说有感》，结论似乎是基于Office的二次开发主要适用于一些内部的小项目，对于商业化开发，似乎就不太合适。

联想到以前自己大概做过三四个Excel的VBA开发，使用自动化以后，确实提高了工作效率，这包括以前给老婆单位做的一个通过，手工统计每个员工的工时，这可能需要很多时间：每个人每天上几个小时（每天可能不同，小时制），一个月来了几天，几百号人。第二个是去年帮一个同学弄的基于Excel的员工自动签到。第三个是一个从该表格中提取数据，按照规定格式，自动打印的功能。

对于第三个，小册子是由上级部门打印好分发下来的，数据则是每个乡镇自己统计的Excel表格数据，现在要做的工作就是把表格里的每一行使用打印机打印到预制好的小册子上。由于工作是在基层乡镇，打印室一半可能就是用Word或Excel设置要每一行数据在纸上的位置（模板），然后把每一行中的各列数据添加到这个Word或Excel中，然后打印到小册子上。这种工作对每一行都要重复一下。如果有几万人，就算每个用2分钟，这个工作量也是很大的。而且小册子上的有些位置的数据可能不是直接从数据表中得到的（例如是多项的求和等）。

所以当时就根据已有的经验，设置了一个Excel表格数据格式化自动打印的程序，其实程序非常简单：只要打印室的人员根据小册子的要求设计好模板（这个是一次性工作，而且是他们非常擅长的），然后按一下“打印”按钮，剩下的工作就是等待打印完成了。

就是这个简单的程序，节省了大量的时间和精力，而且这种工作每年都有，每年都可以用，即使有小改动，也是很轻易完成的。

从这里可以获得一份示意Excel。

Categories: Programming Languages Tags: software, vba

如何从任意分布采样How to sample from arbitrary standard distributions

July 14, 2012 gccheng 1 comment

The topic of this article is about how to produce data following the specified simple distribution. For example, the question we want to answer is like “can you generate 1000 data points from the exponential distribution?”.

Prior to the discussion of today’s topic, let’s first look at the random number generation methods available in prevalent programming languages. Almost every popular languange provides the pseudo-random number generation, most of which have methods to produce data from uniform distribution. For example, to generate an integral number uniformly from {0, 1, 2, …, 9}, we have

/* C style */
int randInt = random (10);
// C++ style
srand((unsigned)time(0));
int randInt = rand() % 10;
// Java style
Random randomGenerator = new Random();
int randInt = randomGenerator.nextInt(10);

Why do we need random numbers in these languages? Part of the reasons is random numbers can be used to simulate different scenaria in which realistic data are hard to obtain. For example, in Mahattan project, it’s almost impossible to use real nuclear explosion data to verify the system so that Markov chain Monte Carlo sampling method was proposed to generate simulation data. For software test, we also need looooooots of different inputs to test the target software, and the random inputs can be used based on random numbers.

Why is uniform distribtion usually implmented? I guess it is because of its simplicity: itself (the distributoin and its implementation) is simple, and it’s relatively simple to generate other distribution using uniform distribution. Interested reader can learn how the implementation are implemented in C/C++/Java.

Here, we will learn how to transform uniform distribution to other simple distributions. (Actullay the idea can be applied to the transformation from one distribution to another) These distribution are “simple” in terms that the inverse of its cumulative distribution function (CDF) can be easily obtained. You will see this later. Suppose we’ve already known how to sample from uniform distribution (*), i.e. we can generate data $z$ over the interval (0,1) from a uniform distribution. Now we want to sample $y$ from the probability distribution of $p(y)$ (f(y) for continuous variable and P(y) for discrete).

Let’s look at a discrete example. Let $y\in \{1,2,3,4\}$ , and $P(y=1)=0.3$ , $P(y=2)=0.4$ , $P(y=3)=0.2$ , $P(y=4)=0.1$ . We want to generate 100 values from $\{1,2,3,4\}$ so that these values have the same (similar) distribution as $y$ . Firstly, the interval (0,1) is divided into 4 intervals [F(0)=0, F(1)], [F(1), F(2)], [F(2), F(3)], [F(3), F(4)] as shown in the figure below where $F(i)=\Sigma_{j=1}^i{P(y=j)}$ :

The width of interval $i$ is $P(y=i)$ . We then draw one point $Z=z$ from (0, 1) uniformly (since we know how to draw samples from uniform distribution (*) ), and search the index $i$ using $F(i-1) \leq z < F(i)$ and assign $y=i$ which has the probability $P(i)$ . Repeat this procedure for 100 times, we can get 100 $y$ values, which have the distribution $P(Y=y)$ .

Think more about this simple example and we can find the procedure is: (1) sample from uniform distribution to get a value $z$ ; (2) solve $y$ from $z$ using the probability distribution funciton $F(y)$ : $z=F(y)$ or its generalized form $Inf_z\{z\leq F(y)\}$ . In this way, we can sample from any discrete distributions whose probability desitribution functions are known.

Extend this to continuous distributions, this basic idea is similar. Suppose we want to sample from distribution $f(y)$ , and its cumulative distribution function (or distribution function) is $h(y)$ as shown in the following figure. Here the vertical axis is three-in-one: $p(y)$ , $h(y)$ and $z$ . Because the range of $h(y)$ is always in (0,1), we can treat it as an interval from which we uniformly draw samples $z$ . Following the same procedure as for discrete case, we draw a point $z$ over (0,1) uniformly, and treat it as the value of $h(y)$ for some unknown $y$ (analogous to $F(i)$ in discrete case), i.e.

$z=h(y)\equiv \int_{-\infty}^y{p(\hat{y})d\hat{y}}$ (1)

Now the problem is “can we get a unique $y$ from $z=h(y)$ ?” or “is $h(y)$ inversible?”. Simply put, “Maybe. There might be several $y$ ‘s which satisfy the equation”.

We know that the CDF ( $h(y)$ ) of a distribution is monotone non-decreasing. Because $p(x)\geq 0$ , it is easy to see from equation (1) that $h(y)$ cannot decrease as $y$ increases.

If CDF $h(y)$ is strict increasing (or $\forall x p(x)>0$ ), then $z=h(y)$ is inversible, and we can get $y=h^{-1}(z)$ ;
If CDF $h(y)$ is not strict increasing (or $\exists x p(x)=0$ ), then $z=h(y)$ is not inversible. One may define a generalized inverse distribution function: $h^{-1}(z)=\inf_{y \in \mathbb{R}} \{ h(y) \geq z \}$ , and then assign it to $y$ .

Repeat this process, we can get as many $y$ values as we want.

One can verify that the probability distribution function of $y$ is $h(y)$ and thus the probability distribution is $p(y)$ . We transform the value of $z$ to $y$ using equation (1), and thus the distribution function of $y$ will be

$P(Y \leq y)=P(h^{-1}(Z) \leq y)=P(Z \leq h(y))=h(y)$ (2)

where the last equation is because $Z$ is a uniform distribution over (0,1), which means $P(Z \leq \alpha)=\alpha (\alpha \in [0,1])$

One can also “feel” from the curve of $p(y)$ and $h(y)$ . For $y$ ‘s whose probability $p(y)$ are big, the corresponding $h(y)$ is steep, which results in large occupation of interval (0, 1) on the vertical axis. Therefore, with large chance the $z$ values will fall in this occupation when we draw samples from (0,1) uniformly, and thus we have larger probability to obtain $y$ s whose probabilities $p(y)$ are large. Readers can try to associate this idea with the discrete example above.

Thus far, we can see the core of the sampling from arbitrary distribution is how to transform from one distribution to the other. The article only talks about “simple” distributions, which mean it’s easy to get the CDF $h(x)$ and its inverse, or it’s impossble to apply the method above.

***********************************************************************************

To answer the question we asked at the very beginning, we want to sample from an exponential distribtuion

$p(y) = \lambda e^{-\lambda y}$ (3)

where $0 \leq y < \infty$ . In this case the lower limit of the integreal in (1) is 0, and so $h(y)=1-e^{-\lambda y}$ . Thus, if we transform our uniformly distributed varable $z$ using $y=-\lambda^{-1}ln(1-z)$ , then $y$ will have an exponential distribution.

***

Image histogram equalization is another application of the distribution transformation. One difference of histogram equalization compared with the discussion above is it tries its best to transform a non-uniform distribution (histogram) to a near-uniform distribution (histogram) while keeping the relative order of pixels at different scales. Therefore,

Compute histogram (probability distribution) of pixels in original image $P_r(z)$
Set $y=\Sigma_{z=0}^{y}P_r(z)$
Convert the range of $y$ from (0, 1) to (0, 255). This step is because it’s the target domain is not within (0,1).

Refer to Histogram Equalization in wikipedia and Gonzalez’s DIP for more information.

One limitation of the method discussed above is the computation of inverse of CDF, which in many cases doesn’t have an easy form or easy to compute. We will see how to solve this problem in next blog.

Categories: 统计 Tags: sampling, statistics

Older Entries