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

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