利用统计方法求π（PI），并可视化显示求解过程（C++&&OpenCV）

  统计方法求π的方式：如果在正方形区域内随机产生大量的均匀分布的点，那么落入内切圆和正方形中的随机点个数的比值等于它们的面积之比。该比值乘以4，即为PI值。这就是统计方法求π的过程。

  可视化求解过程是指：把产生随机点的过程在图像中显示,可视化过程用到了opencv库，利用opencv创建一个图像，并将数据写入到图像数据的内存位置。

具体要求：

1、显示一幅边长为R的正方形（设置为0，全黑色），利用已经给出的画圆函数，显示出此正方形边长为R/2的内切圆（设置为255，白色），圆心位于正方形中心；

2、通过系统函数rand（）产生随机数，将每一个产生的点写入图像的数组，利用上次实验的函数写入图像并在上述的正方形显示区域内显示出来，赋值为255（白色）。每个点的亮度随着时间（新的随机点的产生）减弱至0。

现在我们分步处理：
（一）：显示一幅边长为R的正方形，（设置为0，全黑色）：

先创建一个图片（将P设置为500）：

int width = 500;
int height = 500;
IplImage* img = cvCreateImage(cvSize(width, height), IPL_DEPTH_8U, 1);
unsigned char *img_data = (unsigned char *)(img->imageData);
int width_step = img->widthStep;//width_step在灰度图像中表示一行数据的个数

img_data 是这个图像灰度值存贮的内存首地址。

动态声明一个二维数组，用来存储图像灰度值。

unsigned char **array_of_img; //开辟一个最初始的地址
array_of_img = new unsigned char *[width_step];//开辟的行数
array_of_img[0] = new unsigned char[width_step * width_step];//开辟的总空间

for (int i = 1; i < width_step; ++i)//将每一行的行首地址赋值
{
array_of_img[i] = array_of_img[i-1] + width_step;
}

然后将动态声明的二维数组初始化为0，为下面将图像灰度值赋0做准备；

for (int k = 0; k < width_step; ++k)
{
for (int j = 0; j < width_step; ++j)
{
array_of_img[k][j] = 0;
}
}

调用WriteImageData()函数，将二维数组中的灰度值写入到图像数据中：

void WriteImageData(unsigned char *src, int rows, int cols, int width_step,unsigned char ** &dst, bool flag)
{
//src图像数据在内存的首地址  rows 行，cols列，dst二维数组,flag 是否赋值的参数
int count =0;
for (int i = 0; i < cols; ++i)
{
for (int j = 0; j < rows; ++j)
{
if(flag) src[count] = dst[i][j];
else src[count] = (unsigned char)'a';
count++;
}
}
printf("%d\n", count);
return;
}

(二)利用已经给出的画圆函数，显示出此正方形边长为R/2的内切圆（设置为255，白色），圆心位于正方形中心：

调用opencv库的画圆函数，并显示图像：

cvCircle(img, cvPoint(250,250), 250, cvScalar(255), 1, 8, 0);
cvNamedWindow("Image", CV_WINDOW_AUTOSIZE);
cvShowImage("Image", img);
cvWaitKey(0);

(三)通过系统函数rand（）产生随机数，将每一个产生的点写入图像的数组，利用上次实验的函数写入图像并在上述的正方形显示区域内显示出来，赋值为255（白色）。每个点的亮度随着时间（新的随机点的产生）减弱至0；

int sumSquare = 10000;
int sumCircle = 0;

double PI = 0.0;
//generate 100 point(x,y)
for(int count = 0; count < sumSquare; ++count)
{
//数组中所有非零元素的值减5 使颜色变暗
for (int n = 0; n < width_step; ++n)
{
for(int m = 0; m < width_step; ++m )
{
if (array_of_img[m]
> 0)
{
array_of_img[m]
-= 5;
}
}
}
//产生一个(0,1)之间的 随机点
double numA = (double)(rand())/(double)(RAND_MAX);
double numB = (double)(rand())/(double)(RAND_MAX);

//将[0,1]的点 放大到 [0,500] 之间
int Px = (int) (numA * 500);
int Py = (int) (numB * 500);
//标示 刚生成的点
array_of_img[Px][Py] = 255;
// 把数组写入图像中
WriteImageData(img_data, width_step, width_step,width_step, array_of_img, true);
//重新绘制内切圆
cvCircle(img, cvPoint(250,250), 250, cvScalar(255), 1, 8, 0);

//判断刚生成的点 是否在 园内
if(((numA - 0.5)*(numA - 0.5) + (numB - 0.5)*(numB - 0.5)) <= 0.25)
{
sumCircle++;
}
cvShowImage("Image", img);//显示图像
cvWaitKey(50); //延时50毫秒

PI =(double)sumCircle / (double)count;
PI *= 4;
printf("%d %d\n",sumCircle,count);
printf("The PI is : %lf\n",PI);
}

