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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 波变换 > 查看源码
detect.cpp
资源名称：imagefusionwavelettransform.rar　[点击查看]
上传用户：aqingfeng
上传日期：2014-03-25
资源大小：1839k
文件大小：42k
源码类别：

波变换

开发平台：

Visual C++

detect.cpp：源码内容
  1. #include "GlobalApi.h"
  2. #include "stdafx.h"
  3. #include "math.h"
  4. #include "小波变换Doc.h"
  5. #include <direct.h>
  6. #include <complex>
  7. #include <conio.h>
  8. #include "dib.h"
  11. /*************************************************************************
  12.  *
  13.  * 函数名称：
  14.  *   RegionSegFixThreshold()
  15.  *
  16.  * 输入参数:
  17.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  18.  *   int nThreshold - 区域分割的阈值
  19.  *
  20.  * 返回值:
  21.  *   无
  22.  *
  23.  * 说明:
  24.  *   1（逻辑）表示对应象素为前景区域，0表示背景
  25.  *   阈值分割的关键问题在于阈值的选取。阈值的选取一般应该视实际的应用而
  26.  *   灵活设定。
  27.  *
  28.  *************************************************************************
  29.  */
  30. void RegionSegFixThreshold(CDib * pDib, int nThreshold)
  31. {
  32. //遍历图象的纵坐标
  33. int y;
  35. //遍历图象的横坐标
  36. int x;
  38. //图象的长宽大小
  39. CSize sizeImage = pDib->GetDimensions();
  40. int nWidth = sizeImage.cx ;
  41. int nHeight = sizeImage.cy ;
  43. //图像在计算机在存储中的实际大小
  44. CSize sizeImageSave = pDib->GetDibSaveDim();
  46. //图像在内存中每一行象素占用的实际空间
  47. int nSaveWidth = sizeImageSave.cx;
  50. //图像数据的指针
  51. LPBYTE  pImageData = pDib->m_lpImage;
  53. for(y=0; y<nHeight ; y++ )
  54. for(x=0; x<nWidth ; x++ )
  55. {
  56. if( *(pImageData+y*nSaveWidth+x) < nThreshold)
  57. *(pImageData+y*nSaveWidth+x) = 0;
  58. }
  59. }
  61. /*************************************************************************
  62.  *
  63.  * 函数名称：
  64.  *   RegionSegAdaptive()
  65.  *
  66.  * 输入参数:
  67.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  68.  *
  69.  * 返回值:
  70.  *   无
  71.  *
  72.  * 说明:
  73.  *   1（逻辑）表示对应象素为前景区域，0表示背景
  74.  *   阈值分割的关键问题在于阈值的选取。阈值的选取一般应该视实际的应用而
  75.  *   灵活设定。本函数中，阈值不是固定的，而是根据图象象素的实际性质而设定的。
  76.  *   这个函数把图像分成四个子图象，然后计算每个子图象的均值，根据均值设置阈值
  77.  *   阈值只是应用在对应的子图象
  78.  *
  79.  *************************************************************************
  80.  */
  81. void RegionSegAdaptive(CDib * pDib)
  82. {
  83. //遍历图象的纵坐标
  84. int y;
  86. //遍历图象的横坐标
  87. int x;
  89. //图象的长宽大小
  90. CSize sizeImage = pDib->GetDimensions();
  91. int nWidth = sizeImage.cx ;
  92. int nHeight = sizeImage.cy ;
  94. //图像在计算机在存储中的实际大小
  95. CSize sizeImageSave = pDib->GetDibSaveDim();
  97. //图像在内存中每一行象素占用的实际空间
  98. int nSaveWidth = sizeImageSave.cx;
  100. //图像数据的指针
  101. LPBYTE  lpImage = pDib->m_lpImage;
  103. // 局部阈值
  104. int nThd[2][2] ;
  106. // 子图象的平均值
  107. int nLocAvg ;
  109. // 对左上图像逐点扫描：
  110. nLocAvg = 0 ;
  111. // y方向
  112. for(y=0; y<nHeight/2 ; y++ )
  113. {
  114. // x方向
  115. for(x=0; x<nWidth/2 ; x++ )
  116. {
  117. nLocAvg += lpImage[y*nSaveWidth + x];
  118. }
  119. }
  120. // 计算均值
  121. nLocAvg /= ( (nHeight/2) * (nWidth/2) ) ;
  123. // 设置阈值为子图象的平均值
  124. nThd[0][0] = nLocAvg ;
  126. // 对左上图像逐点扫描进行分割：
  127. // y方向
  128. for(y=0; y<nHeight/2 ; y++ )
  129. {
  130. // x方向
  131. for(x=0; x<nWidth/2 ; x++ )
  132. {
  133. if(lpImage[y*nSaveWidth + x]<nThd[0][0])
  134. lpImage[y*nSaveWidth + x] = 255 ;
  135. else
  136. {
  137. lpImage[y*nSaveWidth + x] = 0 ;
  138. }
  140. }
  141. }
  143. // =============================================
  144. // 对左下图像逐点扫描：
  145. nLocAvg = 0 ;
  146. // y方向
  147. for(y=nHeight/2; y<nHeight ; y++ )
  148. {
  149. // x方向
  150. for(x=0; x<nWidth/2 ; x++ )
  151. {
  152. nLocAvg += lpImage[y*nSaveWidth + x];
  153. }
  154. }
  155. // 计算均值
  156. nLocAvg /= ( (nHeight - nHeight/2) * (nWidth/2) ) ;
  158. // 设置阈值为子图象的平均值
  159. nThd[1][0] = nLocAvg ;
  161. // 对左下图像逐点扫描进行分割：
  162. // y方向
  163. for(y=nHeight/2; y<nHeight ; y++ )
  164. {
  165. // x方向
  166. for(x=0; x<nWidth/2 ; x++ )
  167. {
  168. if(lpImage[y*nSaveWidth + x]<nThd[1][0])
  169. lpImage[y*nSaveWidth + x] = 255 ;
  170. else
  171. {
  172. lpImage[y*nSaveWidth + x] = 0 ;
  173. }
  175. }
  176. }
  178. // =============================================
  179. // 对右上图像逐点扫描：
  180. nLocAvg = 0 ;
  181. // y方向
  182. for(y=0; y<nHeight/2 ; y++ )
  183. {
  184. // x方向
  185. for(x=nWidth/2; x<nWidth ; x++ )
  186. {
  187. nLocAvg += lpImage[y*nSaveWidth + x];
  188. }
  189. }
  190. // 计算均值
  191. nLocAvg /= ( (nHeight/2) * (nWidth - nWidth/2) ) ;
  193. // 设置阈值为子图象的平均值
  194. nThd[0][1] = nLocAvg ;
  196. // 对右上图像逐点扫描进行分割：
  197. // y方向
  198. for(y=0; y<nHeight/2 ; y++ )
  199. {
  200. // x方向
  201. for(x=nWidth/2; x<nWidth ; x++ )
  202. {
  203. if(lpImage[y*nSaveWidth + x]<nThd[0][1])
  204. lpImage[y*nSaveWidth + x] = 255 ;
  205. else
  206. {
  207. lpImage[y*nSaveWidth + x] = 0 ;
  208. }
  210. }
  211. }
  213. // =============================================
  214. // 对右下图像逐点扫描：
  215. nLocAvg = 0 ;
  216. // y方向
  217. for(y=nHeight/2; y<nHeight ; y++ )
  218. {
  219. // x方向
  220. for(x=nWidth/2; x<nWidth ; x++ )
  221. {
  222. nLocAvg += lpImage[y*nSaveWidth + x];
  223. }
  224. }
  225. // 计算均值
  226. nLocAvg /= ( (nHeight - nHeight/2) * (nWidth - nWidth/2) ) ;
  228. // 设置阈值为子图象的平均值
  229. nThd[1][1] = nLocAvg ;
  231. // 对右下图像逐点扫描进行分割：
  232. // y方向
  233. for(y=nHeight/2; y<nHeight ; y++ )
  234. {
  235. // x方向
  236. for(x=nWidth/2; x<nWidth ; x++ )
  237. {
  238. if(lpImage[y*nSaveWidth + x]<nThd[1][1])
  239. lpImage[y*nSaveWidth + x] = 255 ;
  240. else
  241. {
  242. lpImage[y*nSaveWidth + x] = 0 ;
  243. }
  244. }
  245. }
  247. // 为了显示方便显示，逻辑1用黑色显示，逻辑0用白色显示
  248. for(y=0; y<nHeight ; y++ )
  249. {
  250. // x方向
  251. for(x=0; x<nWidth ; x++ )
  252. {
  253. lpImage[y*nSaveWidth + x] = 255 - lpImage[y*nSaveWidth + x] ;
  254. }
  255. }
  256. }
  258. /*************************************************************************
  259.  *
  260.  * 函数名称：
  261.  *   RobertsOperator()
  262.  *
  263.  * 输入参数:
  264.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  265.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  266.  *
  267.  * 返回值:
  268.  *   无
  269.  *
  270.  * 说明:
  271.  *   Roberts算子
  272.  *
  273.  *************************************************************************
  274.  */
  275. void RobertsOperator(CDib * pDib, double * pdGrad)
  276. {
  277. // 遍历图象的纵坐标
  278. int y;
  280. // 遍历图象的横坐标
  281. int x;
  283. // 图象的长宽大小
  284. CSize sizeImage = pDib->GetDimensions();
  285. int nWidth = sizeImage.cx ;
  286. int nHeight = sizeImage.cy ;
  288. // 图像在计算机在存储中的实际大小
  289. CSize sizeImageSave = pDib->GetDibSaveDim();
  291. // 图像在内存中每一行象素占用的实际空间
  292. int nSaveWidth = sizeImageSave.cx;
  294. // 图像数据的指针
  295. LPBYTE  pImageData = pDib->m_lpImage;
  297. // 初始化
  298. for(y=0; y<nHeight ; y++ )
  299. for(x=0 ; x<nWidth ; x++ )
  300. {
  301. *(pdGrad+y*nWidth+x)=0;
  302. }
  305. // 下面开始利用Roberts算子进行计算，为了保证计算所需要的
  306. // 的数据位于图像数据的内部，下面的两重循环的条件是
  307. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  308. // 而不是x<nWidth
  310. //这两个变量用来表示Roberts算子第一个模板的两个象素值
  311. int nUpLeft;
  312. int nDownRight;
  314. // 这两个变量用来表示Roberts算子第二个模板的两个象素值
  315. int nUpRight;
  316. int nDownLeft;
  318. // 这两个变量用来表示Roberts算子计算的结果
  319. int nValueOne;
  320. int nValueTwo;
  322. // 临时变量
  323. double dGrad;
  325. for(y=0; y<nHeight-1 ; y++ )
  326. for(x=0 ; x<nWidth-1 ; x++ )
  327. {
  328. // Roberts算子第一个模板需要的象素值
  329. nUpLeft =*(pImageData+y*nSaveWidth+x) ; 
  330. nDownRight =*( pImageData+(y+1)*nSaveWidth+x+1 );
  331. nDownRight *=-1;
  333. //Roberts算子的第一个模板计算结果
  334. nValueOne =nUpLeft+nDownRight ;
  336. // Roberts算子第二个模板需要的象素值
  337. nUpRight =*( pImageData+y*nSaveWidth+x+1 ) ;
  338. nDownLeft =*( pImageData+(y+1)*nSaveWidth+x );
  339. nDownLeft *=-1;
  341. // Roberts算子的第二个模板计算结果
  342. nValueTwo =nUpRight+nDownLeft;
  344. // 计算两个偏导数的平方和
  345. dGrad=nValueOne*nValueOne + nValueTwo*nValueTwo;
  347. // 开方
  348. dGrad=pow(dGrad,0.5);
  350. // 范数采用欧式距离
  351. *(pdGrad+y*nWidth+x)=dGrad;
  352. }
  353. }
  354. /*************************************************************************
  355.  *
  356.  * 函数名称：
  357.  *   LaplacianOperator()
  358.  *
  359.  * 输入参数:
  360.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  361.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  362.  *
  363.  * 返回值:
  364.  *   无
  365.  *
  366.  * 说明:
  367.  *   LaplacianOperator算子,是二阶算子，不想Roberts算子那样需要两个模板计算
  368.  *   梯度，LaplacianOperator算子只要一个算子就可以计算梯度。但是因为利用了
  369.  *   二阶信息，对噪声比较敏感
  370.  *
  371.  *************************************************************************
  372.  */
  373. void LaplacianOperator(CDib * pDib, double * pdGrad)
  374. {
  375. // 遍历图象的纵坐标
  376. int y;
  378. // 遍历图象的横坐标
  379. int x;
  381. // 图象的长宽大小
  382. CSize sizeImage = pDib->GetDimensions();
  383. int nWidth = sizeImage.cx ;
  384. int nHeight = sizeImage.cy ;
  386. // 图像在计算机在存储中的实际大小
  387. CSize sizeImageSave = pDib->GetDibSaveDim();
  389. // 图像在内存中每一行象素占用的实际空间
  390. int nSaveWidth = sizeImageSave.cx;
  392. // 图像数据的指针
  393. LPBYTE  lpImage = pDib->m_lpImage;
  395. // 初始化
  396. for(y=0; y<nHeight ; y++ )
  397. for(x=0 ; x<nWidth ; x++ )
  398. {
  399. *(pdGrad+y*nWidth+x)=0;
  400. }
  402. // 设置模板系数
  403. static int nWeight[3][3] ;
  404. nWeight[0][0] = -1 ;   
  405. nWeight[0][1] = -1 ;   
  406. nWeight[0][2] = -1 ;   
  407. nWeight[1][0] = -1 ;   
  408. nWeight[1][1] =  8 ;   
  409. nWeight[1][2] = -1 ;   
  410. nWeight[2][0] = -1 ;   
  411. nWeight[2][1] = -1 ;   
  412. nWeight[2][2] = -1 ;   
  416. //这个变量用来表示Laplacian算子象素值
  417. int nTmp[3][3];
  419. // 临时变量
  420. double dGrad;
  422. // 模板循环控制变量
  423. int yy ;
  424. int xx ;
  427. // 下面开始利用Laplacian算子进行计算，为了保证计算所需要的
  428. // 的数据位于图像数据的内部，下面的两重循环的条件是
  429. // y<nHeight-2 而不是y<nHeight，相应的x方向也是x<nWidth-2
  430. // 而不是x<nWidth
  431. for(y=1; y<nHeight-2 ; y++ )
  432. for(x=1 ; x<nWidth-2 ; x++ )
  433. {
  434. dGrad = 0 ; 
  435. // Laplacian算子需要的各点象素值
  437. // 模板第一行
  438. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  439. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  440. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  442. // 模板第二行
  443. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  444. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  445. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  447. // 模板第三行
  448. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  449. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  450. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  452. // 计算梯度
  453. for(yy=0; yy<3; yy++)
  454. for(xx=0; xx<3; xx++)
  455. {
  456. dGrad += nTmp[yy][xx] * nWeight[yy][xx] ;
  457. }
  459. // 梯度值写入内存
  460. *(pdGrad+y*nWidth+x)=dGrad;
  461. }
  463. }
  465. /*************************************************************************
  466.  *
  467.  * 函数名称：
  468.  *   RegionGrow()
  469.  *
  470.  * 输入参数:
  471.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  472.  *   unsigned char * pUnRegion - 指向区域生长结果的指针
  473.  *
  474.  * 返回值:
  475.  *   无
  476.  *
  477.  * 说明:
  478.  *   pUnRegion指针指向的数据区存储了区域生长的结果，其中1（逻辑）表示
  479.  *  对应象素为生长区域，0表示为非生长区域
  480.  *   区域生长一般包含三个比较重要的问题：
  481.  * 1. 种子点的选取
  482.  * 2. 生长准则
  483.  * 3. 终止条件
  484.  *  可以认为，这三个问题需要具体分析，而且每个问题解决的好坏直接关系到
  485.  *  区域生长的结果。
  486.  *  本函数的种子点选取为图像的中心，生长准则是相邻象素的象素值小于
  487.  *  nThreshold, 终止条件是一直进行到再没有满足生长准则需要的象素时为止
  488.  *
  489.  *************************************************************************
  490.  */
  491. void RegionGrow(CDib * pDib, unsigned char * pUnRegion, int nThreshold)
  492. {
  493. static int nDx[]={-1,0,1,0};
  494. static int nDy[]={0,1,0,-1};
  496. // 遍历图象的纵坐标
  497. // int y;
  499. // 遍历图象的横坐标
  500. // int x;
  502. // 图象的长宽大小
  503. CSize sizeImage = pDib->GetDimensions();
  504. int nWidth = sizeImage.cx ;
  505. int nHeight = sizeImage.cy ;
  507. // 图像在计算机在存储中的实际大小
  508. CSize sizeImageSave = pDib->GetDibSaveDim();
  510. // 图像在内存中每一行象素占用的实际空间
  511. int nSaveWidth = sizeImageSave.cx;
  513. // 初始化
  514. memset(pUnRegion,0,sizeof(unsigned char)*nWidth*nHeight);
  516. // 种子点
  517. int nSeedX, nSeedY;
  519. // 设置种子点为图像的中心
  520. nSeedX = nWidth /2 ;
  521. nSeedY = nHeight/2 ;
  523. // 定义堆栈，存储坐标
  524. int * pnGrowQueX ;
  525. int * pnGrowQueY ;
  527. // 分配空间
  528. pnGrowQueX = new int [nWidth*nHeight];
  529. pnGrowQueY = new int [nWidth*nHeight];
  531. // 图像数据的指针
  532. unsigned char *  pUnchInput =(unsigned char * )pDib->m_lpImage;
  534. // 定义堆栈的起点和终点
  535. // 当nStart=nEnd, 表示堆栈中只有一个点
  536. int nStart ;
  537. int nEnd   ;
  539. //初始化
  540. nStart = 0 ;
  541. nEnd   = 0 ;
  543. // 把种子点的坐标压入栈
  544. pnGrowQueX[nEnd] = nSeedX;
  545. pnGrowQueY[nEnd] = nSeedY;
  547. // 当前正在处理的象素
  548. int nCurrX ;
  549. int nCurrY ;
  551. // 循环控制变量
  552. int k ;
  554. // 图象的横纵坐标,用来对当前象素的4邻域进行遍历
  555. int xx;
  556. int yy;
  558. while (nStart<=nEnd)
  559. {
  560. // 当前种子点的坐标
  561. nCurrX = pnGrowQueX[nStart];
  562. nCurrY = pnGrowQueY[nStart];
  564. // 对当前点的4邻域进行遍历
  565. for (k=0; k<4; k++)
  566. {
  567. // 4邻域象素的坐标
  568. xx = nCurrX+nDx[k];
  569. yy = nCurrY+nDy[k];
  571. // 判断象素(xx，yy) 是否在图像内部
  572. // 判断象素(xx，yy) 是否已经处理过
  573. // pUnRegion[yy*nWidth+xx]==0 表示还没有处理
  575. // 生长条件：判断象素(xx，yy)和当前象素(nCurrX,nCurrY) 象素值差的绝对值
  576. if ( (xx < nWidth) && (xx>=0) && (yy<nHeight) && (yy>=0) 
  577.     && (pUnRegion[yy*nWidth+xx]==0) 
  578. && abs(pUnchInput[yy*nSaveWidth+xx] - pUnchInput[nCurrY*nSaveWidth+nCurrX])<nThreshold )
  579. {
  580. // 堆栈的尾部指针后移一位
  581. nEnd++;
  583. // 象素(xx，yy) 压入栈
  584. pnGrowQueX[nEnd] = xx;
  585. pnGrowQueY[nEnd] = yy;
  587. // 把象素(xx，yy)设置成逻辑1（255）
  588. // 同时也表明该象素处理过
  589. pUnRegion[yy*nWidth+xx] = 255 ;
  590. }
  591. }
  592. nStart++;
  593. }
  595. // 释放内存
  596. delete []pnGrowQueX;
  597. delete []pnGrowQueY;
  598.     pnGrowQueX = NULL ;
  599. pnGrowQueY = NULL ;
  600. }
  602. void DFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
  603. {
  604. double PI = 3.14159;
  605. //遍历图象的纵坐标
  606. int y;
  608. //遍历图象的横坐标
  609. int x;
  611. //频域的横坐标
  612. int m;
  614. //频域的纵坐标
  615. int n; 
  617. //图象的长宽大小
  618. CSize sizeImage = pDib->GetDimensions();
  619. int nWidth = sizeImage.cx ;
  620. int nHeight = sizeImage.cy ;
  622. //图像在计算机在存储中的实际大小
  623. CSize sizeImageSave = pDib->GetDibSaveDim();
  625. int nSaveWidth = sizeImageSave.cx;
  627. //图像数据的指针
  628. LPBYTE  pImageData = pDib->m_lpImage;
  630. //初始化结果数据
  631. for(n=0; n<nHeight ; n++ )
  632. for(m=0 ; m<nWidth ; m++ )
  633. {
  634. *( pTrRstRpart + n*nWidth + m ) =0;
  635. *( pTrRstIpart + n*nWidth + m ) =0;
  636. }
  637. double fCosTable;
  638. double fSinTable;
  639. int   nPxValue;
  641. fCosTable=0 ;
  642. nPxValue =0;
  644. double fTmpRstR;
  645. double fTmpRstI;
  646. for(n=0; n<nHeight ; n++ )
  647. for(m=0 ; m<nWidth ; m++ )
  648. {
  649. fTmpRstR=0;
  650. fTmpRstI=0;
  651. for(y=0; y<nHeight ; y++ )
  652. for(x=0 ; x<nWidth ; x++ )
  653. {
  654. fCosTable= cos( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  655. fSinTable= sin( -2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  656. nPxValue = *(pImageData+ y*nSaveWidth + x ) ;
  658. fTmpRstR+=fCosTable* nPxValue ;
  659. fTmpRstI+=fSinTable* nPxValue ;
  660. }
  661. *( pTrRstRpart + nWidth * n + m ) = fTmpRstR;
  662. *( pTrRstIpart + nWidth * n + m ) = fTmpRstI;
  663. }
  664. }
  667. void IDFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
  668. {
  669. double PI = 3.14159;
  670. //遍历图象的纵坐标
  671. int y;
  673. //遍历图象的横坐标
  674. int x;
  676. //频域的横坐标
  677. int m;
  679. //频域的纵坐标
  680. int n; 
  682. //图象的长宽大小
  683. CSize sizeImage = pDib->GetDimensions();
  684. int nWidth = sizeImage.cx ;
  685. int nHeight = sizeImage.cy ;
  687. //图像在计算机在存储中的实际大小
  688. CSize sizeImageSave = pDib->GetDibSaveDim();
  690. int nSaveWidth = sizeImageSave.cx;
  692. //图像数据的指针
  693. LPBYTE  pImageData = pDib->m_lpImage;
  695. double fCosTable;
  696. double fSinTable;
  697. fCosTable=0 ;
  698. fSinTable=0 ;
  700. double fTmpPxValue;
  701. double fRpartValue;
  702. double fIpartValue;
  703. fTmpPxValue=0;
  704. fRpartValue=0;
  705. fIpartValue=0;
  707. for(y=0; y<nHeight ; y++ )
  708. for(x=0 ; x<nWidth ; x++ )
  709. {
  710. fTmpPxValue=0;
  711. for(n=0; n<nHeight ; n++ )
  712. for(m=0 ; m<nWidth ; m++ )
  713. {
  714. fCosTable= cos( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  715. fSinTable= sin( 2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
  716. fRpartValue=*(pTrRstRpart+ n*nHeight + m ) ;
  717. fIpartValue=*(pTrRstIpart+ n*nHeight + m ) ;
  719. fTmpPxValue+=fCosTable* fRpartValue-fSinTable*fIpartValue;
  720. }
  721. fTmpPxValue=fTmpPxValue/(nHeight*nWidth);
  722. *( pImageData + nSaveWidth * y + x) = (unsigned char) fTmpPxValue ;
  723. }
  724. }
  726. /*************************************************************************
  727.  *
  728.  * 函数名称：
  729.  *   SobelOperator()
  730.  *
  731.  * 输入参数:
  732.  *   CDib * pDib   - 指向CDib类的指针，含有原始图象信息
  733.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  734.  *
  735.  * 返回值:
  736.  *   无
  737.  *
  738.  * 说明:
  739.  *   Sobe算子
  740.  *
  741.  *   并行边界分割
  742.  *
  743.  *************************************************************************
  744.  */
  745. void SobelOperator(CDib * pDib, double * pdGrad)
  746. {
  747. // 遍历图象的纵坐标
  748. int y;
  750. // 遍历图象的横坐标
  751. int x;
  753. // 图象的长宽大小
  754. CSize sizeImage = pDib->GetDimensions();
  755. int nWidth = sizeImage.cx ;
  756. int nHeight = sizeImage.cy ;
  758. // 图像在计算机在存储中的实际大小
  759. CSize sizeImageSave = pDib->GetDibSaveDim();
  761. // 图像在内存中每一行象素占用的实际空间
  762. int nSaveWidth = sizeImageSave.cx;
  764. // 图像数据的指针
  765. LPBYTE  lpImage = pDib->m_lpImage;
  767. // 初始化
  768. for(y=0; y<nHeight ; y++ )
  769. for(x=0 ; x<nWidth ; x++ )
  770. {
  771. *(pdGrad+y*nWidth+x)=0;
  772. }
  774. // 设置模板系数
  775. static int nWeight[2][3][3] ;
  776. nWeight[0][0][0] = -1 ;   
  777. nWeight[0][0][1] =  0 ;   
  778. nWeight[0][0][2] =  1 ;   
  779. nWeight[0][1][0] = -2 ;   
  780. nWeight[0][1][1] =  0 ;   
  781. nWeight[0][1][2] =  2 ;   
  782. nWeight[0][2][0] = -1 ;   
  783. nWeight[0][2][1] =  0 ;   
  784. nWeight[0][2][2] =  1 ;   
  786. nWeight[1][0][0] =  1 ;   
  787. nWeight[1][0][1] =  2 ;   
  788. nWeight[1][0][2] =  1 ;   
  789. nWeight[1][1][0] =  0 ;   
  790. nWeight[1][1][1] =  0 ;   
  791. nWeight[1][1][2] =  0 ;   
  792. nWeight[1][2][0] = -1 ;   
  793. nWeight[1][2][1] = -2 ;   
  794. nWeight[1][2][2] = -1 ;   
  798. //这个变量用来表示Laplacian算子象素值
  799. int nTmp[3][3];
  801. // 临时变量
  802. double dGrad   ;
  803. double dGradOne;
  804. double dGradTwo;
  806. // 模板循环控制变量
  807. int yy ;
  808. int xx ;
  811. // 下面开始利用Prewitt算子进行计算，为了保证计算所需要的
  812. // 的数据位于图像数据的内部，下面的两重循环的条件是
  813. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  814. // 而不是x<nWidth
  815. for(y=1; y<nHeight-1 ; y++ )
  816. for(x=1 ; x<nWidth-1 ; x++ )
  817. {
  818. dGrad    = 0 ; 
  819. dGradOne = 0 ;
  820. dGradTwo = 0 ;
  821. // Laplacian算子需要的各点象素值
  823. // 模板第一行
  824. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  825. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  826. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  828. // 模板第二行
  829. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  830. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  831. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  833. // 模板第三行
  834. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  835. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  836. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  838. // 计算梯度
  839. for(yy=0; yy<3; yy++)
  840. for(xx=0; xx<3; xx++)
  841. {
  842. dGradOne += nTmp[yy][xx] * nWeight[0][yy][xx] ;
  843. dGradTwo += nTmp[yy][xx] * nWeight[1][yy][xx] ;
  844. }
  845. dGrad = dGradOne*dGradOne + dGradTwo*dGradTwo  ;
  846. dGrad = sqrt(dGrad) ;
  847. // 梯度值写入内存
  848. *(pdGrad+y*nWidth+x)=dGrad;
  849. }
  850. }
  852. /*************************************************************************
  853.  *
  854.  * 函数名称：
  855.  *   PrewittOperator()
  856.  *
  857.  * 输入参数:
  858.  *   CDib * pDib   - 指向CDib类的指针，含有原始图象信息
  859.  *   double * pdGrad - 指向梯度数据的指针，含有图像的梯度信息
  860.  *
  861.  * 返回值:
  862.  *   无
  863.  *
  864.  * 说明:
  865.  *   Prewitt算子
  866.  *
  867.  *   并行边界分割
  868.  *
  869.  *************************************************************************
  870.  */
  871. void PrewittOperator(CDib * pDib, double * pdGrad)
  872. {
  873. // 遍历图象的纵坐标
  874. int y;
  876. // 遍历图象的横坐标
  877. int x;
  879. // 图象的长宽大小
  880. CSize sizeImage = pDib->GetDimensions();
  881. int nWidth = sizeImage.cx ;
  882. int nHeight = sizeImage.cy ;
  884. // 图像在计算机在存储中的实际大小
  885. CSize sizeImageSave = pDib->GetDibSaveDim();
  887. // 图像在内存中每一行象素占用的实际空间
  888. int nSaveWidth = sizeImageSave.cx;
  890. // 图像数据的指针
  891. LPBYTE  lpImage = pDib->m_lpImage;
  893. // 初始化
  894. for(y=0; y<nHeight ; y++ )
  895. for(x=0 ; x<nWidth ; x++ )
  896. {
  897. *(pdGrad+y*nWidth+x)=0;
  898. }
  900. // 设置模板系数
  901. static int nWeight[2][3][3] ;
  902. nWeight[0][0][0] = -1 ;   
  903. nWeight[0][0][1] =  0 ;   
  904. nWeight[0][0][2] =  1 ;   
  905. nWeight[0][1][0] = -1 ;   
  906. nWeight[0][1][1] =  0 ;   
  907. nWeight[0][1][2] =  1 ;   
  908. nWeight[0][2][0] = -1 ;   
  909. nWeight[0][2][1] =  0 ;   
  910. nWeight[0][2][2] =  1 ;   
  912. nWeight[1][0][0] =  1 ;   
  913. nWeight[1][0][1] =  1 ;   
  914. nWeight[1][0][2] =  1 ;   
  915. nWeight[1][1][0] =  0 ;   
  916. nWeight[1][1][1] =  0 ;   
  917. nWeight[1][1][2] =  0 ;   
  918. nWeight[1][2][0] = -1 ;   
  919. nWeight[1][2][1] = -1 ;   
  920. nWeight[1][2][2] = -1 ;   
  924. //这个变量用来表示Laplacian算子象素值
  925. int nTmp[3][3];
  927. // 临时变量
  928. double dGrad   ;
  929. double dGradOne;
  930. double dGradTwo;
  932. // 模板循环控制变量
  933. int yy ;
  934. int xx ;
  937. // 下面开始利用Prewitt算子进行计算，为了保证计算所需要的
  938. // 的数据位于图像数据的内部，下面的两重循环的条件是
  939. // y<nHeight-1 而不是y<nHeight，相应的x方向也是x<nWidth-1
  940. // 而不是x<nWidth
  941. for(y=1; y<nHeight-1 ; y++ )
  942. for(x=1 ; x<nWidth-1 ; x++ )
  943. {
  944. dGrad    = 0 ; 
  945. dGradOne = 0 ;
  946. dGradTwo = 0 ;
  947. // Laplacian算子需要的各点象素值
  949. // 模板第一行
  950. nTmp[0][0] = lpImage[(y-1)*nSaveWidth + x - 1 ] ; 
  951. nTmp[0][1] = lpImage[(y-1)*nSaveWidth + x     ] ; 
  952. nTmp[0][2] = lpImage[(y-1)*nSaveWidth + x + 1 ] ; 
  954. // 模板第二行
  955. nTmp[1][0] = lpImage[y*nSaveWidth + x - 1 ] ; 
  956. nTmp[1][1] = lpImage[y*nSaveWidth + x     ] ; 
  957. nTmp[1][2] = lpImage[y*nSaveWidth + x + 1 ] ; 
  959. // 模板第三行
  960. nTmp[2][0] = lpImage[(y+1)*nSaveWidth + x - 1 ] ; 
  961. nTmp[2][1] = lpImage[(y+1)*nSaveWidth + x     ] ; 
  962. nTmp[2][2] = lpImage[(y+1)*nSaveWidth + x + 1 ] ; 
  964. // 计算梯度
  965. for(yy=0; yy<3; yy++)
  966. for(xx=0; xx<3; xx++)
  967. {
  968. dGradOne += nTmp[yy][xx] * nWeight[0][yy][xx] ;
  969. dGradTwo += nTmp[yy][xx] * nWeight[1][yy][xx] ;
  970. }
  971. dGrad = dGradOne*dGradOne + dGradTwo*dGradTwo  ;
  972. dGrad = sqrt(dGrad) ;
  973. // 梯度值写入内存
  974. *(pdGrad+y*nWidth+x)=dGrad;
  975. }
  976. }
  978. /*************************************************************************
  979.  *
  980.  * 函数名称：
  981.  *   EdgeTrack()
  982.  *
  983.  * 输入参数:
  984.  *   CDib * pDib - 指向CDib类的指针，含有原始图象信息
  985.  *   unsigned char * pUnEdgeTrack - 指向边界跟踪结果的指针
  986.  *
  987.  * 返回值:
  988.  *   无
  989.  *
  990.  * 说明:
  991.  *   pUnEdgeTrack指针指向的数据区存储了边界跟踪的结果，其中1（逻辑）表示
  992.  *  对应象素为边界点，0表示为非边界点
  993.  *
  994.  *   串行边界分割
  995.  *
  996.  *************************************************************************
  997.  */
  998. void EdgeTrack(CDib * pDib, unsigned char * pUnEdgeTrack)
  999. {
  1000. static int nDx[8]={-1,-1,-1, 0, 0, 1, 1, 1};
  1001. static int nDy[8]={-1, 0, 1,-1, 1,-1, 0, 1};
  1003. // 遍历图象的纵坐标
  1004. int y;
  1006. // 遍历图象的横坐标
  1007. int x;
  1009. // 图象的长宽大小
  1010. CSize sizeImage = pDib->GetDimensions();
  1011. int nWidth = sizeImage.cx ;
  1012. int nHeight = sizeImage.cy ;
  1014. // 指向梯度数据的指针
  1015. double * pdGrad;
  1017. // 按照图像的大小开辟内存空间，存储梯度计算的结果
  1018. pdGrad=new double[nHeight*nWidth];
  1020.     // 调用Roberts算子求梯度
  1021. RobertsOperator(pDib, pdGrad);
  1023. // 定义当前象素梯度值
  1024. double dCurrGrad = 0;
  1026. // 定义最大梯度值
  1027. double dMaxGrad;
  1029. // 设置初值
  1030. dMaxGrad = 0;
  1032. // 最大梯度值对应的象素点坐标
  1033. int nPx;
  1034. int nPy;
  1036. nPx = 0;
  1037. nPy = 0;
  1039. // 求梯度最大值所在的象素点坐标
  1040. for(y=0; y<nHeight; y++)
  1041. {
  1042. for(x=0; x<nWidth; x++)
  1043. {
  1044. dCurrGrad = pdGrad[y*nWidth + x] ;
  1045. if( dMaxGrad< dCurrGrad )
  1046. {
  1047. dMaxGrad = dCurrGrad;
  1048. nPx = x ;
  1049. nPy = y ; 
  1050. }
  1051. }
  1052. }
  1054. // 初始化
  1055. memset(pUnEdgeTrack,0,sizeof(unsigned char)*nWidth*nHeight);
  1056. dCurrGrad = pdGrad[nPy*nWidth + nPx] ;
  1058. // 从（nPx，nPy）点开始进行边界跟踪
  1059. pUnEdgeTrack[nPy*nWidth + nPx] = 255 ;
  1061. // 循环变量，遍历当前象素的8邻域
  1062. int i ;
  1064. int yy;
  1065. int xx;
  1067. int nDetX;
  1068. int nDetY;
  1070. while(dCurrGrad>10)
  1071. {
  1072. // 设置当前点为边界点
  1073. pUnEdgeTrack[nPy*nWidth + nPx] = 255 ;
  1075. dMaxGrad = 0 ;
  1076. for(i=0; i<8; i++)
  1077. {
  1078. nDetX=nDx[i];
  1079. nDetY=nDy[i];
  1080. y = nPy + nDetY;
  1081. x = nPx + nDetX;
  1083. // 判断是否在图像内部
  1084. if(x>=0 && x<nWidth && y>=0 && y<nHeight)
  1085. {
  1086. if( ( pdGrad[y*nWidth + x] > dMaxGrad)  && ( pUnEdgeTrack[y*nWidth + x] == 0) )
  1087. {
  1088. dMaxGrad = pdGrad[y*nWidth + x] ;
  1089. yy = y;
  1090. xx = x;
  1091. }
  1092. }
  1093. }
  1094. // 下一个边界点的梯度，横纵坐标
  1095. dCurrGrad = dMaxGrad ;
  1096. nPy = yy;
  1097. nPx = xx;
  1098. }
  1100. //释放内存
  1101. delete pdGrad;
  1102. pdGrad = NULL;
  1104. }
  1106. /*************************************************************************
  1107.  *
  1108.  * 函数名称：
  1109.  *   MakeGauss()
  1110.  *
  1111.  * 输入参数:
  1112.  *   double sigma         - 高斯函数的标准差
  1113.  *   double **pdKernel - 指向高斯数据数组的指针
  1114.  *   int *pnWindowSize - 数据的长度
  1115.  *
  1116.  * 返回值:
  1117.  *   无
  1118.  *
  1119.  * 说明:
  1120.  *   这个函数可以生成一个一维的高斯函数的数字数据，理论上高斯数据的长度应
  1121.  *   该是无限长的，但是为了计算的简单和速度，实际的高斯数据只能是有限长的
  1122.  *   pnWindowSize就是数据长度
  1123.  *   
  1124.  *************************************************************************
  1125.  */
  1126. void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)
  1127. {
  1128. // 循环控制变量
  1129. int i   ;
  1131. // 数组的中心点
  1132. int nCenter;
  1134. // 数组的某一点到中心点的距离
  1135. double  dDis  ; 
  1137. double PI = 3.14159;
  1138. // 中间变量
  1139. double  dValue; 
  1140. double  dSum  ;
  1141. dSum = 0 ; 
  1143. // 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。
  1144. // 这些数据会覆盖绝大部分的滤波系数
  1145. *pnWindowSize = 1 + 2 * ceil(3 * sigma);
  1147. // 中心
  1148. nCenter = (*pnWindowSize) / 2;
  1150. // 分配内存
  1151. *pdKernel = new double[*pnWindowSize] ;
  1153. for(i=0; i< (*pnWindowSize); i++)
  1154. {
  1155. dDis = (double)(i - nCenter);
  1156. dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma)) / (sqrt(2 * PI) * sigma );
  1157. (*pdKernel)[i] = dValue ;
  1158. dSum += dValue;
  1159. }
  1161. // 归一化
  1162. for(i=0; i<(*pnWindowSize) ; i++)
  1163. {
  1164. (*pdKernel)[i] /= dSum;
  1165. }
  1166. }
  1168. /*************************************************************************
  1169.  *
  1170.  * 函数名称：
  1171.  *   GaussianSmooth()
  1172.  *
  1173.  * 输入参数:
  1174.  *   unsigned char * pUnchImg - 指向图象数据的指针
  1175.  *   int nWidth - 图象数据宽度
  1176.  *   int nHeight - 图象数据高度
  1177.  *   double dSigma - 高斯函数的标准差
  1178.  *   unsigned char * pUnchSmthdImg - 指向经过平滑之后的图象数据
  1179.  *
  1180.  * 返回值:
  1181.  *   无
  1182.  *
  1183.  * 说明:
  1184.  *   为了抑止噪声，采用高斯滤波对图象进行滤波，滤波先对x方向进行，然后对
  1185.  *   y方向进行。
  1186.  *   
  1187.  *************************************************************************
  1188.  */
  1189. void GaussianSmooth(unsigned char *pUnchImg, int nWidth, int nHeight, 
  1190. double sigma, unsigned char * pUnchSmthdImg)
  1191. {
  1192. // 循环控制变量
  1193.   int y;
  1194. int x;
  1196. int i;
  1198. // 高斯滤波器的数组长度
  1200. int nWindowSize;
  1202. //  窗口长度的1/2
  1203. int nHalfLen;   
  1205. // 一维高斯数据滤波器
  1206. double *pdKernel ;
  1208. // 高斯系数与图象数据的点乘
  1209. double  dDotMul     ;
  1211. // 高斯滤波系数的总和
  1212. double  dWeightSum     ;          
  1213.   
  1214. // 中间变量
  1215. double * pdTmp ;
  1217. // 分配内存
  1218. pdTmp = new double[nWidth*nHeight];
  1220. // 产生一维高斯数据滤波器
  1221. // MakeGauss(sigma, &dKernel, &nWindowSize);
  1222. MakeGauss(sigma, &pdKernel, &nWindowSize) ;
  1224. // MakeGauss返回窗口的长度，利用此变量计算窗口的半长
  1225. nHalfLen = nWindowSize / 2;
  1227.   // x方向进行滤波
  1228. for(y=0; y<nHeight; y++)
  1229. {
  1230. for(x=0; x<nWidth; x++)
  1231. {
  1232. dDotMul = 0;
  1233. dWeightSum = 0;
  1234. for(i=(-nHalfLen); i<=nHalfLen; i++)
  1235. {
  1236. // 判断是否在图象内部
  1237. if( (i+x) >= 0  && (i+x) < nWidth )
  1238. {
  1239. dDotMul += (double)pUnchImg[y*nWidth + (i+x)] * pdKernel[nHalfLen+i];
  1240. dWeightSum += pdKernel[nHalfLen+i];
  1241. }
  1242. }
  1243. pdTmp[y*nWidth + x] = dDotMul/dWeightSum ;
  1244. }
  1245. }
  1247. // y方向进行滤波
  1248. for(x=0; x<nWidth; x++)
  1249. {
  1250. for(y=0; y<nHeight; y++)
  1251. {
  1252. dDotMul = 0;
  1253. dWeightSum = 0;
  1254. for(i=(-nHalfLen); i<=nHalfLen; i++)
  1255. {
  1256. // 判断是否在图象内部
  1257. if( (i+y) >= 0  && (i+y) < nHeight )
  1258. {
  1259. dDotMul += (double)pdTmp[(y+i)*nWidth + x] * pdKernel[nHalfLen+i];
  1260. dWeightSum += pdKernel[nHalfLen+i];
  1261. }
  1262. }
  1263. pUnchSmthdImg[y*nWidth + x] = (unsigned char)(int)dDotMul/dWeightSum ;
  1264. }
  1265. }
  1267. // 释放内存
  1268. delete []pdKernel;
  1269. pdKernel = NULL ;
  1271. delete []pdTmp;
  1272. pdTmp = NULL;
  1273. }
  1276. /*************************************************************************
  1277.  *
  1278.  * 函数名称：
  1279.  *   DirGrad()
  1280.  *
  1281.  * 输入参数:
  1282.  *   unsigned char *pUnchSmthdImg - 经过高斯滤波后的图象
  1283.  *   int nWidht - 图象宽度
  1284.  *   int nHeight       - 图象高度
  1285.  *   int *pnGradX                         - x方向的方向导数
  1286.  *   int *pnGradY                         - y方向的方向导数
  1287.  * 返回值:
  1288.  *   无
  1289.  *
  1290.  * 说明:
  1291.  *   这个函数计算方向倒数，采用的微分算子是(-1 0 1) 和 (-1 0 1)'(转置)
  1292.  *   计算的时候对边界象素采用了特殊处理
  1293.  *   
  1294.  *   
  1295.  *************************************************************************
  1296.  */
  1297. void DirGrad(unsigned char *pUnchSmthdImg, int nWidth, int nHeight,
  1298.  int *pnGradX , int *pnGradY)
  1299. {
  1300. // 循环控制变量
  1301. int y ;
  1302. int x ;
  1304. // 计算x方向的方向导数，在边界出进行了处理，防止要访问的象素出界
  1305. for(y=0; y<nHeight; y++)
  1306. {
  1307. for(x=0; x<nWidth; x++)
  1308. {
  1309. pnGradX[y*nWidth+x] = (int) ( pUnchSmthdImg[y*nWidth+min(nWidth-1,x+1)]  
  1310.           - pUnchSmthdImg[y*nWidth+max(0,x-1)]     );
  1311. }
  1312. }
  1314. // 计算y方向的方向导数，在边界出进行了处理，防止要访问的象素出界
  1315. for(x=0; x<nWidth; x++)
  1316. {
  1317. for(y=0; y<nHeight; y++)
  1318. {
  1319. pnGradY[y*nWidth+x] = (int) ( pUnchSmthdImg[min(nHeight-1,y+1)*nWidth + x]  
  1320. - pUnchSmthdImg[max(0,y-1)*nWidth+ x ]     );
  1321. }
  1322. }
  1323. }
  1324. /*************************************************************************
  1325.  *
  1326.  * 函数名称：
  1327.  *   GradMagnitude()
  1328.  *
  1329.  * 输入参数:
  1330.  *   int *pnGradX                         - x方向的方向导数
  1331.  *   int *pnGradY                         - y方向的方向导数
  1332.  *   int nWidht - 图象宽度
  1333.  *   int nHeight       - 图象高度
  1334.  *   int *pnMag                           - 梯度幅度   
  1335.  *
  1336.  * 返回值:
  1337.  *   无
  1338.  *
  1339.  * 说明:
  1340.  *   这个函数利用方向倒数计算梯度幅度，方向倒数是DirGrad函数计算的结果
  1341.  *   
  1342.  *************************************************************************
  1343.  */
  1344. void GradMagnitude(int *pnGradX, int *pnGradY, int nWidth, int nHeight, int *pnMag)
  1345. {
  1347. // 循环控制变量
  1348. int y ;
  1349. int x ;
  1351. // 中间变量
  1352. double dSqtOne;
  1353. double dSqtTwo;
  1355. for(y=0; y<nHeight; y++)
  1356. {
  1357. for(x=0; x<nWidth; x++)
  1358. {
  1359. dSqtOne = pnGradX[y*nWidth + x] * pnGradX[y*nWidth + x];
  1360. dSqtTwo = pnGradY[y*nWidth + x] * pnGradY[y*nWidth + x];
  1361. pnMag[y*nWidth + x] = (int)(sqrt(dSqtOne + dSqtTwo) + 0.5);
  1362. }
  1363. }
  1364. }
  1365. /*************************************************************************
  1366.  *
  1367.  * 函数名称：
  1368.  *   NonmaxSuppress()
  1369.  *
  1370.  * 输入参数:
  1371.  *   int *pnMag                - 梯度图
  1372.  *   int *pnGradX  - x方向的方向导数
  1373.  *   int *pnGradY              - y方向的方向导数
  1374.  *   int nWidth                - 图象数据宽度
  1375.  *   int nHeight               - 图象数据高度
  1376.  *   unsigned char *pUnchRst   - 经过NonmaxSuppress处理后的结果
  1377.  *
  1378.  * 返回值:
  1379.  *   无
  1380.  *
  1381.  * 说明:
  1382.  *   抑止梯度图中非局部极值点的象素。
  1383.  *   
  1384.  *************************************************************************
  1385.  */
  1386. void NonmaxSuppress(int *pnMag, int *pnGradX, int *pnGradY, int nWidth, 
  1387. int nHeight, unsigned char *pUnchRst)
  1388. {
  1389. // 循环控制变量
  1390. int y ;
  1391. int x ;
  1392. int nPos;
  1394. // x方向梯度分量
  1395. int gx  ;
  1396. int gy  ;
  1398. // 临时变量
  1399. int g1, g2, g3, g4 ;
  1400. double weight  ;
  1401. double dTmp1   ;
  1402. double dTmp2   ;
  1403. double dTmp    ;
  1405. // 设置图象边缘部分为不可能的边界点
  1406. for(x=0; x<nWidth; x++)
  1407. {
  1408. pUnchRst[x] = 0 ;
  1409. pUnchRst[nHeight-1+x] = 0;
  1410. }
  1411. for(y=0; y<nHeight; y++)
  1412. {
  1413. pUnchRst[y*nWidth] = 0 ;
  1414. pUnchRst[y*nWidth + nWidth-1] = 0;
  1415. }
  1417. for(y=1; y<nHeight-1; y++)
  1418. {
  1419. for(x=1; x<nWidth-1; x++)
  1420. {
  1421. nPos = y*nWidth + x ;
  1423. // 如果当前象素的梯度幅度为0，则不是边界点
  1424. if(pnMag[nPos] == 0 )
  1425. {
  1426. pUnchRst[nPos] = 0 ;
  1427. }
  1428. else
  1429. {
  1430. // 当前象素的梯度幅度
  1431. dTmp = pnMag[nPos] ;
  1433. // x，y方向导数
  1434. gx = pnGradX[nPos]  ;
  1435. gy = pnGradY[nPos]  ;
  1437. // 如果方向导数y分量比x分量大，说明导数的方向更加“趋向”于y分量。
  1438. if (abs(gy) > abs(gx)) 
  1439. {
  1440. // 计算插值的比例
  1441. weight = fabs(gx)/fabs(gy); 
  1443. g2 = pnMag[nPos-nWidth] ; 
  1444. g4 = pnMag[nPos+nWidth] ;
  1446. // 如果x，y两个方向的方向导数的符号相同
  1447. // C是当前象素，与g1-g4的位置关系为：
  1448. // g1 g2 
  1449. //  C         
  1450. //  g4 g3 
  1451. if (gx*gy > 0) 
  1452. { 
  1453. g1 = pnMag[nPos-nWidth-1] ;
  1454. g3 = pnMag[nPos+nWidth+1] ;
  1455. } 
  1457. // 如果x，y两个方向的方向导数的符号相反
  1458. // C是当前象素，与g1-g4的位置关系为：
  1459. //    g2 g1
  1460. //  C         
  1461. // g3 g4  
  1462. else 
  1463. { 
  1464. g1 = pnMag[nPos-nWidth+1] ;
  1465. g3 = pnMag[nPos+nWidth-1] ;
  1466. } 
  1467. }
  1469. // 如果方向导数x分量比y分量大，说明导数的方向更加“趋向”于x分量
  1470. // 这个判断语句包含了x分量和y分量相等的情况
  1471. else
  1472. {
  1473. // 计算插值的比例
  1474. weight = fabs(gy)/fabs(gx); 
  1476. g2 = pnMag[nPos+1] ; 
  1477. g4 = pnMag[nPos-1] ;
  1479. // 如果x，y两个方向的方向导数的符号相同
  1480. // C是当前象素，与g1-g4的位置关系为：
  1481. // g3   
  1482. // g4 C g2       
  1483. //       g1
  1484. if (gx*gy > 0) 
  1485. {
  1486. g1 = pnMag[nPos+nWidth+1] ;
  1487. g3 = pnMag[nPos-nWidth-1] ;
  1488. } 
  1489. // 如果x，y两个方向的方向导数的符号相反
  1490. // C是当前象素，与g1-g4的位置关系为：
  1491. //      g1
  1492. // g4 C g2       
  1493. //  g3     
  1494. else 
  1495. { 
  1496. g1 = pnMag[nPos-nWidth+1] ;
  1497. g3 = pnMag[nPos+nWidth-1] ;
  1498. }
  1499. }
  1501. // 下面利用g1-g4对梯度进行插值
  1502. {
  1503. dTmp1 = weight*g1 + (1-weight)*g2 ;
  1504. dTmp2 = weight*g3 + (1-weight)*g4 ;
  1506. // 当前象素的梯度是局部的最大值
  1507. // 该点可能是个边界点
  1508. if(dTmp>=dTmp1 && dTmp>=dTmp2)
  1509. {
  1510. pUnchRst[nPos] = 128 ;
  1511. }
  1512. else
  1513. {
  1514. // 不可能是边界点
  1515. pUnchRst[nPos] = 0 ;
  1516. }
  1517. }
  1518. } //else
  1519. } // for
  1521. }
  1522. } 
  1524. /*************************************************************************
  1525.  *
  1526.  * 函数名称：
  1527.  *   TraceEdge()
  1528.  *
  1529.  * 输入参数:
  1530.  *   int    x - 跟踪起点的x坐标 
  1531.  *   int    y - 跟踪起点的y坐标
  1532.  *   int nLowThd - 判断一个点是否为边界点的低阈值
  1533.  *   unsigned char *pUnchEdge - 记录边界点的缓冲区
  1534.  *   int *pnMag               - 梯度幅度图
  1535.  *   int nWidth               - 图象数据宽度
  1536.  *
  1537.  * 返回值:
  1538.  *   无
  1539.  *
  1540.  * 说明:
  1541.  *   递归调用  
  1542.  *   从(x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pUnchEdge中没有处理并且
  1543.  *   可能是边界点的那些象素(=128),象素值为0表明该点不可能是边界点，象素值
  1544.  *   为255表明该点已经被设置为边界点，不必再考虑
  1545.  *   
  1546.  *   
  1547.  *************************************************************************
  1548.  */
  1549. void TraceEdge (int y, int x, int nLowThd, unsigned char *pUnchEdge, int *pnMag, int nWidth) 
  1550. { 
  1551. // 对8邻域象素进行查询
  1552. int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;
  1553. int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;
  1555. int yy ;
  1556. int xx ;
  1558. int k  ;
  1560. for(k=0; k<8; k++)
  1561. {
  1562. yy = y + yNb[k] ;
  1563. xx = x + xNb[k] ;
  1564. // 如果该象素为可能的边界点，又没有处理过
  1565. // 并且梯度大于阈值
  1566. if(pUnchEdge[yy*nWidth+xx] == 128  && pnMag[yy*nWidth+xx]>=nLowThd)
  1567. {
  1568. // 把该点设置成为边界点
  1569. pUnchEdge[yy*nWidth+xx] = 255 ;
  1571. // 以该点为中心进行跟踪
  1572. TraceEdge(yy, xx, nLowThd, pUnchEdge, pnMag, nWidth);
  1573. }
  1574. }
  1575. } 
  1577. /*************************************************************************
  1578.  *
  1579.  * 函数名称：
  1580.  *   EstimateThreshold()
  1581.  *
  1582.  * 输入参数:
  1583.  *   int *pnMag               - 梯度幅度图
  1584.  *  int nWidth               - 图象数据宽度
  1585.  *  int nHeight              - 图象数据高度
  1586.  *   int *pnThdHigh           - 高阈值
  1587.  *   int *pnThdLow            - 低阈值
  1588.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1589.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1590.  *   unsigned char *pUnchEdge - 经过non-maximum处理后的数据
  1591.  *
  1592.  * 返回值:
  1593.  *   无
  1594.  *
  1595.  * 说明:
  1596.  *   经过non-maximum处理后的数据pUnchEdge，统计pnMag的直方图，确定阈值。
  1597.  *   本函数中只是统计pUnchEdge中可能为边界点的那些象素。然后利用直方图，
  1598.  *   根据dRatioHigh设置高阈值，存储到pnThdHigh。利用dRationLow和高阈值，
  1599.  *   设置低阈值，存储到*pnThdLow。dRatioHigh是一种比例：表明梯度小于
  1600.  *   *pnThdHigh的象素数目占象素总数目的比例。dRationLow表明*pnThdHigh
  1601.  *   和*pnThdLow的比例，这个比例在canny算法的原文里，作者给出了一个区间。
  1602.  *
  1603.  *************************************************************************
  1604.  */
  1605. void EstimateThreshold(int *pnMag, int nWidth, int nHeight, int *pnThdHigh,int *pnThdLow, 
  1606.  unsigned char * pUnchEdge, double dRatioHigh, double dRationLow) 
  1607. { 
  1608. // 循环控制变量
  1609. int y;
  1610. int x;
  1611. int k;
  1613. // 该数组的大小和梯度值的范围有关，如果采用本程序的算法，那么梯度的范围不会超过pow(2,10)
  1614. int nHist[1024] ;
  1616. // 可能的边界数目
  1617. int nEdgeNb     ;
  1619. // 最大梯度值
  1620. int nMaxMag     ;
  1622. int nHighCount  ;
  1624. nMaxMag = 0     ; 
  1626. // 初始化
  1627. for(k=0; k<1024; k++) 
  1628. {
  1629. nHist[k] = 0; 
  1630. }
  1632. // 统计直方图，然后利用直方图计算阈值
  1633. for(y=0; y<nHeight; y++)
  1634. {
  1635. for(x=0; x<nWidth; x++)
  1636. {
  1637. // 只是统计那些可能是边界点，并且还没有处理过的象素
  1638. if(pUnchEdge[y*nWidth+x]==128)
  1639. {
  1640. nHist[ pnMag[y*nWidth+x] ]++;
  1641. }
  1642. }
  1643. }
  1645. nEdgeNb = nHist[0]  ;
  1646. nMaxMag = 0         ;
  1647. // 统计经过“非最大值抑止(non-maximum suppression)”后有多少象素
  1648. for(k=1; k<1024; k++)
  1649. {
  1650. if(nHist[k] != 0)
  1651. {
  1652. // 最大梯度值
  1653. nMaxMag = k;
  1654. }
  1656. // 梯度为0的点是不可能为边界点的
  1657. // 经过non-maximum suppression后有多少象素
  1658. nEdgeNb += nHist[k];
  1659. }
  1661. // 梯度比高阈值*pnThdHigh小的象素点总数目
  1662. nHighCount = (int)(dRatioHigh * nEdgeNb +0.5);
  1664. k = 1;
  1665. nEdgeNb = nHist[1];
  1667. // 计算高阈值
  1668. while( (k<(nMaxMag-1)) && (nEdgeNb < nHighCount) )
  1669. {
  1670. k++;
  1671. nEdgeNb += nHist[k];
  1672. }
  1674. // 设置高阈值
  1675. *pnThdHigh = k ;
  1677. // 设置低阈值
  1678. *pnThdLow  = (int)((*pnThdHigh) * dRationLow+ 0.5);
  1679. }
  1681. /*************************************************************************
  1682.  *
  1683.  * 函数名称：
  1684.  *   Hysteresis()
  1685.  *
  1686.  * 输入参数:
  1687.  *   int *pnMag               - 梯度幅度图
  1688.  *  int nWidth               - 图象数据宽度
  1689.  *  int nHeight              - 图象数据高度
  1690.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1691.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1692.  *   unsigned char *pUnchEdge - 记录边界点的缓冲区
  1693.  *
  1694.  * 返回值:
  1695.  *   无
  1696.  *
  1697.  * 说明:
  1698.  *   本函数实现类似“磁滞现象”的一个功能，也就是，先调用EstimateThreshold
  1699.  *   函数对经过non-maximum处理后的数据pUnchSpr估计一个高阈值，然后判断
  1700.  *   pUnchSpr中可能的边界象素(=128)的梯度是不是大于高阈值nThdHigh，如果比
  1701.  *   该阈值大，该点将作为一个边界的起点，调用TraceEdge函数，把对应该边界
  1702.  *   的所有象素找出来。最后，当整个搜索完毕时，如果还有象素没有被标志成
  1703.  *   边界点，那么就一定不是边界点。
  1704.  *   
  1705.  *************************************************************************
  1706.  */
  1707. void Hysteresis(int *pnMag, int nWidth, int nHeight, double dRatioLow, 
  1708. double dRatioHigh, unsigned char *pUnchEdge)
  1709. {
  1710. // 循环控制变量
  1711. int y;
  1712. int x;
  1714. int nThdHigh ;
  1715. int nThdLow  ;
  1717. int nPos;
  1719. // 估计TraceEdge需要的低阈值，以及Hysteresis函数使用的高阈值
  1720. EstimateThreshold(pnMag, nWidth, nHeight, &nThdHigh, 
  1721.                &nThdLow, pUnchEdge,dRatioHigh, dRatioLow);
  1723.   // 这个循环用来寻找大于nThdHigh的点，这些点被用来当作边界点，然后用
  1724. // TraceEdge函数来跟踪该点对应的边界
  1725.    for(y=0; y<nHeight; y++)
  1726.  {
  1727.       for(x=0; x<nWidth; x++)
  1728. {
  1729. nPos = y*nWidth + x ; 
  1731. // 如果该象素是可能的边界点，并且梯度大于高阈值，该象素作为
  1732. // 一个边界的起点
  1733. if((pUnchEdge[nPos] == 128) && (pnMag[nPos] >= nThdHigh))
  1734. {
  1735. // 设置该点为边界点
  1736. pUnchEdge[nPos] = 255;
  1737. TraceEdge(y, x, nThdLow, pUnchEdge, pnMag, nWidth);
  1738. }
  1739.       }
  1740.    }
  1742.  // 那些还没有被设置为边界点的象素已经不可能成为边界点
  1743.    for(y=0; y<nHeight; y++)
  1744.  {
  1745.  for(x=0; x<nWidth; x++)
  1746.  {
  1747.  nPos = y*nWidth + x ; 
  1748.  if(pUnchEdge[nPos] != 255)
  1749.  {
  1750.  // 设置为非边界点
  1751.  pUnchEdge[nPos] = 0 ;
  1752.  }
  1753.  }
  1754.  }
  1755. }
  1758. /*************************************************************************
  1759.  *
  1760.  * 函数名称：
  1761.  *   Canny()
  1762.  *
  1763.  * 输入参数:
  1764.  *   unsigned char *pUnchImage- 图象数据
  1765.  *  int nWidth               - 图象数据宽度
  1766.  *  int nHeight              - 图象数据高度
  1767.  *   double sigma             - 高斯滤波的标准方差
  1768.  *  double dRatioLow         - 低阈值和高阈值之间的比例
  1769.  *  double dRatioHigh        - 高阈值占图象象素总数的比例
  1770.  *   unsigned char *pUnchEdge - canny算子计算后的分割图
  1771.  *
  1772.  * 返回值:
  1773.  *   无
  1774.  *
  1775.  * 说明:
  1776.  *   canny分割算子，计算的结果保存在pUnchEdge中，逻辑1(255)表示该点为
  1777.  *   边界点，逻辑0(0)表示该点为非边界点。该函数的参数sigma，dRatioLow
  1778.  *   dRatioHigh，是需要指定的。这些参数会影响分割后边界点数目的多少
  1779.  *************************************************************************
  1780.  */
  1781. void Canny(unsigned char *pUnchImage, int nWidth, int nHeight, double sigma,
  1782.  double dRatioLow, double dRatioHigh, unsigned char *pUnchEdge)
  1783. {
  1784. // 经过高斯滤波后的图象数据
  1785. unsigned char * pUnchSmooth ;
  1786.   
  1787. // 指向x方向导数的指针
  1788. int * pnGradX ; 
  1790. // 指向y方向导数的指针
  1791. int * pnGradY ;
  1793. // 梯度的幅度
  1794. int * pnGradMag ;
  1796. pUnchSmooth  = new unsigned char[nWidth*nHeight] ;
  1797. pnGradX      = new int [nWidth*nHeight]          ;
  1798. pnGradY      = new int [nWidth*nHeight]          ;
  1799. pnGradMag    = new int [nWidth*nHeight]          ;
  1801. // 对原图象进行滤波
  1802. GaussianSmooth(pUnchImage, nWidth, nHeight, sigma, pUnchSmooth) ;
  1804. // 计算方向导数
  1805. DirGrad(pUnchSmooth, nWidth, nHeight, pnGradX, pnGradY) ;
  1807. // 计算梯度的幅度
  1808. GradMagnitude(pnGradX, pnGradY, nWidth, nHeight, pnGradMag) ;
  1810. // 应用non-maximum 抑制
  1811. NonmaxSuppress(pnGradMag, pnGradX, pnGradY, nWidth, nHeight, pUnchEdge) ;
  1813. // 应用Hysteresis，找到所有的边界
  1814. Hysteresis(pnGradMag, nWidth, nHeight, dRatioLow, dRatioHigh, pUnchEdge);
  1817. // 释放内存
  1818. delete pnGradX      ;
  1819. pnGradX      = NULL ;
  1820. delete pnGradY      ;
  1821. pnGradY      = NULL ;
  1822. delete pnGradMag    ;
  1823. pnGradMag    = NULL ;
  1824. delete pUnchSmooth ;
  1825. pUnchSmooth  = NULL ;
  1826. }