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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 图片显示 > 查看源码
MathEx.cpp
资源名称：QuickImage.rar　[点击查看]
上传用户：gzboli
上传日期：2013-04-10
资源大小：471k
文件大小：14k
源码类别：

图片显示

开发平台：

Visual C++

MathEx.cpp：源码内容
  1. // MathEx.cpp: implementation of the CMathEx class.
  2. //
  3. //////////////////////////////////////////////////////////////////////
  5. #include "stdafx.h"
  6. #include "QuickImage.h"
  7. #include "MathEx.h"
  8. #include <math.h>
  10. #ifdef _DEBUG
  11. #undef THIS_FILE
  12. static char THIS_FILE[]=__FILE__;
  13. #define new DEBUG_NEW
  14. #endif
  16. //////////////////////////////////////////////////////////////////////
  17. // Construction/Destruction
  18. //////////////////////////////////////////////////////////////////////
  19. /*
  20. COMPLEX::COMPLEX()
  21. {
  22. re = 0.0;
  23. im = 0.0;
  24. }
  26. COMPLEX::COMPLEX(double real, double image)
  27. {
  28. re = real;
  29. im = image;
  30. }
  32. COMPLEX::COMPLEX(const COMPLEX &c)
  33. {
  34. re = c.re;
  35. im = c.im;
  36. }
  38. COMPLEX::~COMPLEX()
  39. {
  40. }
  42. COMPLEX COMPLEX::operator +(const COMPLEX &c)
  43. {
  44. COMPLEX ret;
  45. ret.re = re + c.re;
  46. ret.im = im + c.im;
  47. return ret;
  48. }
  51. COMPLEX COMPLEX::operator -(const COMPLEX &c)
  52. {
  53. COMPLEX ret;
  54. ret.re = re - c.re;
  55. ret.im = im - c.im;
  56. return ret;
  57. }
  59. COMPLEX COMPLEX::operator *(const COMPLEX &c)
  60. {
  61. COMPLEX ret;
  62. ret.re = re * c.re - im * c.im;
  63. ret.im = re * c.im + im * c.re;
  64. return ret;
  65. }
  67. COMPLEX& COMPLEX::operator =(const COMPLEX &c)
  68. {
  69. re = c.re;
  70. im = c.im;
  71. return *this;
  72. }
  74. BOOL COMPLEX::operator ==(const COMPLEX &c)
  75. {
  76. return ((re == c.re ) && (im == c.im));
  77. }
  79. COMPLEX COMPLEX::operator *(double d)
  80. {
  81. COMPLEX ret;
  82. ret.re = re * d;
  83. ret.im = im * d;
  84. return ret;
  85. }
  87. COMPLEX operator *(double d, const COMPLEX &c)
  88. {
  89. COMPLEX ret;
  90. ret.re = c.re * d;
  91. ret.im = c.im * d;
  92. return ret;
  93. }
  94. */
  95. /*complex add*/
  96. COMPLEX Add(COMPLEX c1, COMPLEX c2)
  97. {
  98. COMPLEX c;
  99. c.re=c1.re+c2.re;
  100. c.im=c1.im+c2.im;
  101. return c;
  102. }
  104. /*complex substract*/
  105. COMPLEX Sub(COMPLEX c1, COMPLEX c2)
  106. {
  107. COMPLEX c;
  108. c.re=c1.re-c2.re;
  109. c.im=c1.im-c2.im;
  110. return c;
  111. }
  113. /*complex multiple*/
  114. COMPLEX Mul(COMPLEX c1, COMPLEX c2)
  115. {
  116. COMPLEX c;
  117. c.re=c1.re*c2.re-c1.im*c2.im;
  118. c.im=c1.re*c2.im+c2.re*c1.im;
  119. return c;
  120. }
  122. CMathEx::CMathEx()
  123. {
  125. }
  127. CMathEx::~CMathEx()
  128. {
  130. }
  132. void CMathEx::MatrixRotate(double *after, const double *before, int row, int col)
  133. {
  134. ASSERT( NULL != after);
  135. ASSERT( NULL != before);
  137. for(int i=0; i<col; i++)
  138. for(int j=0; j<row; j++)
  139. after[i * row + j] = before[j * col + i];
  141. }
  143. void CMathEx::MatrixMutiply
  144. (double *result, const double *left, const double *right, int row, int coll, int colr)
  145. {
  146. ASSERT(NULL != result);
  147. ASSERT(NULL != left);
  148. ASSERT(NULL != right);
  149. int i,j,k;
  150. double *pRes = result;
  152. for(i = 0; i < row * colr; i++, pRes++)
  153. *pRes = 0.0;
  155. for(i = 0; i < row; i++)
  156. for(j = 0; j< colr;j++)
  157. for(k = 0;k< coll ; k++)
  158. result[colr*i+j] += (left[coll*i+k] * right[colr*k+j]);
  159. }
  161. BOOL CMathEx::MatrixConvert(double *result, const double *before, int size)
  162. {
  163. ASSERT(NULL != result);
  164. ASSERT(NULL != before);
  166. int i = size * size * sizeof(double);
  167. HANDLE hHeap = HeapCreate(HEAP_NO_SERIALIZE, i, 0);
  168. if(hHeap == NULL)
  169. {
  170. return FALSE;
  171. }
  172. HeapLock(hHeap);
  173. double *pSubthis = NULL;
  174. pSubthis=(double*)HeapAlloc(hHeap, HEAP_NO_SERIALIZE, i);
  175. if(pSubthis == NULL)
  176. {
  177. HeapUnlock(hHeap);
  178. HeapDestroy(hHeap);
  179. return FALSE;
  180. }
  181. memcpy((void*)pSubthis, (void*)before, i);
  183. double *pRes = result;
  184. for(i=0; i<size * size; i++, pRes++)
  185. *pRes = 1.0;
  187. int k, k2;
  188. int imax;
  189. double dmax;
  190. for(i=0;i<size;i++)
  191. {
  192. //////////////////////////////////////////////////////////////////
  193. imax=i;
  194. dmax=pSubthis[i*size+i];
  195. for(k=i+1;k<size;k++)
  196. if(fabs(pSubthis[k*size+i])>fabs(dmax))
  197. {
  198. dmax=pSubthis[k*size+i];
  199. imax=k;
  200. }
  201. if(fabs(dmax)<1e-30)
  202. {
  203. AfxMessageBox("Can not convert a ill matrix!",MB_ICONERROR);
  204. HeapFree(hHeap,HEAP_NO_SERIALIZE | HEAP_GENERATE_EXCEPTIONS,pSubthis);
  205. pSubthis = NULL;
  206. HeapUnlock(hHeap);
  207. HeapDestroy(hHeap);
  208. return FALSE;
  209. }
  210. if(imax!=i)
  211. for(k=0;k<size;k++)
  212. {
  213. dmax=pSubthis[i*size+k];
  214. pSubthis[i*size+k]=pSubthis[imax*size+k];
  215. pSubthis[imax*size+k]=dmax;
  217. dmax = result[i*size+k];
  218. result[i*size+k] = result[imax*size+k];
  219. result[imax*size+k] = dmax;
  220. }
  221. //此间为列选主元
  222. ////////////////////////////////////////////////////////////////////
  223. dmax= pSubthis[i*size+i];
  224. for(k=i;k<size;k++)
  225. {
  226. imax = i*size+k;
  227. pSubthis[imax] /= dmax;
  228. result[imax] /= dmax;
  229. }
  230. for(k=0; k<i; k++)
  231. {
  232. dmax=pSubthis[k*size+i];
  233. for(k2 = i; k2 < size; k2++)
  234. {
  235. pSubthis[k*size+k2] -= (pSubthis[i*size+k2]*dmax);
  236. result[k*size+k2] -= (result[i*size+k2]*dmax);
  237. }
  238. }
  239. for(k=i+1; k<size; k++)
  240. {
  241. dmax=pSubthis[k*size+i];
  242. for(k2 = i; k2 < size; k2++)
  243. {
  244. pSubthis[k*size+k2] -= (pSubthis[i*size+k2]*dmax);
  245. result[k*size+k2] -= (result[i*size+k2]*dmax);
  246. }
  247. }
  248. }
  250. try
  251. {
  252. HeapFree(hHeap,HEAP_NO_SERIALIZE | HEAP_GENERATE_EXCEPTIONS,pSubthis);
  253. pSubthis = NULL;
  254. HeapUnlock(hHeap);
  255. HeapDestroy(hHeap);
  256. }
  257. catch(CMemoryException *e)
  258. {
  259. char msg[256];
  260. e->GetErrorMessage(msg,255);
  261. e->ReportError();
  262. e->Delete();
  263. return FALSE;
  264. }
  266. return TRUE;
  267. }
  269. void CMathEx::JieFC(double *xsz, int row, int col, double *result)
  270. {
  271. double b,ab1;
  272. int i,j,k,j1,n1,n2,l;
  273. for( i=0; i<row; i++)
  274. {
  275. result[i]=0.0;
  276. b=0.0;
  277. for(j=i; j<row; j++)
  278. {
  279. if(fabs(b)<=fabs(xsz[j * col + i]))
  280. {
  281. b=xsz[j * col + i];
  282. j1=j;
  283. }
  284. }
  285. for(k=i; k<col; k++)
  286. {
  287. ab1 = xsz[j1 * col + k]/b;
  288. xsz[j1 * col + k] = xsz[i * col + k];
  289. xsz[i * col + k] = ab1;
  290. }
  291. n1=i+1;
  292. if(n1<row)
  293. {
  294. for(j=n1; j<row;j++)
  295. {
  296. for(k=n1;k<col;k++)
  297. {
  298. xsz[j * col + k] -= (xsz[j * col + i]*xsz[i * col + k]);
  299. }
  300. }
  301. }
  302. }
  303. result[row-1] = xsz[row * col -1];//xsz[(row-1) * col + col-1];
  304. n2=row-1;
  305. for(i=0;i<n2;i++)
  306. {
  307. l=n2- 1 - i;
  308. result[l]=xsz[l * col -1];
  309. for(j=l + 1; j<row; j++)
  310. {
  311. result[l] -= result[j]*xsz[l * col + j];
  312. }
  313. }
  314. }
  316. BOOL CMathEx::GSXQ(double *x, double *a, double *b, int size)
  317. {
  318. int i, j, k;
  319. double dMax;
  320. int iMax;
  321. BOOL bIll = TRUE;
  322. for(i = 0; i< size; i++)
  323. {
  324. iMax = i;
  325. dMax = x[i * size + i];
  326. for(j = i +1; j< size; j++)
  327. {
  328. if(fabs(dMax) < fabs(a[j * size + i]))
  329. {
  330. dMax = a[j * size + i];
  331. iMax = j;
  332. }
  333. }
  334. if(dMax > -0.00000000001 && dMax < 0.00000000001)
  335. bIll = FALSE;
  336. if(iMax != i)
  337. {
  338. for(j = i; j < size; j++)
  339. {
  340. dMax = a[i * size + j];
  341. a[i * size + j] = a[iMax * size + j];
  342. a[iMax * size + j] = dMax;
  343. }
  344. dMax = b[i];
  345. b[i] = b[iMax];
  346. b[iMax] = dMax;
  347. }
  348. for(j = i +1; j < size; j++)
  349. {
  350. dMax = a[j * size + i] / a[i * size + i];
  351. a[j * size + i] = 0.0;
  352. for(k = i + 1; k< size; k++)
  353. {
  354. a[j * size + k] -= dMax * a[i * size + k];
  355. }
  356. b[j] -= dMax * b[i];
  357. }
  359. }
  360. x[size -1] = b[size -1] / a[size * size -1];
  361. for(k = size - 2; k > -1; k--)
  362. {
  363. dMax = 0.0;
  364. for(j = k + 1; j< size; j++)
  365. {
  366. dMax += a[k * size + j] * x[j];
  367. }
  368. x[k] = (b[k] - dMax) / a[k * size + k];
  369. }
  370. return bIll;
  371. }
  373. BOOL CMathEx::FFT(const COMPLEX *TD, COMPLEX *FD, int power)
  374. {
  375. ASSERT((NULL != TD) && (NULL != FD));
  376. int i, j, k, bfsize, p, iIndex;;
  377. double angle;
  378. int count = 1 << power;//变换点数
  379. double PAI2 = PAI * 2.0;
  380. COMPLEX *W, *X1, *X2, *X = NULL;
  381. W = new COMPLEX[count / 2];
  382. if(NULL == W)
  383. {
  384. printf("Sorry, insufficent memory available!");
  385. return FALSE;
  386. }
  387. X1 = new COMPLEX[count];
  388. if(NULL == X1)
  389. {
  390. delete W;
  391. printf("Sorry, insufficent memory available!");
  392. return FALSE;
  393. }
  394. X2 = new COMPLEX[count];
  395. if(NULL == X2)
  396. {
  397. printf("Sorry, insufficent memory available!");
  398. delete X1;
  399. delete W;
  400. return FALSE;
  401. }
  402. //计算加权系数
  403. for(i = 0; i < count / 2 ; i++)
  404. {
  405. angle = -i * PAI2 / count;
  406. W[i].re = cos(angle);
  407. W[i].im = sin(angle);
  408. }
  409. memcpy(X1, TD, sizeof(COMPLEX) * count);
  410. //蝶形运算
  411. for(k = 0; k < power; k++)
  412. {
  413. for(j = 0; j < 1<<k; j++)
  414. {
  415. bfsize = (1<<(power - k)) / 2;
  416. for(i = 0; i < bfsize; i++)
  417. {
  418. p = j * bfsize * 2;
  419. iIndex = i + p;
  420. X2[iIndex] = Add(X1[iIndex], X1[iIndex + bfsize]);
  421. X2[iIndex + bfsize] = Mul(Sub(X1[iIndex], X1[iIndex + bfsize]),
  422. (W[i * (1 << k)]));
  423. }
  424. }
  425. X = X1;
  426. X1 = X2;
  427. X2 = X;
  428. }
  429. //重新排续
  430. for(j =0 ; j < count; j++)
  431. {
  432. p = 0;
  433. for(i = 0; i < power; i++)
  434. {
  435. if(j & (1<<i))
  436. {
  437. p += 1<<(power - i -1);
  438. }
  439. FD[j] = X1[p];
  440. }
  441. }
  442. delete W;
  443. delete X1;
  444. delete X2;
  445. return TRUE;
  446. }
  448. BOOL CMathEx::IFFT(const COMPLEX *FD, COMPLEX *TD, int power)
  449. {
  450. ASSERT((NULL != TD) && (NULL != FD));
  451. int count = 1 << power;
  452. COMPLEX *x = new COMPLEX[count];
  453. if(NULL == x)
  454. {
  455. printf("Sorry, insufficent memory available!");
  456. return FALSE;
  457. }
  458. memcpy(x, FD, sizeof(COMPLEX) * count);
  459. for(int i =0; i < count; i ++)
  460. {
  461. x[i].im = - x[i].im;
  462. }
  463. FFT(x, TD, power);
  465. for(i =0; i < count; i++)
  466. {
  467. TD[i].re /= count;
  468. TD[i].im = -TD[i].im / count;
  469. }
  470. delete x;
  471. return TRUE;
  472. }
  475. /*******************************************************
  476. DCT()
  478. 参数：
  480. f为时域值
  481. F为频域值
  482. power为2的幂数
  484. 返回值：
  486. 无
  488. 说明：
  490. 本函数利用快速傅立叶变换实现快速离散余弦变换
  491. ********************************************************/
  492. void CMathEx::DCT(double *f, double *F, int power)
  493. {
  494. int i,count;
  495. COMPLEX *X;
  496. double s;
  498. /*计算离散余弦变换点数*/
  499. count=1<<power;
  501. /*分配运算所需存储器*/
  502. X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
  504. /*延拓*/
  505. memset(X,0,sizeof(COMPLEX)*count*2);
  507. /*将时域点写入存储器*/
  508. for(i=0;i<count;i++)
  509. {
  510. X[i].re=f[i];
  511. }
  513. /*调用快速傅立叶变换*/
  514. FFT(X,X,power+1);
  516. /*调整系数*/
  517. s=1/sqrt(count);
  518. F[0]=X[0].re*s;
  519. s*=sqrt(2);
  520. for(i=1;i<count;i++)
  521. {
  522. F[i]=(X[i].re*cos(i*PAI/(count*2))+X[i].im*sin(i*PAI/(count*2)))*s;
  523. }
  525. /*释放存储器*/
  526. free(X);
  527. }
  529. /************************************************************
  530. IDCT()
  532. 参数：
  534. F为频域值
  535. f为时域值
  536. power为2的幂数
  538. 返回值：
  540. 无
  542. 说明：
  544. 本函数利用快速傅立叶反变换实现快速离散反余弦变换
  545. *************************************************************/
  546. void CMathEx::IDCT(double *F, double *f, int power)
  547. {
  548. int i,count;
  549. COMPLEX *X;
  550. double s;
  552. /*计算离散反余弦变换点数*/
  553. count=1<<power;
  555. /*分配运算所需存储器*/
  556. X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2);
  558. /*延拓*/
  559. memset(X,0,sizeof(COMPLEX)*count*2);
  561. /*调整频域点，写入存储器*/
  562. for(i=0;i<count;i++)
  563. {
  564. X[i].re=F[i]*cos(i*PAI/(count*2));
  565. X[i].im=F[i]*sin(i*PAI/(count*2));
  566. }
  568. /*调用快速傅立叶反变换*/
  569. IFFT(X,X,power+1);
  571. /*调整系数*/
  572. s=1/sqrt(count);
  573. for(i=1;i<count;i++)
  574. {
  575. f[i]=(1-sqrt(2))*s*F[0]+sqrt(2)*s*X[i].re*count*2;
  576. }
  578. /*释放存储器*/
  579. free(X);
  580. }
  582. /**********************************************************
  583. WALh()
  585. 参数：
  587. f为时域值
  588. W为频域值
  589. power为2的幂数
  591. 返回值：
  593. 无
  595. 说明：
  597. 本函数利用快速傅立叶变换实现快速沃尔什-哈达玛变换
  598. *************************************************************/
  599. void CMathEx::WALh(double *f, double *W, int power)
  600. {
  601. int count;
  602. int i,j,k,bfsize,p;
  603. double *X1,*X2,*X;
  605. /*计算快速沃尔什变换点数*/
  607. count=1<<power;
  609. /*分配运算所需存储器*/
  610. X1=(double *)malloc(sizeof(double)*count);
  611. X2=(double *)malloc(sizeof(double)*count);
  613. /*将时域点写入存储器*/
  614. memcpy(X1,f,sizeof(double)*count);
  616. /*蝶形运算*/
  617. for(k=0;k<power;k++)
  618. {
  619. for(j=0;j<1<<k;j++)
  620. {
  621. bfsize=1<<(power-k);
  622. for(i=0;i<bfsize/2;i++)
  623. {
  624. p=j*bfsize;
  625. X2[i+p]=X1[i+p]+X1[i+p+bfsize/2];
  626. X2[i+p+bfsize/2]=X1[i+p]-X1[i+p+bfsize/2];
  627. }
  628. }
  629. X=X1;
  630. X1=X2;
  631. X2=X;
  632. }
  633. /*调整系数*/
  634. // for(i=0;i<count;i++)
  635. // {
  636. // W[i]=X1[i]/count;
  637. // }
  638. for(j=0;j<count;j++)
  639. {
  640. p=0;
  641. for(i=0;i<power;i++)
  642. {
  643. if (j&(1<<i)) p+=1<<(power-i-1);
  644. }
  645. W[j]=X1[p]/count;
  646. }
  648. /*释放存储器*/
  649. free(X1);
  650. free(X2);
  651. }
  653. /*********************************************************************
  654. IWALh()
  656. 参数：
  658. W为频域值
  659. f为时域值
  660. power为2的幂数
  662. 返回值：
  664. 无
  666. 说明：
  668. 本函数利用快速沃尔什-哈达玛变换实现快速沃尔什-哈达玛反变换
  669. **********************************************************************/
  670. void CMathEx::IWALh(double *W, double *f, int power)
  671. {
  672. int i, count;
  674. /*计算快速沃尔什反变换点数*/
  675. count=1<<power;
  677. /*调用快速沃尔什-哈达玛变换*/
  678. WALh(W, f, power);
  680. /*调整系数*/
  681. for(i=0;i<count;i++)
  682. {
  683. f[i] *= count;
  684. }
  685. }
  687. double CMathEx::GetArea(POINT vert, POINT from, POINT to)
  688. {
  689. int dx, dy;
  690. dx = from.x - vert.x;
  691. dy = from.y - vert.y;
  692. double sf = sqrt(double(dx * dx + dy * dy));
  693. double sinf = (0 == dx && 0 == dy) ? 0 : dy / sf;
  694. double cosf = (0 == dx && 0 == dy) ? 0 : dx / sf;
  695. dx = to.x - vert.x;
  696. dy = to.y - vert.y;
  697. double st = sqrt(double(dx * dx + dy * dy));
  698. double sint = (0 == dx && 0 == dy) ? 0 : dy / st;
  699. double cost = (0 == dx && 0 == dy) ? 0 : dx / st;
  701. return sf * (sinf * cost - cosf * sint) * st * 0.5;
  702. }
  704. int CMathEx::compare(const void *arg1, const void *arg2)
  705. {
  706. if(*((float*)arg1) < *((float*)arg2))
  707. return -1;
  708. if(*((float*)arg1) > *((float*)arg2))
  709. return 1;
  710. return 0;
  711. }
  713. BOOL CMathEx::MedianFilter(float *pData, int iWidth, int iHeight)
  714. {
  715. float *pSub = new float[iWidth * 2];
  716. if(NULL == pSub)
  717. {
  718. return FALSE;
  719. }
  720. float *pWin = new float[9];
  721. if(NULL == pWin)
  722. {
  723. delete pWin;
  724. return FALSE;
  725. }
  726. int x, y;
  727. int iRowSize = iWidth * sizeof(float);
  728. int row;
  730. memcpy(&pSub[iWidth], pData, iRowSize);
  731. for(y = 1; y < iHeight -1; y++)
  732. {
  733. memmove(pSub, &pSub[iWidth], iRowSize);
  734. row = y * iWidth;
  735. memcpy(&pSub[iWidth], &pData[row], iRowSize);
  737. for(x = 1; x <  iWidth-1; x++)
  738. {
  739. pWin[0] = pSub[x -1];
  740. pWin[1] = pSub[x];
  741. pWin[2] = pSub[x +1];
  742. pWin[3] = pSub[iWidth + x -1];
  743. pWin[4] = pSub[iWidth + x];
  744. pWin[5] = pSub[iWidth + x +1];
  745. pWin[6] = pData[row + iWidth + x -1];
  746. pWin[7] = pData[row + iWidth + x];
  747. pWin[8] = pData[row + iWidth + x +1];
  748. qsort(pWin, 9, sizeof(float), compare);
  750. pData[row + x] = pWin[4];
  751. }
  752. }
  754. delete pWin;
  755. delete pSub;
  756. return TRUE;
  757. }