这个可视化过程其实就是在求π的循环内对图像数据进行赋值的过程，一个write函数就可。

下面是完整的程序源代码：

#include "cv.h"
#include "highgui.h"
#include <cstdio>

//read the img data from the Memory to own array
void ReadImageData(unsigned char *src, int rows, int cols, int width_step, unsigned char ** &dst)
{
int count = 0;
for (int i = 0; i < cols; ++i)
{
for (int j = 0; j < rows; ++j)
{
// printf("%c\n", *((unsigned char*)dst + rows * i + j));
// *((unsigned char*)dst + rows * i + j) = src[count];
dst[i][j] = src[count];
count++;
}
}
printf("%d\n", count);
return;
}
//将用户自己数组存储的输入拷贝到图像相应的内存位置
//write the img data ( my array) to the memory
void WriteImageData(unsigned char *src, int rows, int cols, int width_step,unsigned char ** &dst, bool flag)
{
int count =0;
for (int i = 0; i < cols; ++i)
{
for (int j = 0; j < rows; ++j)
{
if(flag) src[count] = dst[i][j];
else src[count] = (unsigned char)'a';
count++;
}
}
printf("%d\n", count);

return;
}

int main(int args,char * argv[])
{
//	double PI = 0.0;
//设定图片的大小
int width = 500;
int height = 500;
//利用opencv创建一个 空的图片图像
IplImage* img = cvCreateImage(cvSize(width, height), IPL_DEPTH_8U, 1);
unsigned char *img_data = (unsigned char *)(img->imageData);
//获得 一行 大小
int width_step = img->widthStep;

//动态创建 正方形的 二维数组  边长为width
unsigned char **array_of_img;
array_of_img = new unsigned char *[width_step];
array_of_img[0] = new unsigned char[width_step * width_step];

//为每一行 分配地址
for (int i = 1; i < width_step; ++i)
{
array_of_img[i] = array_of_img[i-1] + width_step;
}
//初始化动态 声明的二维数组 0，即黑色
for (int k = 0; k < width_step; ++k)
{
for (int j = 0; j < width_step; ++j)
{
array_of_img[k][j] = 0;
}
}
//将数组的数据写入到内存相应的位置
WriteImageData(img_data, width_step, width_step,width_step, array_of_img, true);
//画圆
cvCircle(img, cvPoint(250,250), 250, cvScalar(255), 1, 8, 0);

cvNamedWindow("Image", CV_WINDOW_AUTOSIZE);//创建窗口
cvShowImage("Image", img);//显示图像
cvWaitKey(0); //等待按键，开始采样随机点

int sumSquare = 10000;
int sumCircle = 0;

double PI = 0.0;
//generate 100 point(x,y)
for(int count = 0; count < sumSquare; ++count)
{
//数组中所有非零元素的值减5 使颜色变暗
for (int n = 0; n < width_step; ++n)
{
for(int m = 0; m < width_step; ++m )
{
if (array_of_img[m]
> 0)
{
array_of_img[m]
-= 5;
}
}
}
//产生一个(0,1)之间的 随机点
double numA = (double)(rand())/(double)(RAND_MAX);
double numB = (double)(rand())/(double)(RAND_MAX);

//将[0,1]的点 放大到 [0,500] 之间
int Px = (int) (numA * 500);
int Py = (int) (numB * 500);
//标示 刚生成的点
array_of_img[Px][Py] = 255;
// 把数组写入图像中
WriteImageData(img_data, width_step, width_step,width_step, array_of_img, true);
//重新绘制内切圆
cvCircle(img, cvPoint(250,250), 250, cvScalar(255), 1, 8, 0);

//判断刚生成的点 是否在 园内
if(((numA - 0.5)*(numA - 0.5) + (numB - 0.5)*(numB - 0.5)) <= 0.25)
{
sumCircle++;
}
cvShowImage("Image", img);//显示图像
cvWaitKey(50); //延时50毫秒

PI =(double)sumCircle / (double)count;
PI *= 4;
printf("%d %d\n",sumCircle,count);
printf("The PI is : %lf\n",PI);
}
PI =(double)sumCircle / (double)sumSquare;
PI *= 4;
printf("The PI is : %lf\n",PI);
cvDestroyWindow("Image");//销毁窗口
cvReleaseImage(&img); //释放图像

delete [] array_of_img[0];
delete [] array_of_img;

return 0;
}