开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
SegApi.cpp
资源名称：图像识别.rar　[点击查看]
上传用户：xjt2008yy
上传日期：2010-01-18
资源大小：272k
文件大小：42k
源码类别：

生物技术

开发平台：

Visual C++

SegApi.cpp：源码内容
  1. #include "GlobalApi.h"
  2. #include "stdafx.h"
  3. #include "cdib.h"
  4. #include "math.h"
  5. #include <direct.h>
  6. #include <complex>
  7. using namespace std;
  9. /*************************************************************************
  10.  *
  11.  * 函数名称：
  12.  *   RegionSegFixThreshold()
  13.  *
  14.  * 输入参数:
  15.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  16.  *   int nThreshold - 区域分割的阈值
  17.  *
  18.  * 返回值:
  19.  *   无
  20.  *
  21.  * 说明:
  22.  *   1（逻辑）表示对应象素为前景区域，0表示背景
  23.  *   阈值分割的关键问题在于阈值的选取。阈值的选取一般应该视实际的应用而
  24.  *   灵活设定。
  25.  *
  26.  *************************************************************************
  27.  */
  28. void RegionSegFixThreshold(CDib * pDib, int nThreshold)
  29. {
  30. //遍历图象的纵坐标
  31. int y;
  33. //遍历图象的横坐标
  34. int x;
  36. //图象的长宽大小
  37. CSize sizeImage = pDib->GetDimensions();
  38. int nWidth = sizeImage.cx ;
  39. int nHeight = sizeImage.cy ;
  41. //图像在计算机在存储中的实际大小
  42. CSize sizeImageSave = pDib->GetDibSaveDim();
  44. //图像在内存中每一行象素占用的实际空间
  45. int nSaveWidth = sizeImageSave.cx;
  48. //图像数据的指针
  49. LPBYTE  pImageData = pDib->m_lpImage;
  51. for(y=0; y<nHeight ; y++ )
  52. for(x=0; x<nWidth ; x++ )
  53. {
  54. if( *(pImageData+y*nSaveWidth+x) < nThreshold)
  55. *(pImageData+y*nSaveWidth+x) = 0;
  56. }
  57. }
  59. /*************************************************************************
  60.  *
  61.  * 函数名称：
  62.  *   RegionSegAdaptive()
  63.  *
  64.  * 输入参数:
  65.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  66.  *
  67.  * 返回值:
  68.  *   无
  69.  *
  70.  * 说明:
  71.  *   1（逻辑）表示对应象素为前景区域，0表示背景
  72.  *   阈值分割的关键问题在于阈值的选取。阈值的选取一般应该视实际的应用而
  73.  *   灵活设定。本函数中，阈值不是固定的，而是根据图象象素的实际性质而设定的。
  74.  *   这个函数把图像分成四个子图象，然后计算每个子图象的均值，根据均值设置阈值
  75.  *   阈值只是应用在对应的子图象
  76.  *
  77.  *************************************************************************
  78.  */
  79. void RegionSegAdaptive(CDib * pDib)
  80. {
  81. //遍历图象的纵坐标
  82. int y;
  84. //遍历图象的横坐标
  85. int x;
  87. //图象的长宽大小
  88. CSize sizeImage = pDib->GetDimensions();
  89. int nWidth = sizeImage.cx ;
  90. int nHeight = sizeImage.cy ;
  92. //图像在计算机在存储中的实际大小
  93. CSize sizeImageSave = pDib->GetDibSaveDim();
  95. //图像在内存中每一行象素占用的实际空间
  96. int nSaveWidth = sizeImageSave.cx;
  98. //图像数据的指针
  99. LPBYTE  lpImage = pDib->m_lpImage;
  101. // 局部阈值
  102. int nThd[2][2] ;
  104. // 子图象的平均值
  105. int nLocAvg ;
  107. // 对左上图像逐点扫描：
  108. nLocAvg = 0 ;
  109. // y方向
  110. for(y=0; y<nHeight/2 ; y++ )
  111. {
  112. // x方向
  113. for(x=0; x<nWidth/2 ; x++ )
  114. {
  115. nLocAvg += lpImage[y*nSaveWidth + x];
  116. }
  117. }
  118. // 计算均值
  119. nLocAvg /= ( (nHeight/2) * (nWidth/2) ) ;
  121. // 设置阈值为子图象的平均值
  122. nThd[0][0] = nLocAvg ;
  124. // 对左上图像逐点扫描进行分割：
  125. // y方向
  126. for(y=0; y<nHeight/2 ; y++ )
  127. {
  128. // x方向
  129. for(x=0; x<nWidth/2 ; x++ )
  130. {
  131. if(lpImage[y*nSaveWidth + x]<nThd[0][0])
  132. lpImage[y*nSaveWidth + x] = 255 ;
  133. else
  134. {
  135. lpImage[y*nSaveWidth + x] = 0 ;
  136. }
  138. }
  139. }
  141. // =============================================
  142. // 对左下图像逐点扫描：
  143. nLocAvg = 0 ;
  144. // y方向
  145. for(y=nHeight/2; y<nHeight ; y++ )
  146. {
  147. // x方向
  148. for(x=0; x<nWidth/2 ; x++ )
  149. {
  150. nLocAvg += lpImage[y*nSaveWidth + x];
  151. }
  152. }
  153. // 计算均值
  154. nLocAvg /= ( (nHeight - nHeight/2) * (nWidth/2) ) ;
  156. // 设置阈值为子图象的平均值
  157. nThd[1][0] = nLocAvg ;
  159. // 对左下图像逐点扫描进行分割：
  160. // y方向
  161. for(y=nHeight/2; y<nHeight ; y++ )
  162. {
  163. // x方向
  164. for(x=0; x<nWidth/2 ; x++ )
  165. {
  166. if(lpImage[y*nSaveWidth + x]<nThd[1][0])
  167. lpImage[y*nSaveWidth + x] = 255 ;
  168. else
  169. {
  170. lpImage[y*nSaveWidth + x] = 0 ;
  171. }
  173. }
  174. }
  176. // =============================================
  177. // 对右上图像逐点扫描：
  178. nLocAvg = 0 ;
  179. // y方向
  180. for(y=0; y<nHeight/2 ; y++ )
  181. {
  182. // x方向
  183. for(x=nWidth/2; x<nWidth ; x++ )
  184. {
  185. nLocAvg += lpImage[y*nSaveWidth + x];
  186. }
  187. }
  188. // 计算均值
  189. nLocAvg /= ( (nHeight/2) * (nWidth - nWidth/2) ) ;
  191. // 设置阈值为子图象的平均值
  192. nThd[0][1] = nLocAvg ;
  194. // 对右上图像逐点扫描进行分割：
  195. // y方向
  196. for(y=0; y<nHeight/2 ; y++ )
  197. {
  198. // x方向
  199. for(x=nWidth/2; x<nWidth ; x++ )
  200. {
  201. if(lpImage[y*nSaveWidth + x]<nThd[0][1])
  202. lpImage[y*nSaveWidth + x] = 255 ;
  203. else
  204. {
  205. lpImage[y*nSaveWidth + x] = 0 ;
  206. }
  208. }
  209. }
  211. // =============================================
  212. // 对右下图像逐点扫描：
  213. nLocAvg = 0 ;
  214. // y方向
  215. for(y=nHeight/2; y<nHeight ; y++ )
  216. {
  217. // x方向
  218. for(x=nWidth/2; x<nWidth ; x++ )
  219. {
  220. nLocAvg += lpImage[y*nSaveWidth + x];
  221. }
  222. }
  223. // 计算均值
  224. nLocAvg /= ( (nHeight - nHeight/2) * (nWidth - nWidth/2) ) ;
  226. // 设置阈值为子图象的平均值
  227. nThd[1][1] = nLocAvg ;
  229. // 对右下图像逐点扫描进行分割：
  230. // y方向
  231. for(y=nHeight/2; y<nHeight ; y++ )
  232. {
  233. // x方向
  234. for(x=nWidth/2; x<nWidth ; x++ )
  235. {
  236. if(lpImage[y*nSaveWidth + x]<nThd[1][1])
  237. lpImage[y*nSaveWidth + x] = 255 ;
  238. else
  239. {
  240. lpImage[y*nSaveWidth + x] = 0 ;
  241. }
  242. }
  243. }
  245. // 为了显示方便显示，逻辑1用黑色显示，逻辑0用白色显示
  246. for(y=0; y<nHeight ; y++ )
  247. {
  248. // x方向
  249. for(x=0; x<nWidth ; x++ )
  250. {
  251. lpImage[y*nSaveWidth + x] = 255 - lpImage[y*nSaveWidth + x] ;
  252. }
  253. }
  254. }
  256. /*************************************************************************
  257.  *
  258.  * 函数名称：
  259.  *   RobertsOperator()
  260.  *
  261.  * 输入参数:
  262.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  263.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  264.  *
  265.  * 返回值:
  266.  *   无
  267.  *
  268.  * 说明:
  269.  *   Roberts算子
  270.  *
  271.  *************************************************************************
  272.  */
  273. void RobertsOperator(CDib * pDib, double * pdGrad)
  274. {
  275. // 遍历图象的纵坐标
  276. int y;
  278. // 遍历图象的横坐标
  279. int x;
  281. // 图象的长宽大小
  282. CSize sizeImage = pDib->GetDimensions();
  283. int nWidth = sizeImage.cx ;
  284. int nHeight = sizeImage.cy ;
  286. // 图像在计算机在存储中的实际大小
  287. CSize sizeImageSave = pDib->GetDibSaveDim();
  289. // 图像在内存中每一行象素占用的实际空间
  290. int nSaveWidth = sizeImageSave.cx;
  292. // 图像数据的指针
  293. LPBYTE  pImageData = pDib->m_lpImage;
  295. // 初始化
  296. for(y=0; y<nHeight ; y++ )
  297. for(x=0 ; x<nWidth ; x++ )
  298. {
  299. *(pdGrad+y*nWidth+x)=0;
  300. }
  303. // 下面开始利用Roberts算子进行计算，为了保证计算所需要的
  304. // 的数据位于图像数据的内部，下面的两重循环的条件是
  305. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  306. // 而不是x<nWidth
  308. //这两个变量用来表示Roberts算子第一个模板的两个象素值
  309. int nUpLeft;
  310. int nDownRight;
  312. // 这两个变量用来表示Roberts算子第二个模板的两个象素值
  313. int nUpRight;
  314. int nDownLeft;
  316. // 这两个变量用来表示Roberts算子计算的结果
  317. int nValueOne;
  318. int nValueTwo;
  320. // 临时变量
  321. double dGrad;
  323. for(y=0; y<nHeight-1 ; y++ )
  324. for(x=0 ; x<nWidth-1 ; x++ )
  325. {
  326. // Roberts算子第一个模板需要的象素值
  327. nUpLeft =*(pImageData+y*nSaveWidth+x) ; 
  328. nDownRight =*( pImageData+(y+1)*nSaveWidth+x+1 );
  329. nDownRight *=-1;
  331. //Roberts算子的第一个模板计算结果
  332. nValueOne =nUpLeft+nDownRight ;
  334. // Roberts算子第二个模板需要的象素值
  335. nUpRight =*( pImageData+y*nSaveWidth+x+1 ) ;
  336. nDownLeft =*( pImageData+(y+1)*nSaveWidth+x );
  337. nDownLeft *=-1;
  339. // Roberts算子的第二个模板计算结果
  340. nValueTwo =nUpRight+nDownLeft;
  342. // 计算两个偏导数的平方和
  343. dGrad=nValueOne*nValueOne + nValueTwo*nValueTwo;
  345. // 开方
  346. dGrad=pow(dGrad,0.5);
  348. // 范数采用欧式距离
  349. *(pdGrad+y*nWidth+x)=dGrad;
  350. }
  351. }
  352. /*************************************************************************
  353.  *
  354.  * 函数名称：
  355.  *   LaplacianOperator()
  356.  *
  357.  * 输入参数:
  358.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  359.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  360.  *
  361.  * 返回值:
  362.  *   无
  363.  *
  364.  * 说明:
  365.  *   LaplacianOperator算子,是二阶算子，不想Roberts算子那样需要两个模板计算
  366.  *   梯度，LaplacianOperator算子只要一个算子就可以计算梯度。但是因为利用了
  367.  *   二阶信息，对噪声比较敏感
  368.  *
  369.  *************************************************************************
  370.  */
  371. void LaplacianOperator(CDib * pDib, double * pdGrad)
  372. {
  373. // 遍历图象的纵坐标
  374. int y;
  376. // 遍历图象的横坐标
  377. int x;
  379. // 图象的长宽大小
  380. CSize sizeImage = pDib->GetDimensions();
  381. int nWidth = sizeImage.cx ;
  382. int nHeight = sizeImage.cy ;
  384. // 图像在计算机在存储中的实际大小
  385. CSize sizeImageSave = pDib->GetDibSaveDim();
  387. // 图像在内存中每一行象素占用的实际空间
  388. int nSaveWidth = sizeImageSave.cx;
  390. // 图像数据的指针
  391. LPBYTE  lpImage = pDib->m_lpImage;
  393. // 初始化
  394. for(y=0; y<nHeight ; y++ )
  395. for(x=0 ; x<nWidth ; x++ )
  396. {
  397. *(pdGrad+y*nWidth+x)=0;
  398. }
  400. // 设置模板系数
  401. static int nWeight[3][3] ;
  402. nWeight[0][0] = -1 ;   
  403. nWeight[0][1] = -1 ;   
  404. nWeight[0][2] = -1 ;   
  405. nWeight[1][0] = -1 ;   
  406. nWeight[1][1] =  8 ;   
  407. nWeight[1][2] = -1 ;   
  408. nWeight[2][0] = -1 ;   
  409. nWeight[2][1] = -1 ;   
  410. nWeight[2][2] = -1 ;   
  414. //这个变量用来表示Laplacian算子象素值
  415. int nTmp[3][3];
  417. // 临时变量
  418. double dGrad;
  420. // 模板循环控制变量
  421. int yy ;
  422. int xx ;
  425. // 下面开始利用Laplacian算子进行计算，为了保证计算所需要的
  426. // 的数据位于图像数据的内部，下面的两重循环的条件是
  427. // y<nHeight-2 而不是y<nHeight，相应的x方向也是x<nWidth-2
  428. // 而不是x<nWidth
  429. for(y=1; y<nHeight-2 ; y++ )
  430. for(x=1 ; x<nWidth-2 ; x++ )
  431. {
  432. dGrad = 0 ; 
  433. // Laplacian算子需要的各点象素值
  435. // 模板第一行
  436. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  437. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  438. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  440. // 模板第二行
  441. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  442. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  443. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  445. // 模板第三行
  446. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  447. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  448. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  450. // 计算梯度
  451. for(yy=0; yy<3; yy++)
  452. for(xx=0; xx<3; xx++)
  453. {
  454. dGrad += nTmp[yy][xx] * nWeight[yy][xx] ;
  455. }
  457. // 梯度值写入内存
  458. *(pdGrad+y*nWidth+x)=dGrad;
  459. }
  461. }
  463. /*************************************************************************
  464.  *
  465.  * 函数名称：
  466.  *   RegionGrow()
  467.  *
  468.  * 输入参数:
  469.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  470.  *   unsigned char * pUnRegion - 指向区域生长结果的指针
  471.  *
  472.  * 返回值:
  473.  *   无
  474.  *
  475.  * 说明:
  476.  *   pUnRegion指针指向的数据区存储了区域生长的结果，其中1（逻辑）表示
  477.  *  对应象素为生长区域，0表示为非生长区域
  478.  *   区域生长一般包含三个比较重要的问题：
  479.  * 1. 种子点的选取
  480.  * 2. 生长准则
  481.  * 3. 终止条件
  482.  *  可以认为，这三个问题需要具体分析，而且每个问题解决的好坏直接关系到
  483.  *  区域生长的结果。
  484.  *  本函数的种子点选取为图像的中心，生长准则是相邻象素的象素值小于
  485.  *  nThreshold, 终止条件是一直进行到再没有满足生长准则需要的象素时为止
  486.  *
  487.  *************************************************************************
  488.  */
  489. void RegionGrow(CDib * pDib, unsigned char * pUnRegion, int nThreshold)
  490. {
  491. static int nDx[]={-1,0,1,0};
  492. static int nDy[]={0,1,0,-1};
  494. // 遍历图象的纵坐标
  495. // int y;
  497. // 遍历图象的横坐标
  498. // int x;
  500. // 图象的长宽大小
  501. CSize sizeImage = pDib->GetDimensions();
  502. int nWidth = sizeImage.cx ;
  503. int nHeight = sizeImage.cy ;
  505. // 图像在计算机在存储中的实际大小
  506. CSize sizeImageSave = pDib->GetDibSaveDim();
  508. // 图像在内存中每一行象素占用的实际空间
  509. int nSaveWidth = sizeImageSave.cx;
  511. // 初始化
  512. memset(pUnRegion,0,sizeof(unsigned char)*nWidth*nHeight);
  514. // 种子点
  515. int nSeedX, nSeedY;
  517. // 设置种子点为图像的中心
  518. nSeedX = nWidth /2 ;
  519. nSeedY = nHeight/2 ;
  521. // 定义堆栈，存储坐标
  522. int * pnGrowQueX ;
  523. int * pnGrowQueY ;
  525. // 分配空间
  526. pnGrowQueX = new int [nWidth*nHeight];
  527. pnGrowQueY = new int [nWidth*nHeight];
  529. // 图像数据的指针
  530. unsigned char *  pUnchInput =(unsigned char * )pDib->m_lpImage;
  532. // 定义堆栈的起点和终点
  533. // 当nStart=nEnd, 表示堆栈中只有一个点
  534. int nStart ;
  535. int nEnd   ;
  537. //初始化
  538. nStart = 0 ;
  539. nEnd   = 0 ;
  541. // 把种子点的坐标压入栈
  542. pnGrowQueX[nEnd] = nSeedX;
  543. pnGrowQueY[nEnd] = nSeedY;
  545. // 当前正在处理的象素
  546. int nCurrX ;
  547. int nCurrY ;
  549. // 循环控制变量
  550. int k ;
  552. // 图象的横纵坐标,用来对当前象素的4邻域进行遍历
  553. int xx;
  554. int yy;
  556. while (nStart<=nEnd)
  557. {
  558. // 当前种子点的坐标
  559. nCurrX = pnGrowQueX[nStart];
  560. nCurrY = pnGrowQueY[nStart];
  562. // 对当前点的4邻域进行遍历
  563. for (k=0; k<4; k++)
  564. {
  565. // 4邻域象素的坐标
  566. xx = nCurrX+nDx[k];
  567. yy = nCurrY+nDy[k];
  569. // 判断象素(xx，yy) 是否在图像内部
  570. // 判断象素(xx，yy) 是否已经处理过
  571. // pUnRegion[yy*nWidth+xx]==0 表示还没有处理
  573. // 生长条件：判断象素(xx，yy)和当前象素(nCurrX,nCurrY) 象素值差的绝对值
  574. if ( (xx < nWidth) && (xx>=0) && (yy<nHeight) && (yy>=0) 
  575.     && (pUnRegion[yy*nWidth+xx]==0) 
  576. && abs(pUnchInput[yy*nSaveWidth+xx] - pUnchInput[nCurrY*nSaveWidth+nCurrX])<nThreshold )
  577. {
  578. // 堆栈的尾部指针后移一位
  579. nEnd++;
  581. // 象素(xx，yy) 压入栈
  582. pnGrowQueX[nEnd] = xx;
  583. pnGrowQueY[nEnd] = yy;
  585. // 把象素(xx，yy)设置成逻辑1（255）
  586. // 同时也表明该象素处理过
  587. pUnRegion[yy*nWidth+xx] = 255 ;
  588. }
  589. }
  590. nStart++;
  591. }
  593. // 释放内存
  594. delete []pnGrowQueX;
  595. delete []pnGrowQueY;
  596.     pnGrowQueX = NULL ;
  597. pnGrowQueY = NULL ;
  598. }
  600. void DFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
  601. {
  602. double PI = 3.14159;
  603. //遍历图象的纵坐标
  604. int y;
  606. //遍历图象的横坐标
  607. int x;
  609. //频域的横坐标
  610. int m;
  612. //频域的纵坐标
  613. int n; 
  615. //图象的长宽大小
  616. CSize sizeImage = pDib->GetDimensions();
  617. int nWidth = sizeImage.cx ;
  618. int nHeight = sizeImage.cy ;
  620. //图像在计算机在存储中的实际大小
  621. CSize sizeImageSave = pDib->GetDibSaveDim();
  623. int nSaveWidth = sizeImageSave.cx;
  625. //图像数据的指针
  626. LPBYTE  pImageData = pDib->m_lpImage;
  628. //初始化结果数据
  629. for(n=0; n<nHeight ; n++ )
  630. for(m=0 ; m<nWidth ; m++ )
  631. {
  632. *( pTrRstRpart + n*nWidth + m ) =0;
  633. *( pTrRstIpart + n*nWidth + m ) =0;
  634. }
  635. double fCosTable;
  636. double fSinTable;
  637. int   nPxValue;
  639. fCosTable=0 ;
  640. nPxValue =0;
  642. double fTmpRstR;
  643. double fTmpRstI;
  644. for(n=0; n<nHeight ; n++ )
  645. for(m=0 ; m<nWidth ; m++ )
  646. {
  647. fTmpRstR=0;
  648. fTmpRstI=0;
  649. for(y=0; y<nHeight ; y++ )
  650. for(x=0 ; x<nWidth ; x++ )
  651. {
  652. fCosTable= cos( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  653. fSinTable= sin( -2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  654. nPxValue = *(pImageData+ y*nSaveWidth + x ) ;
  656. fTmpRstR+=fCosTable* nPxValue ;
  657. fTmpRstI+=fSinTable* nPxValue ;
  658. }
  659. *( pTrRstRpart + nWidth * n + m ) = fTmpRstR;
  660. *( pTrRstIpart + nWidth * n + m ) = fTmpRstI;
  661. }
  662. }
  665. void IDFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
  666. {
  667. double PI = 3.14159;
  668. //遍历图象的纵坐标
  669. int y;
  671. //遍历图象的横坐标
  672. int x;
  674. //频域的横坐标
  675. int m;
  677. //频域的纵坐标
  678. int n; 
  680. //图象的长宽大小
  681. CSize sizeImage = pDib->GetDimensions();
  682. int nWidth = sizeImage.cx ;
  683. int nHeight = sizeImage.cy ;
  685. //图像在计算机在存储中的实际大小
  686. CSize sizeImageSave = pDib->GetDibSaveDim();
  688. int nSaveWidth = sizeImageSave.cx;
  690. //图像数据的指针
  691. LPBYTE  pImageData = pDib->m_lpImage;
  693. double fCosTable;
  694. double fSinTable;
  695. fCosTable=0 ;
  696. fSinTable=0 ;
  698. double fTmpPxValue;
  699. double fRpartValue;
  700. double fIpartValue;
  701. fTmpPxValue=0;
  702. fRpartValue=0;
  703. fIpartValue=0;
  705. for(y=0; y<nHeight ; y++ )
  706. for(x=0 ; x<nWidth ; x++ )
  707. {
  708. fTmpPxValue=0;
  709. for(n=0; n<nHeight ; n++ )
  710. for(m=0 ; m<nWidth ; m++ )
  711. {
  712. fCosTable= cos( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  713. fSinTable= sin( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  714. fRpartValue=*(pTrRstRpart+ n*nHeight + m ) ;
  715. fIpartValue=*(pTrRstIpart+ n*nHeight + m ) ;
  717. fTmpPxValue+=fCosTable* fRpartValue-fSinTable*fIpartValue;
  718. }
  719. fTmpPxValue=fTmpPxValue/(nHeight*nWidth);
  720. *( pImageData + nSaveWidth * y + x) = (unsigned char) fTmpPxValue ;
  721. }
  722. }
  724. /*************************************************************************
  725.  *
  726.  * 函数名称：
  727.  *   SobelOperator()
  728.  *
  729.  * 输入参数:
  730.  *   CDib * pDib   - 指向CDib类的指针，含有原始图象信息
  731.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  732.  *
  733.  * 返回值:
  734.  *   无
  735.  *
  736.  * 说明:
  737.  *   Sobe算子
  738.  *
  739.  *   并行边界分割
  740.  *
  741.  *************************************************************************
  742.  */
  743. void SobelOperator(CDib * pDib, double * pdGrad)
  744. {
  745. // 遍历图象的纵坐标
  746. int y;
  748. // 遍历图象的横坐标
  749. int x;
  751. // 图象的长宽大小
  752. CSize sizeImage = pDib->GetDimensions();
  753. int nWidth = sizeImage.cx ;
  754. int nHeight = sizeImage.cy ;
  756. // 图像在计算机在存储中的实际大小
  757. CSize sizeImageSave = pDib->GetDibSaveDim();
  759. // 图像在内存中每一行象素占用的实际空间
  760. int nSaveWidth = sizeImageSave.cx;
  762. // 图像数据的指针
  763. LPBYTE  lpImage = pDib->m_lpImage;
  765. // 初始化
  766. for(y=0; y<nHeight ; y++ )
  767. for(x=0 ; x<nWidth ; x++ )
  768. {
  769. *(pdGrad+y*nWidth+x)=0;
  770. }
  772. // 设置模板系数
  773. static int nWeight[2][3][3] ;
  774. nWeight[0][0][0] = -1 ;   
  775. nWeight[0][0][1] =  0 ;   
  776. nWeight[0][0][2] =  1 ;   
  777. nWeight[0][1][0] = -2 ;   
  778. nWeight[0][1][1] =  0 ;   
  779. nWeight[0][1][2] =  2 ;   
  780. nWeight[0][2][0] = -1 ;   
  781. nWeight[0][2][1] =  0 ;   
  782. nWeight[0][2][2] =  1 ;   
  784. nWeight[1][0][0] =  1 ;   
  785. nWeight[1][0][1] =  2 ;   
  786. nWeight[1][0][2] =  1 ;   
  787. nWeight[1][1][0] =  0 ;   
  788. nWeight[1][1][1] =  0 ;   
  789. nWeight[1][1][2] =  0 ;   
  790. nWeight[1][2][0] = -1 ;   
  791. nWeight[1][2][1] = -2 ;   
  792. nWeight[1][2][2] = -1 ;   
  796. //这个变量用来表示Laplacian算子象素值
  797. int nTmp[3][3];
  799. // 临时变量
  800. double dGrad   ;
  801. double dGradOne;
  802. double dGradTwo;
  804. // 模板循环控制变量
  805. int yy ;
  806. int xx ;
  809. // 下面开始利用Prewitt算子进行计算，为了保证计算所需要的
  810. // 的数据位于图像数据的内部，下面的两重循环的条件是
  811. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  812. // 而不是x<nWidth
  813. for(y=1; y<nHeight-1 ; y++ )
  814. for(x=1 ; x<nWidth-1 ; x++ )
  815. {
  816. dGrad    = 0 ; 
  817. dGradOne = 0 ;
  818. dGradTwo = 0 ;
  819. // Laplacian算子需要的各点象素值
  821. // 模板第一行
  822. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  823. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  824. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  826. // 模板第二行
  827. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  828. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  829. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  831. // 模板第三行
  832. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  833. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  834. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  836. // 计算梯度
  837. for(yy=0; yy<3; yy++)
  838. for(xx=0; xx<3; xx++)
  839. {
  840. dGradOne += nTmp[yy][xx] * nWeight[0][yy][xx] ;
  841. dGradTwo += nTmp[yy][xx] * nWeight[1][yy][xx] ;
  842. }
  843. dGrad = dGradOne*dGradOne + dGradTwo*dGradTwo  ;
  844. dGrad = sqrt(dGrad) ;
  845. // 梯度值写入内存
  846. *(pdGrad+y*nWidth+x)=dGrad;
  847. }
  848. }
  850. /*************************************************************************
  851.  *
  852.  * 函数名称：
  853.  *   PrewittOperator()
  854.  *
  855.  * 输入参数:
  856.  *   CDib * pDib   - 指向CDib类的指针，含有原始图象信息
  857.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  858.  *
  859.  * 返回值:
  860.  *   无
  861.  *
  862.  * 说明:
  863.  *   Prewitt算子
  864.  *
  865.  *   并行边界分割
  866.  *
  867.  *************************************************************************
  868.  */
  869. void PrewittOperator(CDib * pDib, double * pdGrad)
  870. {
  871. // 遍历图象的纵坐标
  872. int y;
  874. // 遍历图象的横坐标
  875. int x;
  877. // 图象的长宽大小
  878. CSize sizeImage = pDib->GetDimensions();
  879. int nWidth = sizeImage.cx ;
  880. int nHeight = sizeImage.cy ;
  882. // 图像在计算机在存储中的实际大小
  883. CSize sizeImageSave = pDib->GetDibSaveDim();
  885. // 图像在内存中每一行象素占用的实际空间
  886. int nSaveWidth = sizeImageSave.cx;
  888. // 图像数据的指针
  889. LPBYTE  lpImage = pDib->m_lpImage;
  891. // 初始化
  892. for(y=0; y<nHeight ; y++ )
  893. for(x=0 ; x<nWidth ; x++ )
  894. {
  895. *(pdGrad+y*nWidth+x)=0;
  896. }
  898. // 设置模板系数
  899. static int nWeight[2][3][3] ;
  900. nWeight[0][0][0] = -1 ;   
  901. nWeight[0][0][1] =  0 ;   
  902. nWeight[0][0][2] =  1 ;   
  903. nWeight[0][1][0] = -1 ;   
  904. nWeight[0][1][1] =  0 ;   
  905. nWeight[0][1][2] =  1 ;   
  906. nWeight[0][2][0] = -1 ;   
  907. nWeight[0][2][1] =  0 ;   
  908. nWeight[0][2][2] =  1 ;   
  910. nWeight[1][0][0] =  1 ;   
  911. nWeight[1][0][1] =  1 ;   
  912. nWeight[1][0][2] =  1 ;   
  913. nWeight[1][1][0] =  0 ;   
  914. nWeight[1][1][1] =  0 ;   
  915. nWeight[1][1][2] =  0 ;   
  916. nWeight[1][2][0] = -1 ;   
  917. nWeight[1][2][1] = -1 ;   
  918. nWeight[1][2][2] = -1 ;   
  922. //这个变量用来表示Laplacian算子象素值
  923. int nTmp[3][3];
  925. // 临时变量
  926. double dGrad   ;
  927. double dGradOne;
  928. double dGradTwo;
  930. // 模板循环控制变量
  931. int yy ;
  932. int xx ;
  935. // 下面开始利用Prewitt算子进行计算，为了保证计算所需要的
  936. // 的数据位于图像数据的内部，下面的两重循环的条件是
  937. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  938. // 而不是x<nWidth
  939. for(y=1; y<nHeight-1 ; y++ )
  940. for(x=1 ; x<nWidth-1 ; x++ )
  941. {
  942. dGrad    = 0 ; 
  943. dGradOne = 0 ;
  944. dGradTwo = 0 ;
  945. // Laplacian算子需要的各点象素值
  947. // 模板第一行
  948. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  949. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  950. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  952. // 模板第二行
  953. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  954. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  955. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  957. // 模板第三行
  958. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  959. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  960. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  962. // 计算梯度
  963. for(yy=0; yy<3; yy++)
  964. for(xx=0; xx<3; xx++)
  965. {
  966. dGradOne += nTmp[yy][xx] * nWeight[0][yy][xx] ;
  967. dGradTwo += nTmp[yy][xx] * nWeight[1][yy][xx] ;
  968. }
  969. dGrad = dGradOne*dGradOne + dGradTwo*dGradTwo  ;
  970. dGrad = sqrt(dGrad) ;
  971. // 梯度值写入内存
  972. *(pdGrad+y*nWidth+x)=dGrad;
  973. }
  974. }
  976. /*************************************************************************
  977.  *
  978.  * 函数名称：
  979.  *   EdgeTrack()
  980.  *
  981.  * 输入参数:
  982.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  983.  *   unsigned char * pUnEdgeTrack - 指向边界跟踪结果的指针
  984.  *
  985.  * 返回值:
  986.  *   无
  987.  *
  988.  * 说明:
  989.  *   pUnEdgeTrack指针指向的数据区存储了边界跟踪的结果，其中1（逻辑）表示
  990.  *  对应象素为边界点，0表示为非边界点
  991.  *
  992.  *   串行边界分割
  993.  *
  994.  *************************************************************************
  995.  */
  996. void EdgeTrack(CDib * pDib, unsigned char * pUnEdgeTrack)
  997. {
  998. static int nDx[8]={-1,-1,-1, 0, 0, 1, 1, 1};
  999. static int nDy[8]={-1, 0, 1,-1, 1,-1, 0, 1};
  1001. // 遍历图象的纵坐标
  1002. int y;
  1004. // 遍历图象的横坐标
  1005. int x;
  1007. // 图象的长宽大小
  1008. CSize sizeImage = pDib->GetDimensions();
  1009. int nWidth = sizeImage.cx ;
  1010. int nHeight = sizeImage.cy ;
  1012. // 指向梯度数据的指针
  1013. double * pdGrad;
  1015. // 按照图像的大小开辟内存空间，存储梯度计算的结果
  1016. pdGrad=new double[nHeight*nWidth];
  1018.     // 调用Roberts算子求梯度
  1019. RobertsOperator(pDib, pdGrad);
  1021. // 定义当前象素梯度值
  1022. double dCurrGrad = 0;
  1024. // 定义最大梯度值
  1025. double dMaxGrad;
  1027. // 设置初值
  1028. dMaxGrad = 0;
  1030. // 最大梯度值对应的象素点坐标
  1031. int nPx;
  1032. int nPy;
  1034. nPx = 0;
  1035. nPy = 0;
  1037. // 求梯度最大值所在的象素点坐标
  1038. for(y=0; y<nHeight; y++)
  1039. {
  1040. for(x=0; x<nWidth; x++)
  1041. {
  1042. dCurrGrad = pdGrad[y*nWidth + x] ;
  1043. if( dMaxGrad< dCurrGrad )
  1044. {
  1045. dMaxGrad = dCurrGrad;
  1046. nPx = x ;
  1047. nPy = y ; 
  1048. }
  1049. }
  1050. }
  1052. // 初始化
  1053. memset(pUnEdgeTrack,0,sizeof(unsigned char)*nWidth*nHeight);
  1054. dCurrGrad = pdGrad[nPy*nWidth + nPx] ;
  1056. // 从（nPx，nPy）点开始进行边界跟踪
  1057. pUnEdgeTrack[nPy*nWidth + nPx] = 255 ;
  1059. // 循环变量，遍历当前象素的8邻域
  1060. int i ;
  1062. int yy;
  1063. int xx;
  1065. int nDetX;
  1066. int nDetY;
  1068. while(dCurrGrad>10)
  1069. {
  1070. // 设置当前点为边界点
  1071. pUnEdgeTrack[nPy*nWidth + nPx] = 255 ;
  1073. dMaxGrad = 0 ;
  1074. for(i=0; i<8; i++)
  1075. {
  1076. nDetX=nDx[i];
  1077. nDetY=nDy[i];
  1078. y = nPy + nDetY;
  1079. x = nPx + nDetX;
  1081. // 判断是否在图像内部
  1082. if(x>=0 && x<nWidth && y>=0 && y<nHeight)
  1083. {
  1084. if( ( pdGrad[y*nWidth + x] > dMaxGrad)  && ( pUnEdgeTrack[y*nWidth + x] == 0) )
  1085. {
  1086. dMaxGrad = pdGrad[y*nWidth + x] ;
  1087. yy = y;
  1088. xx = x;
  1089. }
  1090. }
  1091. }
  1092. // 下一个边界点的梯度，横纵坐标
  1093. dCurrGrad = dMaxGrad ;
  1094. nPy = yy;
  1095. nPx = xx;
  1096. }
  1098. //释放内存
  1099. delete pdGrad;
  1100. pdGrad = NULL;
  1102. }
  1104. /*************************************************************************
  1105.  *
  1106.  * 函数名称：
  1107.  *   MakeGauss()
  1108.  *
  1109.  * 输入参数:
  1110.  *   double sigma         - 高斯函数的标准差
  1111.  *   double **pdKernel - 指向高斯数据数组的指针
  1112.  *   int *pnWindowSize - 数据的长度
  1113.  *
  1114.  * 返回值:
  1115.  *   无
  1116.  *
  1117.  * 说明:
  1118.  *   这个函数可以生成一个一维的高斯函数的数字数据，理论上高斯数据的长度应
  1119.  *   该是无限长的，但是为了计算的简单和速度，实际的高斯数据只能是有限长的
  1120.  *   pnWindowSize就是数据长度
  1121.  *   
  1122.  *************************************************************************
  1123.  */
  1124. void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)
  1125. {
  1126. // 循环控制变量
  1127. int i   ;
  1129. // 数组的中心点
  1130. int nCenter;
  1132. // 数组的某一点到中心点的距离
  1133. double  dDis  ; 
  1135. double PI = 3.14159;
  1136. // 中间变量
  1137. double  dValue; 
  1138. double  dSum  ;
  1139. dSum = 0 ; 
  1141. // 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。
  1142. // 这些数据会覆盖绝大部分的滤波系数
  1143. *pnWindowSize = 1 + 2 * ceil(3 * sigma);
  1145. // 中心
  1146. nCenter = (*pnWindowSize) / 2;
  1148. // 分配内存
  1149. *pdKernel = new double[*pnWindowSize] ;
  1151. for(i=0; i< (*pnWindowSize); i++)
  1152. {
  1153. dDis = (double)(i - nCenter);
  1154. dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma)) / (sqrt(2 * PI) * sigma );
  1155. (*pdKernel)[i] = dValue ;
  1156. dSum += dValue;
  1157. }
  1159. // 归一化
  1160. for(i=0; i<(*pnWindowSize) ; i++)
  1161. {
  1162. (*pdKernel)[i] /= dSum;
  1163. }
  1164. }
  1166. /*************************************************************************
  1167.  *
  1168.  * 函数名称：
  1169.  *   GaussianSmooth()
  1170.  *
  1171.  * 输入参数:
  1172.  *   unsigned char * pUnchImg - 指向图象数据的指针
  1173.  *   int nWidth - 图象数据宽度
  1174.  *   int nHeight - 图象数据高度
  1175.  *   double dSigma - 高斯函数的标准差
  1176.  *   unsigned char * pUnchSmthdImg - 指向经过平滑之后的图象数据
  1177.  *
  1178.  * 返回值:
  1179.  *   无
  1180.  *
  1181.  * 说明:
  1182.  *   为了抑止噪声，采用高斯滤波对图象进行滤波，滤波先对x方向进行，然后对
  1183.  *   y方向进行。
  1184.  *   
  1185.  *************************************************************************
  1186.  */
  1187. void GaussianSmooth(unsigned char *pUnchImg, int nWidth, int nHeight, 
  1188. double sigma, unsigned char * pUnchSmthdImg)
  1189. {
  1190. // 循环控制变量
  1191.   int y;
  1192. int x;
  1194. int i;
  1196. // 高斯滤波器的数组长度
  1198. int nWindowSize;
  1200. //  窗口长度的1/2
  1201. int nHalfLen;   
  1203. // 一维高斯数据滤波器
  1204. double *pdKernel ;
  1206. // 高斯系数与图象数据的点乘
  1207. double  dDotMul     ;
  1209. // 高斯滤波系数的总和
  1210. double  dWeightSum     ;          
  1211.   
  1212. // 中间变量
  1213. double * pdTmp ;
  1215. // 分配内存
  1216. pdTmp = new double[nWidth*nHeight];
  1218. // 产生一维高斯数据滤波器
  1219. // MakeGauss(sigma, &dKernel, &nWindowSize);
  1220. MakeGauss(sigma, &pdKernel, &nWindowSize) ;
  1222. // MakeGauss返回窗口的长度，利用此变量计算窗口的半长
  1223. nHalfLen = nWindowSize / 2;
  1225.   // x方向进行滤波
  1226. for(y=0; y<nHeight; y++)
  1227. {
  1228. for(x=0; x<nWidth; x++)
  1229. {
  1230. dDotMul = 0;
  1231. dWeightSum = 0;
  1232. for(i=(-nHalfLen); i<=nHalfLen; i++)
  1233. {
  1234. // 判断是否在图象内部
  1235. if( (i+x) >= 0  && (i+x) < nWidth )
  1236. {
  1237. dDotMul += (double)pUnchImg[y*nWidth + (i+x)] * pdKernel[nHalfLen+i];
  1238. dWeightSum += pdKernel[nHalfLen+i];
  1239. }
  1240. }
  1241. pdTmp[y*nWidth + x] = dDotMul/dWeightSum ;
  1242. }
  1243. }
  1245. // y方向进行滤波
  1246. for(x=0; x<nWidth; x++)
  1247. {
  1248. for(y=0; y<nHeight; y++)
  1249. {
  1250. dDotMul = 0;
  1251. dWeightSum = 0;
  1252. for(i=(-nHalfLen); i<=nHalfLen; i++)
  1253. {
  1254. // 判断是否在图象内部
  1255. if( (i+y) >= 0  && (i+y) < nHeight )
  1256. {
  1257. dDotMul += (double)pdTmp[(y+i)*nWidth + x] * pdKernel[nHalfLen+i];
  1258. dWeightSum += pdKernel[nHalfLen+i];
  1259. }
  1260. }
  1261. pUnchSmthdImg[y*nWidth + x] = (unsigned char)(int)dDotMul/dWeightSum ;
  1262. }
  1263. }
  1265. // 释放内存
  1266. delete []pdKernel;
  1267. pdKernel = NULL ;
  1269. delete []pdTmp;
  1270. pdTmp = NULL;
  1271. }
  1274. /*************************************************************************
  1275.  *
  1276.  * 函数名称：
  1277.  *   DirGrad()
  1278.  *
  1279.  * 输入参数:
  1280.  *   unsigned char *pUnchSmthdImg - 经过高斯滤波后的图象
  1281.  *   int nWidht - 图象宽度
  1282.  *   int nHeight       - 图象高度
  1283.  *   int *pnGradX                         - x方向的方向导数
  1284.  *   int *pnGradY                         - y方向的方向导数
  1285.  * 返回值:
  1286.  *   无
  1287.  *
  1288.  * 说明:
  1289.  *   这个函数计算方向倒数，采用的微分算子是(-1 0 1) 和 (-1 0 1)'(转置)
  1290.  *   计算的时候对边界象素采用了特殊处理
  1291.  *   
  1292.  *   
  1293.  *************************************************************************
  1294.  */
  1295. void DirGrad(unsigned char *pUnchSmthdImg, int nWidth, int nHeight,
  1296.  int *pnGradX , int *pnGradY)
  1297. {
  1298. // 循环控制变量
  1299. int y ;
  1300. int x ;
  1302. // 计算x方向的方向导数，在边界出进行了处理，防止要访问的象素出界
  1303. for(y=0; y<nHeight; y++)
  1304. {
  1305. for(x=0; x<nWidth; x++)
  1306. {
  1307. pnGradX[y*nWidth+x] = (int) ( pUnchSmthdImg[y*nWidth+min(nWidth-1,x+1)]  
  1308.           - pUnchSmthdImg[y*nWidth+max(0,x-1)]     );
  1309. }
  1310. }
  1312. // 计算y方向的方向导数，在边界出进行了处理，防止要访问的象素出界
  1313. for(x=0; x<nWidth; x++)
  1314. {
  1315. for(y=0; y<nHeight; y++)
  1316. {
  1317. pnGradY[y*nWidth+x] = (int) ( pUnchSmthdImg[min(nHeight-1,y+1)*nWidth + x]  
  1318. - pUnchSmthdImg[max(0,y-1)*nWidth+ x ]     );
  1319. }
  1320. }
  1321. }
  1322. /*************************************************************************
  1323.  *
  1324.  * 函数名称：
  1325.  *   GradMagnitude()
  1326.  *
  1327.  * 输入参数:
  1328.  *   int *pnGradX                         - x方向的方向导数
  1329.  *   int *pnGradY                         - y方向的方向导数
  1330.  *   int nWidht - 图象宽度
  1331.  *   int nHeight       - 图象高度
  1332.  *   int *pnMag                           - 梯度幅度   
  1333.  *
  1334.  * 返回值:
  1335.  *   无
  1336.  *
  1337.  * 说明:
  1338.  *   这个函数利用方向倒数计算梯度幅度，方向倒数是DirGrad函数计算的结果
  1339.  *   
  1340.  *************************************************************************
  1341.  */
  1342. void GradMagnitude(int *pnGradX, int *pnGradY, int nWidth, int nHeight, int *pnMag)
  1343. {
  1345. // 循环控制变量
  1346. int y ;
  1347. int x ;
  1349. // 中间变量
  1350. double dSqtOne;
  1351. double dSqtTwo;
  1353. for(y=0; y<nHeight; y++)
  1354. {
  1355. for(x=0; x<nWidth; x++)
  1356. {
  1357. dSqtOne = pnGradX[y*nWidth + x] * pnGradX[y*nWidth + x];
  1358. dSqtTwo = pnGradY[y*nWidth + x] * pnGradY[y*nWidth + x];
  1359. pnMag[y*nWidth + x] = (int)(sqrt(dSqtOne + dSqtTwo) + 0.5);
  1360. }
  1361. }
  1362. }
  1363. /*************************************************************************
  1364.  *
  1365.  * 函数名称：
  1366.  *   NonmaxSuppress()
  1367.  *
  1368.  * 输入参数:
  1369.  *   int *pnMag                - 梯度图
  1370.  *   int *pnGradX  - x方向的方向导数
  1371.  *   int *pnGradY              - y方向的方向导数
  1372.  *   int nWidth                - 图象数据宽度
  1373.  *   int nHeight               - 图象数据高度
  1374.  *   unsigned char *pUnchRst   - 经过NonmaxSuppress处理后的结果
  1375.  *
  1376.  * 返回值:
  1377.  *   无
  1378.  *
  1379.  * 说明:
  1380.  *   抑止梯度图中非局部极值点的象素。
  1381.  *   
  1382.  *************************************************************************
  1383.  */
  1384. void NonmaxSuppress(int *pnMag, int *pnGradX, int *pnGradY, int nWidth, 
  1385. int nHeight, unsigned char *pUnchRst)
  1386. {
  1387. // 循环控制变量
  1388. int y ;
  1389. int x ;
  1390. int nPos;
  1392. // x方向梯度分量
  1393. int gx  ;
  1394. int gy  ;
  1396. // 临时变量
  1397. int g1, g2, g3, g4 ;
  1398. double weight  ;
  1399. double dTmp1   ;
  1400. double dTmp2   ;
  1401. double dTmp    ;
  1403. // 设置图象边缘部分为不可能的边界点
  1404. for(x=0; x<nWidth; x++)
  1405. {
  1406. pUnchRst[x] = 0 ;
  1407. pUnchRst[nHeight-1+x] = 0;
  1408. }
  1409. for(y=0; y<nHeight; y++)
  1410. {
  1411. pUnchRst[y*nWidth] = 0 ;
  1412. pUnchRst[y*nWidth + nWidth-1] = 0;
  1413. }
  1415. for(y=1; y<nHeight-1; y++)
  1416. {
  1417. for(x=1; x<nWidth-1; x++)
  1418. {
  1419. nPos = y*nWidth + x ;
  1421. // 如果当前象素的梯度幅度为0，则不是边界点
  1422. if(pnMag[nPos] == 0 )
  1423. {
  1424. pUnchRst[nPos] = 0 ;
  1425. }
  1426. else
  1427. {
  1428. // 当前象素的梯度幅度
  1429. dTmp = pnMag[nPos] ;
  1431. // x，y方向导数
  1432. gx = pnGradX[nPos]  ;
  1433. gy = pnGradY[nPos]  ;
  1435. // 如果方向导数y分量比x分量大，说明导数的方向更加“趋向”于y分量。
  1436. if (abs(gy) > abs(gx)) 
  1437. {
  1438. // 计算插值的比例
  1439. weight = fabs(gx)/fabs(gy); 
  1441. g2 = pnMag[nPos-nWidth] ; 
  1442. g4 = pnMag[nPos+nWidth] ;
  1444. // 如果x，y两个方向的方向导数的符号相同
  1445. // C是当前象素，与g1-g4的位置关系为：
  1446. // g1 g2 
  1447. //  C         
  1448. //  g4 g3 
  1449. if (gx*gy > 0) 
  1450. { 
  1451. g1 = pnMag[nPos-nWidth-1] ;
  1452. g3 = pnMag[nPos+nWidth+1] ;
  1453. } 
  1455. // 如果x，y两个方向的方向导数的符号相反
  1456. // C是当前象素，与g1-g4的位置关系为：
  1457. //    g2 g1
  1458. //  C         
  1459. // g3 g4  
  1460. else 
  1461. { 
  1462. g1 = pnMag[nPos-nWidth+1] ;
  1463. g3 = pnMag[nPos+nWidth-1] ;
  1464. } 
  1465. }
  1467. // 如果方向导数x分量比y分量大，说明导数的方向更加“趋向”于x分量
  1468. // 这个判断语句包含了x分量和y分量相等的情况
  1469. else
  1470. {
  1471. // 计算插值的比例
  1472. weight = fabs(gy)/fabs(gx); 
  1474. g2 = pnMag[nPos+1] ; 
  1475. g4 = pnMag[nPos-1] ;
  1477. // 如果x，y两个方向的方向导数的符号相同
  1478. // C是当前象素，与g1-g4的位置关系为：
  1479. // g3   
  1480. // g4 C g2       
  1481. //       g1
  1482. if (gx*gy > 0) 
  1483. {
  1484. g1 = pnMag[nPos+nWidth+1] ;
  1485. g3 = pnMag[nPos-nWidth-1] ;
  1486. } 
  1487. // 如果x，y两个方向的方向导数的符号相反
  1488. // C是当前象素，与g1-g4的位置关系为：
  1489. //      g1
  1490. // g4 C g2       
  1491. //  g3     
  1492. else 
  1493. { 
  1494. g1 = pnMag[nPos-nWidth+1] ;
  1495. g3 = pnMag[nPos+nWidth-1] ;
  1496. }
  1497. }
  1499. // 下面利用g1-g4对梯度进行插值
  1500. {
  1501. dTmp1 = weight*g1 + (1-weight)*g2 ;
  1502. dTmp2 = weight*g3 + (1-weight)*g4 ;
  1504. // 当前象素的梯度是局部的最大值
  1505. // 该点可能是个边界点
  1506. if(dTmp>=dTmp1 && dTmp>=dTmp2)
  1507. {
  1508. pUnchRst[nPos] = 128 ;
  1509. }
  1510. else
  1511. {
  1512. // 不可能是边界点
  1513. pUnchRst[nPos] = 0 ;
  1514. }
  1515. }
  1516. } //else
  1517. } // for
  1519. }
  1520. } 
  1522. /*************************************************************************
  1523.  *
  1524.  * 函数名称：
  1525.  *   TraceEdge()
  1526.  *
  1527.  * 输入参数:
  1528.  *   int    x - 跟踪起点的x坐标 
  1529.  *   int    y - 跟踪起点的y坐标
  1530.  *   int nLowThd - 判断一个点是否为边界点的低阈值
  1531.  *   unsigned char *pUnchEdge - 记录边界点的缓冲区
  1532.  *   int *pnMag               - 梯度幅度图
  1533.  *   int nWidth               - 图象数据宽度
  1534.  *
  1535.  * 返回值:
  1536.  *   无
  1537.  *
  1538.  * 说明:
  1539.  *   递归调用  
  1540.  *   从(x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pUnchEdge中没有处理并且
  1541.  *   可能是边界点的那些象素(=128),象素值为0表明该点不可能是边界点，象素值
  1542.  *   为255表明该点已经被设置为边界点，不必再考虑
  1543.  *   
  1544.  *   
  1545.  *************************************************************************
  1546.  */
  1547. void TraceEdge (int y, int x, int nLowThd, unsigned char *pUnchEdge, int *pnMag, int nWidth) 
  1548. { 
  1549. // 对8邻域象素进行查询
  1550. int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;
  1551. int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;
  1553. int yy ;
  1554. int xx ;
  1556. int k  ;
  1558. for(k=0; k<8; k++)
  1559. {
  1560. yy = y + yNb[k] ;
  1561. xx = x + xNb[k] ;
  1562. // 如果该象素为可能的边界点，又没有处理过
  1563. // 并且梯度大于阈值
  1564. if(pUnchEdge[yy*nWidth+xx] == 128  && pnMag[yy*nWidth+xx]>=nLowThd)
  1565. {
  1566. // 把该点设置成为边界点
  1567. pUnchEdge[yy*nWidth+xx] = 255 ;
  1569. // 以该点为中心进行跟踪
  1570. TraceEdge(yy, xx, nLowThd, pUnchEdge, pnMag, nWidth);
  1571. }
  1572. }
  1573. } 
  1575. /*************************************************************************
  1576.  *
  1577.  * 函数名称：
  1578.  *   EstimateThreshold()
  1579.  *
  1580.  * 输入参数:
  1581.  *   int *pnMag               - 梯度幅度图
  1582.  *  int nWidth               - 图象数据宽度
  1583.  *  int nHeight              - 图象数据高度
  1584.  *   int *pnThdHigh           - 高阈值
  1585.  *   int *pnThdLow            - 低阈值
  1586.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1587.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1588.  *   unsigned char *pUnchEdge - 经过non-maximum处理后的数据
  1589.  *
  1590.  * 返回值:
  1591.  *   无
  1592.  *
  1593.  * 说明:
  1594.  *   经过non-maximum处理后的数据pUnchEdge，统计pnMag的直方图，确定阈值。
  1595.  *   本函数中只是统计pUnchEdge中可能为边界点的那些象素。然后利用直方图，
  1596.  *   根据dRatioHigh设置高阈值，存储到pnThdHigh。利用dRationLow和高阈值，
  1597.  *   设置低阈值，存储到*pnThdLow。dRatioHigh是一种比例：表明梯度小于
  1598.  *   *pnThdHigh的象素数目占象素总数目的比例。dRationLow表明*pnThdHigh
  1599.  *   和*pnThdLow的比例，这个比例在canny算法的原文里，作者给出了一个区间。
  1600.  *
  1601.  *************************************************************************
  1602.  */
  1603. void EstimateThreshold(int *pnMag, int nWidth, int nHeight, int *pnThdHigh,int *pnThdLow, 
  1604.  unsigned char * pUnchEdge, double dRatioHigh, double dRationLow) 
  1605. { 
  1606. // 循环控制变量
  1607. int y;
  1608. int x;
  1609. int k;
  1611. // 该数组的大小和梯度值的范围有关，如果采用本程序的算法，那么梯度的范围不会超过pow(2,10)
  1612. int nHist[1024] ;
  1614. // 可能的边界数目
  1615. int nEdgeNb     ;
  1617. // 最大梯度值
  1618. int nMaxMag     ;
  1620. int nHighCount  ;
  1622. nMaxMag = 0     ; 
  1624. // 初始化
  1625. for(k=0; k<1024; k++) 
  1626. {
  1627. nHist[k] = 0; 
  1628. }
  1630. // 统计直方图，然后利用直方图计算阈值
  1631. for(y=0; y<nHeight; y++)
  1632. {
  1633. for(x=0; x<nWidth; x++)
  1634. {
  1635. // 只是统计那些可能是边界点，并且还没有处理过的象素
  1636. if(pUnchEdge[y*nWidth+x]==128)
  1637. {
  1638. nHist[ pnMag[y*nWidth+x] ]++;
  1639. }
  1640. }
  1641. }
  1643. nEdgeNb = nHist[0]  ;
  1644. nMaxMag = 0         ;
  1645. // 统计经过“非最大值抑止(non-maximum suppression)”后有多少象素
  1646. for(k=1; k<1024; k++)
  1647. {
  1648. if(nHist[k] != 0)
  1649. {
  1650. // 最大梯度值
  1651. nMaxMag = k;
  1652. }
  1654. // 梯度为0的点是不可能为边界点的
  1655. // 经过non-maximum suppression后有多少象素
  1656. nEdgeNb += nHist[k];
  1657. }
  1659. // 梯度比高阈值*pnThdHigh小的象素点总数目
  1660. nHighCount = (int)(dRatioHigh * nEdgeNb +0.5);
  1662. k = 1;
  1663. nEdgeNb = nHist[1];
  1665. // 计算高阈值
  1666. while( (k<(nMaxMag-1)) && (nEdgeNb < nHighCount) )
  1667. {
  1668. k++;
  1669. nEdgeNb += nHist[k];
  1670. }
  1672. // 设置高阈值
  1673. *pnThdHigh = k ;
  1675. // 设置低阈值
  1676. *pnThdLow  = (int)((*pnThdHigh) * dRationLow+ 0.5);
  1677. }
  1679. /*************************************************************************
  1680.  *
  1681.  * 函数名称：
  1682.  *   Hysteresis()
  1683.  *
  1684.  * 输入参数:
  1685.  *   int *pnMag               - 梯度幅度图
  1686.  *  int nWidth               - 图象数据宽度
  1687.  *  int nHeight              - 图象数据高度
  1688.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1689.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1690.  *   unsigned char *pUnchEdge - 记录边界点的缓冲区
  1691.  *
  1692.  * 返回值:
  1693.  *   无
  1694.  *
  1695.  * 说明:
  1696.  *   本函数实现类似“磁滞现象”的一个功能，也就是，先调用EstimateThreshold
  1697.  *   函数对经过non-maximum处理后的数据pUnchSpr估计一个高阈值，然后判断
  1698.  *   pUnchSpr中可能的边界象素(=128)的梯度是不是大于高阈值nThdHigh，如果比
  1699.  *   该阈值大，该点将作为一个边界的起点，调用TraceEdge函数，把对应该边界
  1700.  *   的所有象素找出来。最后，当整个搜索完毕时，如果还有象素没有被标志成
  1701.  *   边界点，那么就一定不是边界点。
  1702.  *   
  1703.  *************************************************************************
  1704.  */
  1705. void Hysteresis(int *pnMag, int nWidth, int nHeight, double dRatioLow, 
  1706. double dRatioHigh, unsigned char *pUnchEdge)
  1707. {
  1708. // 循环控制变量
  1709. int y;
  1710. int x;
  1712. int nThdHigh ;
  1713. int nThdLow  ;
  1715. int nPos;
  1717. // 估计TraceEdge需要的低阈值，以及Hysteresis函数使用的高阈值
  1718. EstimateThreshold(pnMag, nWidth, nHeight, &nThdHigh, 
  1719.                &nThdLow, pUnchEdge,dRatioHigh, dRatioLow);
  1721.   // 这个循环用来寻找大于nThdHigh的点，这些点被用来当作边界点，然后用
  1722. // TraceEdge函数来跟踪该点对应的边界
  1723.    for(y=0; y<nHeight; y++)
  1724.  {
  1725.       for(x=0; x<nWidth; x++)
  1726. {
  1727. nPos = y*nWidth + x ; 
  1729. // 如果该象素是可能的边界点，并且梯度大于高阈值，该象素作为
  1730. // 一个边界的起点
  1731. if((pUnchEdge[nPos] == 128) && (pnMag[nPos] >= nThdHigh))
  1732. {
  1733. // 设置该点为边界点
  1734. pUnchEdge[nPos] = 255;
  1735. TraceEdge(y, x, nThdLow, pUnchEdge, pnMag, nWidth);
  1736. }
  1737.       }
  1738.    }
  1740.  // 那些还没有被设置为边界点的象素已经不可能成为边界点
  1741.    for(y=0; y<nHeight; y++)
  1742.  {
  1743.  for(x=0; x<nWidth; x++)
  1744.  {
  1745.  nPos = y*nWidth + x ; 
  1746.  if(pUnchEdge[nPos] != 255)
  1747.  {
  1748.  // 设置为非边界点
  1749.  pUnchEdge[nPos] = 0 ;
  1750.  }
  1751.  }
  1752.  }
  1753. }
  1756. /*************************************************************************
  1757.  *
  1758.  * 函数名称：
  1759.  *   Canny()
  1760.  *
  1761.  * 输入参数:
  1762.  *   unsigned char *pUnchImage- 图象数据
  1763.  *  int nWidth               - 图象数据宽度
  1764.  *  int nHeight              - 图象数据高度
  1765.  *   double sigma             - 高斯滤波的标准方差
  1766.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1767.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1768.  *   unsigned char *pUnchEdge - canny算子计算后的分割图
  1769.  *
  1770.  * 返回值:
  1771.  *   无
  1772.  *
  1773.  * 说明:
  1774.  *   canny分割算子，计算的结果保存在pUnchEdge中，逻辑1(255)表示该点为
  1775.  *   边界点，逻辑0(0)表示该点为非边界点。该函数的参数sigma，dRatioLow
  1776.  *   dRatioHigh，是需要指定的。这些参数会影响分割后边界点数目的多少
  1777.  *************************************************************************
  1778.  */
  1779. void Canny(unsigned char *pUnchImage, int nWidth, int nHeight, double sigma,
  1780.  double dRatioLow, double dRatioHigh, unsigned char *pUnchEdge)
  1781. {
  1782. // 经过高斯滤波后的图象数据
  1783. unsigned char * pUnchSmooth ;
  1784.   
  1785. // 指向x方向导数的指针
  1786. int * pnGradX ; 
  1788. // 指向y方向导数的指针
  1789. int * pnGradY ;
  1791. // 梯度的幅度
  1792. int * pnGradMag ;
  1794. pUnchSmooth  = new unsigned char[nWidth*nHeight] ;
  1795. pnGradX      = new int [nWidth*nHeight]          ;
  1796. pnGradY      = new int [nWidth*nHeight]          ;
  1797. pnGradMag    = new int [nWidth*nHeight]          ;
  1799. // 对原图象进行滤波
  1800. GaussianSmooth(pUnchImage, nWidth, nHeight, sigma, pUnchSmooth) ;
  1802. // 计算方向导数
  1803. DirGrad(pUnchSmooth, nWidth, nHeight, pnGradX, pnGradY) ;
  1805. // 计算梯度的幅度
  1806. GradMagnitude(pnGradX, pnGradY, nWidth, nHeight, pnGradMag) ;
  1808. // 应用non-maximum 抑制
  1809. NonmaxSuppress(pnGradMag, pnGradX, pnGradY, nWidth, nHeight, pUnchEdge) ;
  1811. // 应用Hysteresis，找到所有的边界
  1812. Hysteresis(pnGradMag, nWidth, nHeight, dRatioLow, dRatioHigh, pUnchEdge);
  1815. // 释放内存
  1816. delete pnGradX      ;
  1817. pnGradX      = NULL ;
  1818. delete pnGradY      ;
  1819. pnGradY      = NULL ;
  1820. delete pnGradMag    ;
  1821. pnGradMag    = NULL ;
  1822. delete pUnchSmooth ;
  1823. pUnchSmooth  = NULL ;
  1824. }