Retinex图像增强算法代码

 

 

http://www.cnblogs.com/sleepwalker/p/3676600.html?utm_source=tuicool

http://blog.csdn.net/carson2005/article/details/9502053

 

 

Retinex理论

Retinex理论始于Land和McCann于20世纪60年代作出的一系列贡献,其基本思想是人感知到某点的颜色和亮度并不仅仅取决于该点进入人眼的绝对光线，还和其周围的颜色和亮度有关。Retinex这个词是由视网膜(Retina)和大脑皮层(Cortex)两个词组合构成的.Land之所以设计这个词，是为了表明他不清楚视觉系统的特性究竟取决于此两个生理结构中的哪一个，抑或是与两者都有关系。

Land的Retinex模型是建立在以下的基础之上的：

一、真实世界是无颜色的，我们所感知的颜色是光与物质的相互作用的结果。我们见到的水是无色的，但是水膜—肥皂膜却是显现五彩缤纷，那是薄膜表面光干涉的结果；

二、每一颜色区域由给定波长的红、绿、蓝三原色构成的；

三、三原色决定了每个单位区域的颜色。

Retinex 理论的基本内容是物体的颜色是由物体对长波（红）、中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定的；物体的色彩不受光照非均性的影响，具有一致性，即Retinex理论是以色感一致性（颜色恒常性）为基础的。如下图所示，观察者所看到的物体的图像S是由物体表面对入射光L反射得到的，反射率R由物体本身决定，不受入射光L变化。

 

 

本来想把下面的代码改成OpenCV版本的，但是不太会把读出来的mat里面数据改成BYTE*里面，在主函数里面写的一点都注释了

 

 

[cpp]  view plain  copy

 print ?

1. // Retinex.cpp : Defines the entry point for the console application.  
2. //  
3.   
4. #include "stdafx.h"  
5.   
6. #include <stdio.h>  
7. #include <string.h>  
8. #include <windows.h>  
9. #include <cmath>  
10. #include <time.h>  
11. #include <iostream>  
12.   
13. #include "opencv2/core/core.hpp"  
14. #include "opencv2/highgui/highgui.hpp"  
15. #include "opencv2/imgproc/imgproc.hpp"  
16.   
17. #pragma comment(lib,"opencv_core2410d.lib")                    
18. #pragma comment(lib,"opencv_highgui2410d.lib")                    
19. #pragma comment(lib,"opencv_imgproc2410d.lib")    
20.   
21.   
22. using namespace std;   
23. using namespace cv;  
24.   
25. #define EPSILON 1  
26. #define DELTA 1  
27. #define GAMMA 0.9     
28. #define PI 3.1415926  
29. #define ALPHA_c 1  
30.   
31. //  
32. //  Read a 8 bit bmp File  
33. //  
34. BYTE *Read8BitBmpFile2Img(const char * filename,int *width,int *height)  
35. {   FILE * BinFile;  
36. BITMAPFILEHEADER FileHeader;  
37. BITMAPINFOHEADER BmpHeader;  
38. BYTE *img;  
39. int size;  
40. int Suc=1;  
41. // Open File  
42. *width=*height=0;  
43. if((BinFile=fopen(filename,"rb"))==NULL) return NULL;  
44. // Read Struct Info  
45. if (fread((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
46. if (fread((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
47. if (Suc==-1) { fclose(BinFile); return NULL; }  
48. // Read Image Data  
49. *width=(BmpHeader.biWidth+3)/4*4;  
50. *height=BmpHeader.biHeight;  
51. size=(BmpHeader.biWidth+3)/4*4*BmpHeader.biHeight;  
52. fseek(BinFile,FileHeader.bfOffBits,SEEK_SET);  
53. if ( (img=new BYTE[size+8])!=NULL)  
54. {   if(fread(img+8-int(img)%8,sizeof(BYTE),size,BinFile)!=(unsigned int)size)  
55. { fclose(BinFile);  
56. delete img;  
57. img=NULL;  
58. return NULL;  
59. }  
60. }  
61. fclose(BinFile);  
62. return img;  
63. }  
64.   
65. /  
66. // Write a 8 bit bmp File  
67. /  
68. int Write8BitImg2BmpFile(BYTE *img,int width,int height,const char * filename)  
69. {   FILE * BinFile;  
70. BITMAPFILEHEADER FileHeader;  
71. BITMAPINFOHEADER BmpHeader;  
72. BYTE p[4];  
73. int i,Suc=1;  
74.   
75. // Open File  
76. if((BinFile=fopen(filename,"w+b"))==NULL) {  return -1; }  
77. //  Fill the FileHeade)  
78. FileHeader.bfType= ((WORD) ('M' << 8) | 'B');  
79. FileHeader.bfOffBits=sizeof(BITMAPFILEHEADER)+sizeof(BmpHeader)+256*4L;  
80. FileHeader.bfSize=FileHeader.bfOffBits+width*height ;  
81. FileHeader.bfReserved1=0;  
82. FileHeader.bfReserved2=0;  
83. if (fwrite((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
84. // Fill the ImgHeader  
85. BmpHeader.biSize = 40;  
86. BmpHeader.biWidth = width;  
87. BmpHeader.biHeight = height;  
88. BmpHeader.biPlanes = 1 ;  
89. BmpHeader.biBitCount = 8 ;  
90. BmpHeader.biCompression = 0 ;  
91. BmpHeader.biSizeImage = 0 ;  
92. BmpHeader.biXPelsPerMeter = 0;  
93. BmpHeader.biYPelsPerMeter = 0;  
94. BmpHeader.biClrUsed = 0;  
95. BmpHeader.biClrImportant = 0;  
96. if (fwrite((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
97. // write Pallete  
98. for (i=0,p[3]=0;i<256;i++)   
99. {    
100.     p[3]=0;  
101.     p[0]=p[1]=p[2]=i; // blue,green,red;  
102.     if (fwrite((void *)p,1,4,BinFile)!=4) { Suc=-1; break; }  
103. }  
104. // write image data  
105. if (fwrite((void *)img,1,width*height,BinFile)!=(unsigned int) width*height) Suc=-1;  
106. // return;  
107. fclose(BinFile);  
108. return Suc;  
109. }  
110.   
111. //  
112. //  Read a 24 bit bmp File  
113. //  
114. BYTE *Read24BitBmpFile2Img(const char * filename,int *width,int *height)  
115. {    
116.     FILE * BinFile;  
117.     BITMAPFILEHEADER FileHeader;  
118.     BITMAPINFOHEADER BmpHeader;  
119.     BYTE *img;  
120.     int size;  
121.     int Suc=1;  
122.   
123.     // Open File  
124.     *width=*height=0;  
125.     if((BinFile=fopen(filename,"rb"))==NULL) return NULL;  
126.     // Read Struct Info  
127.     if (fread((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
128.     if (fread((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
129.     if (Suc==-1) { fclose(BinFile); return NULL; }  
130.     // Read Image Data  
131.     *width=(BmpHeader.biWidth+3)/4*4;  
132.     *height=BmpHeader.biHeight;  
133.     size=(*width)*(*height)*3;  
134.     fseek(BinFile,FileHeader.bfOffBits,SEEK_SET);  
135.     if ( (img=new BYTE[size+8])!=NULL)  
136.     {     
137.         if(fread(img+8-int(img)%8,sizeof(BYTE),size,BinFile)!=(unsigned int)size)  
138.         {   
139.             fclose(BinFile);  
140.             delete img;  
141.             img=NULL;  
142.             return NULL;  
143.         }  
144.     }  
145.     fclose(BinFile);  
146.     return img;  
147. }  
148. //  
149. //  write a 24 bit bmp File  
150. //  
151. bool Write24BitImg2BmpFile(BYTE *img,int width,int height,const char * filename)  
152. {     
153.     FILE * BinFile;  
154.     BITMAPFILEHEADER FileHeader;  
155.     BITMAPINFOHEADER BmpHeader;  
156.     bool Suc=true;  
157.     int y,i,extend;  
158.     BYTE *pCur;  
159.   
160.     // Open File  
161.     if((BinFile=fopen(filename,"w+b"))==NULL) {  return false; }  
162.     // Fill the FileHeader  
163.     FileHeader.bfType= ((WORD) ('M' << 8) | 'B');  
164.     FileHeader.bfOffBits=sizeof(BITMAPFILEHEADER)+sizeof(BmpHeader);  
165.     FileHeader.bfSize=FileHeader.bfOffBits+width*height*3L ;  
166.     FileHeader.bfReserved1=0;  
167.     FileHeader.bfReserved2=0;  
168.     if (fwrite((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=false;  
169.     // Fill the ImgHeader  
170.     BmpHeader.biSize = 40;  
171.     BmpHeader.biWidth = width;  
172.     BmpHeader.biHeight = height;  
173.     BmpHeader.biPlanes = 1 ;  
174.     BmpHeader.biBitCount = 24 ;  
175.     BmpHeader.biCompression = 0 ;  
176.     BmpHeader.biSizeImage = 0 ;  
177.     BmpHeader.biXPelsPerMeter = 0;  
178.     BmpHeader.biYPelsPerMeter = 0;  
179.     BmpHeader.biClrUsed = 0;  
180.     BmpHeader.biClrImportant = 0;  
181.     if (fwrite((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=false;  
182.     // write image data  
183.     extend=(width+3)/4*4-width;  
184.     if (extend==0)  
185.     {     
186.         if (fwrite((void *)img,1,width*height*3,BinFile)!=(unsigned int)3*width*height) Suc=false;  
187.     }  
188.     else  
189.     {     
190.         for(y=0,pCur=img;y<height;y++,pCur+=3*width)  
191.         {     
192.             if (fwrite((void *)pCur,1,width*3,BinFile)!=(unsigned int)3*width) Suc=false; // 真实的数据  
193.             for(i=0;i<extend;i++) // 扩充的数据  
194.             {   
195.                 if (fwrite((void *)(pCur+3*(width-1)+0),1,1,BinFile)!=1) Suc=false;  
196.                 if (fwrite((void *)(pCur+3*(width-1)+1),1,1,BinFile)!=1) Suc=false;  
197.                 if (fwrite((void *)(pCur+3*(width-1)+2),1,1,BinFile)!=1) Suc=false;  
198.             }  
199.         }  
200.     }  
201.     // return;  
202.     fclose(BinFile);  
203.     return Suc;  
204. }  
205.   
206. //  
207. // Logarithm Transform  
208. // OrgImg: point to original image  
209. // widht: the width of the image  
210. // height:the height of the image  
211. // LogImg: point to logarithm transform  of the image  
212. //  
213. void LogarithmTransform(int *OrgImg,int width,int height,int *LogImg)  
214. {  
215.     int *pLog=LogImg,*pCur=OrgImg;  
216.     int i,size,temp;  
217.     size=width*height;  
218.     for(i=0;i<size;i++)  
219.     {  
220.         temp=*pCur++;  
221.         if(temp==0) *pLog++=0;  
222.         else *pLog++=log(float(temp))*2048;  
223.     }  
224.     return;  
225. }  
226.   
227. //  计算窗口内像素灰度和  
228.   
229. void Ini(BYTE *pOrgImg,int width,int height,int *sum)  
230. {  
231.     int i,j;  
232.     BYTE *pCur;  
233.     pCur=pOrgImg;  
234.     *sum=*pCur;  
235.     for(i=1;i<height;i++)  
236.         *(sum+i*width)=*(sum+(i-1)*width)+*(pCur+i*width);  
237.     for(j=1;j<width;j++)  
238.         *(sum+j)=*(sum+j-1)+*(pCur+j);  
239.     for(i=1;i<height;i++)  
240.     {  
241.         for(j=1;j<width;j++)   
242.         {  
243.             *(sum+i*width+j)=*(sum+(i-1)*width+j)+*(sum+i*width+j-1)-*(sum+(i-1)*width+j-1)+*(pCur+i*width+j); //卷积计算  
244.         }  
245.     }  
246.     return;  
247. }  
248.   
249. //  
250. // 局部非线性对比度增强  
251. //  
252. void LocalNonlinearStretch(BYTE *OrgImg,int width,int height,int *ResData)  
253. {  
254.     int i,j,k,s,size=width*height;  
255.     int *pData,*sum,sum1;  
256.     double avg,min,max,nor;  
257.     BYTE *pCur=OrgImg;  
258.     sum=new int[width*height];  
259.     pData=ResData+width+1;  
260.     //  Ini(OrgImg,width,height,sum);  
261.     for(i=0;i<height-2;i++,pCur+=2,pData+=2)  
262.     {  
263.         min=*pCur;  
264.         max=*pCur;  
265.         for(k=0;k<3;k++)  
266.         {  
267.             for(s=0;s<3;s++)  
268.             {  
269.                 if(*(pCur+k*width+s)<min) min=*(pCur+k*width+s);  
270.                 else if(*(pCur+k*width+s)>max) max=*(pCur+k*width+s);  
271.             }  
272.         }     
273.         sum1=(*pCur+*(pCur+1)+*(pCur+2)+*(pCur+width)+*(pCur+width+1)+*(pCur+width+2)+*(pCur+width*2)+*(pCur+width*2+1)+*(pCur+width*2+2));  
274.         avg=(sum1-*(pCur+width+1))/8.0;  
275.         nor=(*(pCur+width+1)-min+1)/double(max-min+1);  
276.         *pData=(*(pCur+width+1)+(*(pCur+width+1)-avg)*pow((nor+EPSILON),DELTA)+0.5)*2048;  
277.         pCur++;  
278.         pData++;  
279.         for(j=1;j<width-2;j++,pCur++,pData++)  
280.         {  
281.             min=*pCur;  
282.             max=*pCur;  
283.             for(k=0;k<3;k++)  
284.             {  
285.                 for(s=0;s<3;s++)  
286.                 {  
287.                     if(*(pCur+k*width+s)<min) min=*(pCur+k*width+s);  
288.                     else if(*(pCur+k*width+s)>max) max=*(pCur+k*width+s);  
289.                 }  
290.                 sum1=sum1-*(pCur+k*width-1)+*(pCur+k*width+2);  
291.             }  
292.             //  avg=((*(sum+(i+3)*width+j+3)-*(sum+i*width+j+3)-*(sum+(i+3)*width+j)+*(sum+i*width+j))-*(pCur+width+1))/8.0;  
293.             //  avg=(*pCur+*(pCur+1)+*(pCur+2)+*(pCur+width)+*(pCur+width+2)+*(pCur+width*2)+*(pCur+width*2+1)+*(pCur+width*2+2))/8.0; //  s/8.0*256  
294.             avg=(sum1-*(pCur+width+1))/8.0;  
295.             nor=(*(pCur+width+1)-min+1)/double(max-min+1);  
296.             *pData=(*(pCur+width+1)+(*(pCur+width+1)-avg)*pow((nor+EPSILON),DELTA)+0.5)*2048;  
297.         }  
298.     }  
299.     delete sum;  
300.     return;  
301. }  
302. //  
303. //  Gaussian Template  
304. //  
305. void GaussianTemplate(int *Template,int Tsize,double c)  
306. {  
307.     int *pCur;  
308.     double Lemda,c1=c*c;  
309.     int i,j;  
310.     Lemda=0;  
311.     for(pCur=Template,i=-((Tsize-1)>>1);i<=((Tsize-1)>>1);i++)  
312.     {  
313.         for(j=-((Tsize-1)>>1);j<=((Tsize-1)>>1);j++,pCur++)  
314.         {  
315.             *pCur=(exp(-(i*i+j*j)/c1))*2048;  
316.             Lemda+=*pCur;  
317.         }  
318.     }  
319.     Lemda=2048.0/Lemda;  
320.     for(pCur=Template,i=0;i<Tsize*Tsize;i++,pCur++)  
321.     {  
322.         *pCur=Lemda*(*pCur);  
323.     }  
324.   
325.     return;  
326. }  
327.   
328. //  
329. //  3*3 Gaussian Template  
330. //  
331. void GaussianTemplate2(int *Template,double c)  
332. {  
333.     int *pCur;  
334.     double Lemda,c1=c*c;  
335.     int i,j;  
336.     Lemda=1.0/sqrt(c*c*PI)*0.7*2048;  
337.     for(pCur=Template,i=-1;i<=1;i++)  
338.     {  
339.         for(j=-1;j<=1;j++)  
340.         {  
341.             *pCur++=Lemda*exp(-(i*i+j*j)/c1);  
342.         }  
343.   
344.     }  
345.     return;  
346. }  
347.   
348. //  
349. // 单尺度Retinex  
350. // OrgImg: point to  original image  
351. // widht: the width of the image  
352. // height: the height of the image  
353. // ResImg: point to the result image  
354. //  
355. void SSR(int *LogImg,BYTE *OrgImg,int width,int height,int *ResData,int *Template,int Tsize)  
356. {  
357.     BYTE *pCur=OrgImg;  
358.     int i,j,k,s,size=width*height;  
359.     int temp,*pData,*pCtr,*ptmp,*pRes,temp2;  
360.     double r=1.0/GAMMA;  
361.     memset(ResData,0,sizeof(int)*width*height);  
362.     pRes=ResData+((Tsize-1)/2)*width+((Tsize-1)/2);  
363.     pCtr=LogImg+((Tsize-1)/2)*width+((Tsize-1)/2);  
364.     ptmp=Template;  
365.     for(i=(Tsize-1)/2;i<height-((Tsize-1)/2);i++,pRes+=Tsize-1,pCtr+=Tsize-1,pCur+=Tsize-1)  
366.     {  
367.         for(j=(Tsize-1)/2;j<width-((Tsize-1)/2);j++,pRes++,pCtr++,pCur++)  
368.         {  
369.             temp=0;  
370.             ptmp=Template;  
371.             for(k=0;k<Tsize;k++)  
372.             {  
373.                 for(s=0;s<Tsize;s++)  
374.                 {  
375.                     temp+=(*(pCur+k*width+s)*(*ptmp++));  
376.                 }  
377.             }  
378.             if(temp==0) *pRes=exp(pow(*pCtr>>11,r));  
379.             else   
380.             {  
381.                 temp2=(*pCtr)-(log(float(temp>>22)))*2048;  
382.                 if(temp2>0) *pRes=(exp(pow((temp2>>11),r)))*2048+(temp>>11);  
383.                 else if(temp2<0) *pRes=exp(0-pow(0-(temp2>>11),r))*2048+(temp>>11);  
384.                 else *pRes=(temp>>11);  
385.             }  
386.         }  
387.     }  
388.     //四边不处理  
389.     for(i=0,pRes=ResData,pCur=OrgImg;i<width*(Tsize-1)/2;i++,pCur++,pRes++)  
390.     {  
391.         *pRes=*pCur;  
392.     }  
393.     for(i=(Tsize-1)/2;i<height-(Tsize-1)/2;i++)  
394.     {  
395.         for(j=0;j<(Tsize-1)/2;j++)  
396.         {  
397.             *pRes++=*pCur++;  
398.         }  
399.         pRes+=width-(Tsize-1);  
400.         pCur+=width-(Tsize-1);  
401.         for(j=0;j<(Tsize-1)/2;j++)  
402.         {  
403.             *pRes++=*pCur++;  
404.         }  
405.     }  
406.     for(i=0;i<width*(Tsize-1)/2;i++)  
407.     {  
408.         *pRes++=*pCur++;  
409.     }  
410.     return;  
411. }  
412. /  
413. //  Get Mean And Deviance  
414. /  
415. void GetMeanAndDeviance(int *Temp,int width,int height,int Tsize,int *mean,int *dev)  
416. {  
417.     int i,j,size;  
418.     size=(width-(Tsize-1))*(height-(Tsize-1));  
419.     int *t;  
420.     long double sum;  
421.     for(t=Temp+(Tsize-1)/2*width+(Tsize-1)/2,sum=0,i=(Tsize-1)/2;i<(height-(Tsize-1)/2);i++,t+=Tsize-1)  
422.     {  
423.         for(j=(Tsize-1)/2;j<width-(Tsize-1)/2;j++)  
424.         {  
425.             sum+=*t++;            
426.         }  
427.     }  
428.     *mean=sum/size;  
429.     for(t=Temp+(Tsize-1)/2*width+(Tsize-1)/2,sum=0,i=(Tsize-1)/2;i<height-(Tsize-1)/2;i++,t+=Tsize-1)  
430.     {  
431.         for(j=(Tsize-1)/2;j<width-(Tsize-1)/2;j++)  
432.         {  
433.             sum+=pow(float(*t-*mean),2);  
434.         }  
435.     }  
436.     *dev=sqrt(sum/size);  
437.     return;  
438. }  
439.   
440. /  
441. //  Linear Stretch  
442. // Temp: point to the image before stratching  
443. // widht: the width of the image  
444. // height: the height of the image  
445. // ResImg: point to the resultant image  
446. /  
447. void LinearStretch(int *Temp,int width,int height,int *mean,int *dev,BYTE *ResImg)  
448. {  
449.     BYTE *pRes;  
450.     int *t,min,max,temp,c;  
451.     int i,size=width*height;  
452.     min=*mean-3*(*dev);  
453.     max=*mean+3*(*dev);  
454.     c=255.0/(max-min)*2048;  
455.     for(pRes=ResImg,t=Temp,i=0;i<size;i++,t++,pRes++)  
456.     {  
457.         temp=((*t-min)*c)>>11;  
458.         if(temp>255) *pRes=255;  
459.         else if(temp<0) *pRes=0;  
460.         else *pRes=temp;  
461.     }  
462.   
463.     return;  
464. }  
465.   
466. /* 
467. // 
468. // GAMMA Correction 
469. // 
470. void GammaCorrection(double *OrgData,int width,int heihgt,double *ResData) 
471. { 
472. int i,size; 
473. double *pOrg,*pRes; 
474. for(i=0,pOrg=OrgDat,pRes=ResData;i<size;i++) 
475. { 
476. *pRes++=pow(*pOrg++,1.0/GAMMA); 
477. }  
478. } 
479. */  
480. //  
481. // Contrast  
482. //  
483. void Contrast(BYTE *OrgImg,int x,int y,int width,int height,int blockw,int blockh,double &wcontrast,double &mcontrast)  
484. {  
485.     BYTE *pCur=OrgImg+x*blockw+y*width*blockh;  
486.     double min=10000,max=-1;  
487.     int i,j;  
488.   
489.     for(i=0;i<blockh;i++,pCur=pCur+width-blockw)  
490.     {  
491.         for(j=0;j<blockw;j++,pCur++)  
492.         {  
493.             if(*pCur<min) min=*pCur;  
494.             if(*pCur>max) max=*pCur;  
495.         }  
496.     }  
497.     mcontrast=(max-min+1)/(max+min+1);  
498.     if(min==0) wcontrast=(max+5)/(min+5);  
499.     else wcontrast=max/min;  
500.   
501.     return;  
502. }  
503.   
504. //  
505. //  Measure of Performance  
506. //  
507. void Measure(BYTE *ResImg,int width,int height,double &emee,double &ame)  
508. {  
509.     int k1=8,k2=8,i,j,blockw,blockh;  
510.     double wcontrast,mcontrast;  
511.     blockw=width/k2;  
512.     blockh=height/k1;  
513.     emee=0;  
514.     ame=0;  
515.     for(i=0;i<k1;i++)  
516.     {  
517.         for(j=0;j<k2;j++)  
518.         {  
519.             Contrast(ResImg,i,j,width,height,blockw,blockh,wcontrast,mcontrast);  
520.             emee+=pow(wcontrast,ALPHA_c)*log(wcontrast);  
521.             ame+=pow(mcontrast,ALPHA_c)*log(mcontrast);  
522.         }  
523.     }  
524.     emee=ALPHA_c*emee/(k1*k2);  
525.     ame=ALPHA_c*ame/(k1*k2);  
526.     return ;  
527. }  
528.   
529. /  
530. // Gray Image Process  
531. /  
532. void GrayImageProcess(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
533. {  
534.     int *Data,*LogImg,*Template,mean,dev;  
535.     int Tsize;  
536.     double c;  
537.     Tsize=5;  
538.     c=20;  
539.     Template=new int[Tsize*Tsize];  
540.     Data=new int[width*height];  
541.     LogImg=new int[width*height];  
542.     LocalNonlinearStretch(OrgImg,width,height,Data);      
543.     LogarithmTransform(Data,width,height,LogImg);  
544.     //GaussianTemplate(Template,Tsize,c);  
545.     GaussianTemplate2(Template,0.5);Tsize=3;  
546.     SSR(LogImg,OrgImg,width,height,Data,Template,Tsize);  
547.     GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
548.     LinearStretch(Data,width,height,&mean,&dev,ResImg);  
549.     delete Template;  
550.     delete Data;  
551.     delete LogImg;  
552.     return;  
553. }  
554.   
555. /  
556. // Color Image Process  
557. /  
558. void ColorImageProcess(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
559. {  
560.     BYTE *Value;  
561.     int i,j,Tsize,temp;  
562.     int *Data,*LogImg,*Template,*Percent,mean,dev;  
563.     double c,emee,ame;  
564.     Tsize=5;  
565.     c=0.75;  
566.     Template=new int[Tsize*Tsize];  
567.     Data=new int[width*height];  
568.     LogImg=new int[width*height];  
569.     Percent=new int[width*height*3];  
570.     Value=new BYTE[width*height];  
571.     memset(Value,0,sizeof(BYTE)*width*height);  
572.     for(j=0;j<width*height;j++)  
573.     {    
574.         temp=0;  
575.         for(i=0;i<3;i++)  
576.         {   
577.             temp+=*(OrgImg+j*3+i);  
578.         }  
579.         for(i=0;i<3;i++)  
580.         {   
581.             *(Percent+j*3+i)=*(OrgImg+j*3+i)/double(temp)*2048;  
582.             *(Value+j)+=(*(OrgImg+j*3+i)*(*(Percent+j*3+i)))>>11;  
583.         }  
584.     }  
585.     LocalNonlinearStretch(Value,width,height,Data);  
586.     LogarithmTransform(Data,width,height,LogImg);  
587.     //  GaussianTemplate(Template,Tsize,c);  
588.     GaussianTemplate2(Template,0.5);Tsize=3;  
589.     SSR(LogImg,Value,width,height,Data,Template,Tsize);  
590.     GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
591.     LinearStretch(Data,width,height,&mean,&dev,Value);  
592.     for(j=0;j<width*height;j++)  
593.     {  
594.         for(i=0;i<3;i++)  
595.         {   
596.             temp=(*(Percent+j*3+i)*(*(Value+j))*3)>>11;  
597.             if(temp>255) *(ResImg+j*3+i)=255;  
598.             else *(ResImg+j*3+i)=temp;  
599.         }  
600.     }  
601.     printf("*****彩色图像三通道联合处理结果******************\n");  
602.     Measure(Value,width,height,emee,ame);  
603.     cout<<"EMEE:"<<emee<<endl;  
604.     delete Template;  
605.     delete Data;  
606.     delete LogImg;  
607.     delete Value;  
608.     delete Percent;  
609.     return;  
610. }  
611.   
612.   
613. //  
614. //  Color Image process 2  
615. //  三通道分别处理  
616. //  
617. void ColorImageProcess2(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
618. {  
619.     BYTE *Value;  
620.     int i,j,Tsize;  
621.     int *Data,*LogImg,*Template,mean,dev;  
622.     double c,emee=0,ame=0,x1,x2;  
623.     Tsize=5;  
624.     c=0.5;  
625.     Template=new int[Tsize*Tsize];  
626.     Value=new BYTE[width*height];     
627.     Data=new int[width*height];  
628.     LogImg=new int[width*height];  
629.     GaussianTemplate(Template,Tsize,c);  
630.     //      GaussianTemplate2(Template,0.5);    Tsize=3;  
631.     for(i=0;i<3;i++)  
632.     {  
633.         for(j=0;j<width*height;j++)  
634.         {  
635.             *(Value+j)=*(OrgImg+j*3+i);  
636.         }  
637.         LocalNonlinearStretch(Value,width,height,Data);  
638.         LogarithmTransform(Data,width,height,LogImg);  
639.         SSR(LogImg,Value,width,height,Data,Template,Tsize);   
640.         GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
641.         LinearStretch(Data,width,height,&mean,&dev,Value);  
642.         Measure(Value,width,height,x1,x2);  
643.         emee+=x1;  
644.         ame+=x2;  
645.         for(j=0;j<width*height;j++)  
646.         {  
647.             *(ResImg+j*3+i)=*(Value+j);  
648.         }  
649.     }  
650.     printf("*****彩色图像三通道分别处理结果******************\n");  
651.     cout<<"EMEE:"<<emee/3.0<<endl;          
652.     delete Template;  
653.     delete Value;  
654.     delete Data;  
655.     delete LogImg;  
656.     return;  
657. }  
658. //  
659. //  main  
660. //  
661. int main()  
662. {     
663.     BYTE *OrgImg,*ResImg,*p1,*p2;  
664.     int width,height,i,j,suc,area1,area2,Tsize;  
665.     bool isRGB,ret;  
666.     clock_t t1,t2;  
667.     double emee,ame;  
668.     char ch;  
669.     int *offsetdata=new int[2];  
670.     for(i=0;i<2;i++)  
671.     {  
672.         *(offsetdata+i)=0X80808080;  
673.     }      
674.     system( "cls" );      
675.     printf("******中心/环绕Retienx算法******************\n");  
676.     printf("       1.处理灰度图像\n");  
677.     printf("       2.处理彩色图像\n");  
678.     printf("*********************************************\n");  
679.     printf("请选择(1或2): ");     
680.     do  
681.     {  
682.         cin>>ch;   
683.     }while( ch != '1' && ch != '2');  
684.   
685.     //system ( "cls" );  
686.   
687.   
688.     if ( ch == '1')  
689.         isRGB=false;  
690.     else if ( ch == '2')  
691.         isRGB=true;  
692.     // open file  
693.   
694.     string image_name;  
695.     cout<<endl<<"输入图像名：";  
696.     cin>>image_name;  
697.   
698.     //Mat image_dst;  
699.   
700.     //if (!isRGB)  
701.     //{  
702.     //  image_dst = imread(image_name,0);  
703.     //}  
704.     //else  
705.     //{  
706.     //  image_dst = imread(image_name,1);  
707.     //}  
708.   
709.     //  
710.     //  
711.     //if(image_dst.empty())  
712.     //{  
713.     //  return -1;  
714.     //} //是否加载成功  
715.   
716.     //imshow(image_name,image_dst);  
717.     //  
718.   
719.     // width = image_dst.cols;    
720.     // height = image_dst.rows;    
721.     //int channel = image_dst.channels();   
722.     // int step = width * channel* 1;  
723.     //uchar* ps = NULL;  
724.     //  
725.     //p1 = new BYTE[width*height*channel+8];  
726.   
727.     //for (int i = 0; i < height; i++)    
728.     //{    
729.     //  ps = image_dst.ptr<uchar>(i);    
730.     //  for (int j = 0; j < width; j++)    
731.     //  {    
732.     //      if (1 == channel)    
733.     //      {    
734.     //          *(p1 + i*step + j) = ps[j];    
735.     //            
736.     //      }    
737.     //      else if (3 == channel)    
738.     //      {    
739.     //          *(p1 + i*step + j*3) = ps[j*3 + 2];    
740.     //          *(p1 + i*step + j*3 + 1) = ps[j*3 + 1];    
741.     //          *(p1 + i*step + j*3 + 2) = ps[j*3];    
742.     //      }    
743.     //  }    
744.     //}    
745.   
746.   
747.     if(!isRGB)  
748.     p1=Read8BitBmpFile2Img(image_name.c_str(),&width,&height);  
749.     else  
750.     p1=Read24BitBmpFile2Img(image_name.c_str(),&width,&height);  
751.       
752.   
753.     if (p1==NULL)  
754.     {    
755.         printf("*fopen err!\n");  
756.         delete p1;  
757.         return 0;  
758.     }  
759.     if(width%64!=0||height%64!=0)  
760.     {  
761.         cout<<"图像大小需是64的倍数"<<endl;  
762.         delete p1;  
763.         return 0;  
764.     }  
765.   
766.     area1=width*height;  
767.   
768.     if(!isRGB)  
769.     {  
770.         OrgImg=p1+8-int(p1)%8;    
771.         p2=new BYTE[width*height+8];  
772.         ResImg=p2+8-int(p2)%8;  
773.         t2=clock();  
774.         GrayImageProcess(OrgImg,width,height,ResImg);  
775.         t1=clock();  
776.         printf("*****灰度图像处理结果******************\n");  
777.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
778.         Measure(ResImg,width,height,emee,ame);  
779.         cout<<"EMEE:"<<emee<<endl;  
780.         suc=Write8BitImg2BmpFile(ResImg,width,height,"result_1.bmp");  
781.     }  
782.     else  
783.     {     
784.         OrgImg=p1+8-int(p1)%8;    
785.         p2=new BYTE[width*height*3+8];        
786.         ResImg=p2+8-int(p2)%8;  
787.         t2=clock();  
788.         ColorImageProcess(OrgImg,width,height,ResImg);  
789.         t1=clock();  
790.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
791.         t2=clock();  
792.         suc=Write24BitImg2BmpFile(ResImg,width,height,"result_1.bmp");  
793.         ColorImageProcess2(OrgImg,width,height,ResImg);  
794.         t1=clock();  
795.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
796.         suc=Write24BitImg2BmpFile(ResImg,width,height,"result_2.bmp");  
797.     }  
798.     if(suc==-1)  
799.     {  
800.         printf("*fwrite err!\n");    
801.     }  
802.   
803.     Mat result1 = imread("result_1.bmp",1);  
804.     Mat result2 = imread("result_2.bmp",1);  
805.   
806.     imshow("result1",result1);  
807.     imshow("result2",result2);  
808.     //release mem  
809.     delete p1;  
810.     delete p2;  
811.       
812.   
813.     waitKey(0);  
814.     return 0;  
815. }  

 

 

后面用opencv改写了一下主要想把像素数据写到BYTE *指向的内存空间中，这样的话可以加载其他格式的图像文件了:但是出现了一些问题，可能跟之前作者代码里面的限制有关，必须是64的倍数。还写了一个功能实现把像素数据写到txt中保存下来，各位凑活看。

参见帖子：

http://bbs.csdn.net/topics/391005171

http://bbs.csdn.net/topics/391004682

 

代码如下：

 

[cpp]  view plain  copy

 print ?

1. #include <stdio.h>  
2. #include <string.h>  
3. #include <windows.h>  
4. #include <cmath>  
5. #include <time.h>  
6. #include <iostream>  
7.   
8. #include "opencv2/core/core.hpp"  
9. #include "opencv2/highgui/highgui.hpp"  
10. #include "opencv2/imgproc/imgproc.hpp"  
11.   
12. #pragma comment(lib,"opencv_core2410d.lib")                    
13. #pragma comment(lib,"opencv_highgui2410d.lib")                    
14. #pragma comment(lib,"opencv_imgproc2410d.lib")    
15.   
16.   
17. using namespace std;   
18. using namespace cv;  
19.   
20. #define EPSILON 1  
21. #define DELTA 1  
22. #define GAMMA 0.9     
23. #define PI 3.1415926  
24. #define ALPHA_c 1  
25.   
26. int size_mem = 0;  
27.   
28. //  
29. //  Read a 8 bit bmp File  
30. //  
31. BYTE *Read8BitBmpFile2Img(const char * filename,int *width,int *height)  
32. {   FILE * BinFile;  
33. BITMAPFILEHEADER FileHeader;  
34. BITMAPINFOHEADER BmpHeader;  
35. BYTE *img;  
36. int size;  
37. int Suc=1;  
38. // Open File  
39. *width=*height=0;  
40. if((BinFile=fopen(filename,"rb"))==NULL) return NULL;  
41. // Read Struct Info  
42. if (fread((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
43. if (fread((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
44. if (Suc==-1) { fclose(BinFile); return NULL; }  
45. // Read Image Data  
46. *width=(BmpHeader.biWidth+3)/4*4;  
47. *height=BmpHeader.biHeight;  
48. size=(BmpHeader.biWidth+3)/4*4*BmpHeader.biHeight;  
49. fseek(BinFile,FileHeader.bfOffBits,SEEK_SET);  
50. if ( (img=new BYTE[size+8])!=NULL)  
51. {   if(fread(img+8-int(img)%8,sizeof(BYTE),size,BinFile)!=(unsigned int)size)  
52. { fclose(BinFile);  
53. delete img;  
54. img=NULL;  
55. return NULL;  
56. }  
57. }  
58. fclose(BinFile);  
59. return img;  
60. }  
61.   
62. /  
63. // Write a 8 bit bmp File  
64. /  
65. int Write8BitImg2BmpFile(BYTE *img,int width,int height,const char * filename)  
66. {   FILE * BinFile;  
67. BITMAPFILEHEADER FileHeader;  
68. BITMAPINFOHEADER BmpHeader;  
69. BYTE p[4];  
70. int i,Suc=1;  
71.   
72. // Open File  
73. if((BinFile=fopen(filename,"w+b"))==NULL) {  return -1; }  
74. //  Fill the FileHeade)  
75. FileHeader.bfType= ((WORD) ('M' << 8) | 'B');  
76. FileHeader.bfOffBits=sizeof(BITMAPFILEHEADER)+sizeof(BmpHeader)+256*4L;  
77. FileHeader.bfSize=FileHeader.bfOffBits+width*height ;  
78. FileHeader.bfReserved1=0;  
79. FileHeader.bfReserved2=0;  
80. if (fwrite((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
81. // Fill the ImgHeader  
82. BmpHeader.biSize = 40;  
83. BmpHeader.biWidth = width;  
84. BmpHeader.biHeight = height;  
85. BmpHeader.biPlanes = 1 ;  
86. BmpHeader.biBitCount = 8 ;  
87. BmpHeader.biCompression = 0 ;  
88. BmpHeader.biSizeImage = 0 ;  
89. BmpHeader.biXPelsPerMeter = 0;  
90. BmpHeader.biYPelsPerMeter = 0;  
91. BmpHeader.biClrUsed = 0;  
92. BmpHeader.biClrImportant = 0;  
93. if (fwrite((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
94. // write Pallete  
95. for (i=0,p[3]=0;i<256;i++)   
96. {    
97.     p[3]=0;  
98.     p[0]=p[1]=p[2]=i; // blue,green,red;  
99.     if (fwrite((void *)p,1,4,BinFile)!=4) { Suc=-1; break; }  
100. }  
101. // write image data  
102. if (fwrite((void *)img,1,width*height,BinFile)!=(unsigned int) width*height) Suc=-1;  
103. // return;  
104. fclose(BinFile);  
105. return Suc;  
106. }  
107.   
108. //  
109. //  Read a 24 bit bmp File  
110. //  
111. BYTE *Read24BitBmpFile2Img(const char * filename,int *width,int *height)  
112. {    
113.     FILE * BinFile;  
114.     BITMAPFILEHEADER FileHeader;  
115.     BITMAPINFOHEADER BmpHeader;  
116.     BYTE *img;  
117.     int size;  
118.     int Suc=1;  
119.   
120.     // Open File  
121.     *width=*height=0;  
122.     if((BinFile=fopen(filename,"rb"))==NULL) return NULL;  
123.     // Read Struct Info  
124.     if (fread((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=-1;  
125.     if (fread((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=-1;  
126.     if (Suc==-1) { fclose(BinFile); return NULL; }  
127.     // Read Image Data  
128.     *width=(BmpHeader.biWidth+3)/4*4;  
129.     *height=BmpHeader.biHeight;  
130.     size=(*width)*(*height)*3;  
131.     fseek(BinFile,FileHeader.bfOffBits,SEEK_SET);  
132.     if ( (img=new BYTE[size+8])!=NULL)  
133.     {     
134.         //size_mem = *((int *)img - 1);  
135.         if(fread(img+8-int(img)%8,sizeof(BYTE),size,BinFile)!=(unsigned int)size)  
136.         {   
137.             fclose(BinFile);  
138.             delete img;  
139.             img=NULL;  
140.             return NULL;  
141.         }  
142.     }  
143.     fclose(BinFile);  
144.     return img;  
145. }  
146. //  
147. //  write a 24 bit bmp File  
148. //  
149. bool Write24BitImg2BmpFile(BYTE *img,int width,int height,const char * filename)  
150. {     
151.     FILE * BinFile;  
152.     BITMAPFILEHEADER FileHeader;  
153.     BITMAPINFOHEADER BmpHeader;  
154.     bool Suc=true;  
155.     int y,i,extend;  
156.     BYTE *pCur;  
157.   
158.     // Open File  
159.     if((BinFile=fopen(filename,"w+b"))==NULL) {  return false; }  
160.     // Fill the FileHeader  
161.     FileHeader.bfType= ((WORD) ('M' << 8) | 'B');  
162.     FileHeader.bfOffBits=sizeof(BITMAPFILEHEADER)+sizeof(BmpHeader);  
163.     FileHeader.bfSize=FileHeader.bfOffBits+width*height*3L ;  
164.     FileHeader.bfReserved1=0;  
165.     FileHeader.bfReserved2=0;  
166.     if (fwrite((void *)&FileHeader,1,sizeof(FileHeader),BinFile)!=sizeof(FileHeader)) Suc=false;  
167.     // Fill the ImgHeader  
168.     BmpHeader.biSize = 40;  
169.     BmpHeader.biWidth = width;  
170.     BmpHeader.biHeight = height;  
171.     BmpHeader.biPlanes = 1 ;  
172.     BmpHeader.biBitCount = 24 ;  
173.     BmpHeader.biCompression = 0 ;  
174.     BmpHeader.biSizeImage = 0 ;  
175.     BmpHeader.biXPelsPerMeter = 0;  
176.     BmpHeader.biYPelsPerMeter = 0;  
177.     BmpHeader.biClrUsed = 0;  
178.     BmpHeader.biClrImportant = 0;  
179.     if (fwrite((void *)&BmpHeader,1,sizeof(BmpHeader),BinFile)!=sizeof(BmpHeader)) Suc=false;  
180.     // write image data  
181.     extend=(width+3)/4*4-width;  
182.     if (extend==0)  
183.     {     
184.         if (fwrite((void *)img,1,width*height*3,BinFile)!=(unsigned int)3*width*height) Suc=false;  
185.     }  
186.     else  
187.     {     
188.         for(y=0,pCur=img;y<height;y++,pCur+=3*width)  
189.         {     
190.             if (fwrite((void *)pCur,1,width*3,BinFile)!=(unsigned int)3*width) Suc=false; // 真实的数据  
191.             for(i=0;i<extend;i++) // 扩充的数据  
192.             {   
193.                 if (fwrite((void *)(pCur+3*(width-1)+0),1,1,BinFile)!=1) Suc=false;  
194.                 if (fwrite((void *)(pCur+3*(width-1)+1),1,1,BinFile)!=1) Suc=false;  
195.                 if (fwrite((void *)(pCur+3*(width-1)+2),1,1,BinFile)!=1) Suc=false;  
196.             }  
197.         }  
198.     }  
199.     // return;  
200.     fclose(BinFile);  
201.     return Suc;  
202. }  
203.   
204. //  
205. // Logarithm Transform  
206. // OrgImg: point to original image  
207. // widht: the width of the image  
208. // height:the height of the image  
209. // LogImg: point to logarithm transform  of the image  
210. //  
211. void LogarithmTransform(int *OrgImg,int width,int height,int *LogImg)  
212. {  
213.     int *pLog=LogImg,*pCur=OrgImg;  
214.     int i,size,temp;  
215.     size=width*height;  
216.     for(i=0;i<size;i++)  
217.     {  
218.         temp=*pCur++;  
219.         if(temp==0) *pLog++=0;  
220.         else *pLog++=log(float(temp))*2048;  
221.     }  
222.     return;  
223. }  
224.   
225. //  计算窗口内像素灰度和  
226.   
227. void Ini(BYTE *pOrgImg,int width,int height,int *sum)  
228. {  
229.     int i,j;  
230.     BYTE *pCur;  
231.     pCur=pOrgImg;  
232.     *sum=*pCur;  
233.     for(i=1;i<height;i++)  
234.         *(sum+i*width)=*(sum+(i-1)*width)+*(pCur+i*width);  
235.     for(j=1;j<width;j++)  
236.         *(sum+j)=*(sum+j-1)+*(pCur+j);  
237.     for(i=1;i<height;i++)  
238.     {  
239.         for(j=1;j<width;j++)   
240.         {  
241.             *(sum+i*width+j)=*(sum+(i-1)*width+j)+*(sum+i*width+j-1)-*(sum+(i-1)*width+j-1)+*(pCur+i*width+j); //卷积计算  
242.         }  
243.     }  
244.     return;  
245. }  
246.   
247. //  
248. // 局部非线性对比度增强  
249. //  
250. void LocalNonlinearStretch(BYTE *OrgImg,int width,int height,int *ResData)  
251. {  
252.     int i,j,k,s,size=width*height;  
253.     int *pData,*sum,sum1;  
254.     double avg,min,max,nor;  
255.     BYTE *pCur=OrgImg;  
256.     sum=new int[width*height];  
257.     pData=ResData+width+1;  
258.     //  Ini(OrgImg,width,height,sum);  
259.     for(i=0;i<height-2;i++,pCur+=2,pData+=2)  
260.     {  
261.         min=*pCur;  
262.         max=*pCur;  
263.         for(k=0;k<3;k++)  
264.         {  
265.             for(s=0;s<3;s++)  
266.             {  
267.                 if(*(pCur+k*width+s)<min) min=*(pCur+k*width+s);  
268.                 else if(*(pCur+k*width+s)>max) max=*(pCur+k*width+s);  
269.             }  
270.         }     
271.         sum1=(*pCur+*(pCur+1)+*(pCur+2)+*(pCur+width)+*(pCur+width+1)+*(pCur+width+2)+*(pCur+width*2)+*(pCur+width*2+1)+*(pCur+width*2+2));  
272.         avg=(sum1-*(pCur+width+1))/8.0;  
273.         nor=(*(pCur+width+1)-min+1)/double(max-min+1);  
274.         *pData=(*(pCur+width+1)+(*(pCur+width+1)-avg)*pow((nor+EPSILON),DELTA)+0.5)*2048;  
275.         pCur++;  
276.         pData++;  
277.         for(j=1;j<width-2;j++,pCur++,pData++)  
278.         {  
279.             min=*pCur;  
280.             max=*pCur;  
281.             for(k=0;k<3;k++)  
282.             {  
283.                 for(s=0;s<3;s++)  
284.                 {  
285.                     if(*(pCur+k*width+s)<min) min=*(pCur+k*width+s);  
286.                     else if(*(pCur+k*width+s)>max) max=*(pCur+k*width+s);  
287.                 }  
288.                 sum1=sum1-*(pCur+k*width-1)+*(pCur+k*width+2);  
289.             }  
290.             //  avg=((*(sum+(i+3)*width+j+3)-*(sum+i*width+j+3)-*(sum+(i+3)*width+j)+*(sum+i*width+j))-*(pCur+width+1))/8.0;  
291.             //  avg=(*pCur+*(pCur+1)+*(pCur+2)+*(pCur+width)+*(pCur+width+2)+*(pCur+width*2)+*(pCur+width*2+1)+*(pCur+width*2+2))/8.0; //  s/8.0*256  
292.             avg=(sum1-*(pCur+width+1))/8.0;  
293.             nor=(*(pCur+width+1)-min+1)/double(max-min+1);  
294.             *pData=(*(pCur+width+1)+(*(pCur+width+1)-avg)*pow((nor+EPSILON),DELTA)+0.5)*2048;  
295.         }  
296.     }  
297.     delete sum;  
298.     return;  
299. }  
300. //  
301. //  Gaussian Template  
302. //  
303. void GaussianTemplate(int *Template,int Tsize,double c)  
304. {  
305.     int *pCur;  
306.     double Lemda,c1=c*c;  
307.     int i,j;  
308.     Lemda=0;  
309.     for(pCur=Template,i=-((Tsize-1)>>1);i<=((Tsize-1)>>1);i++)  
310.     {  
311.         for(j=-((Tsize-1)>>1);j<=((Tsize-1)>>1);j++,pCur++)  
312.         {  
313.             *pCur=(exp(-(i*i+j*j)/c1))*2048;  
314.             Lemda+=*pCur;  
315.         }  
316.     }  
317.     Lemda=2048.0/Lemda;  
318.     for(pCur=Template,i=0;i<Tsize*Tsize;i++,pCur++)  
319.     {  
320.         *pCur=Lemda*(*pCur);  
321.     }  
322.   
323.     return;  
324. }  
325.   
326. //  
327. //  3*3 Gaussian Template  
328. //  
329. void GaussianTemplate2(int *Template,double c)  
330. {  
331.     int *pCur;  
332.     double Lemda,c1=c*c;  
333.     int i,j;  
334.     Lemda=1.0/sqrt(c*c*PI)*0.7*2048;  
335.     for(pCur=Template,i=-1;i<=1;i++)  
336.     {  
337.         for(j=-1;j<=1;j++)  
338.         {  
339.             *pCur++=Lemda*exp(-(i*i+j*j)/c1);  
340.         }  
341.   
342.     }  
343.     return;  
344. }  
345.   
346. //  
347. // 单尺度Retinex  
348. // OrgImg: point to  original image  
349. // widht: the width of the image  
350. // height: the height of the image  
351. // ResImg: point to the result image  
352. //  
353. void SSR(int *LogImg,BYTE *OrgImg,int width,int height,int *ResData,int *Template,int Tsize)  
354. {  
355.     BYTE *pCur=OrgImg;  
356.     int i,j,k,s,size=width*height;  
357.     int temp,*pData,*pCtr,*ptmp,*pRes,temp2;  
358.     double r=1.0/GAMMA;  
359.     memset(ResData,0,sizeof(int)*width*height);  
360.     pRes=ResData+((Tsize-1)/2)*width+((Tsize-1)/2);  
361.     pCtr=LogImg+((Tsize-1)/2)*width+((Tsize-1)/2);  
362.     ptmp=Template;  
363.     for(i=(Tsize-1)/2;i<height-((Tsize-1)/2);i++,pRes+=Tsize-1,pCtr+=Tsize-1,pCur+=Tsize-1)  
364.     {  
365.         for(j=(Tsize-1)/2;j<width-((Tsize-1)/2);j++,pRes++,pCtr++,pCur++)  
366.         {  
367.             temp=0;  
368.             ptmp=Template;  
369.             for(k=0;k<Tsize;k++)  
370.             {  
371.                 for(s=0;s<Tsize;s++)  
372.                 {  
373.                     temp+=(*(pCur+k*width+s)*(*ptmp++));  
374.                 }  
375.             }  
376.             if(temp==0) *pRes=exp(pow(*pCtr>>11,r));  
377.             else   
378.             {  
379.                 temp2=(*pCtr)-(log(float(temp>>22)))*2048;  
380.                 if(temp2>0) *pRes=(exp(pow((temp2>>11),r)))*2048+(temp>>11);  
381.                 else if(temp2<0) *pRes=exp(0-pow(0-(temp2>>11),r))*2048+(temp>>11);  
382.                 else *pRes=(temp>>11);  
383.             }  
384.         }  
385.     }  
386.     //四边不处理  
387.     for(i=0,pRes=ResData,pCur=OrgImg;i<width*(Tsize-1)/2;i++,pCur++,pRes++)  
388.     {  
389.         *pRes=*pCur;  
390.     }  
391.     for(i=(Tsize-1)/2;i<height-(Tsize-1)/2;i++)  
392.     {  
393.         for(j=0;j<(Tsize-1)/2;j++)  
394.         {  
395.             *pRes++=*pCur++;  
396.         }  
397.         pRes+=width-(Tsize-1);  
398.         pCur+=width-(Tsize-1);  
399.         for(j=0;j<(Tsize-1)/2;j++)  
400.         {  
401.             *pRes++=*pCur++;  
402.         }  
403.     }  
404.     for(i=0;i<width*(Tsize-1)/2;i++)  
405.     {  
406.         *pRes++=*pCur++;  
407.     }  
408.     return;  
409. }  
410. /  
411. //  Get Mean And Deviance  
412. /  
413. void GetMeanAndDeviance(int *Temp,int width,int height,int Tsize,int *mean,int *dev)  
414. {  
415.     int i,j,size;  
416.     size=(width-(Tsize-1))*(height-(Tsize-1));  
417.     int *t;  
418.     long double sum;  
419.     for(t=Temp+(Tsize-1)/2*width+(Tsize-1)/2,sum=0,i=(Tsize-1)/2;i<(height-(Tsize-1)/2);i++,t+=Tsize-1)  
420.     {  
421.         for(j=(Tsize-1)/2;j<width-(Tsize-1)/2;j++)  
422.         {  
423.             sum+=*t++;            
424.         }  
425.     }  
426.     *mean=sum/size;  
427.     for(t=Temp+(Tsize-1)/2*width+(Tsize-1)/2,sum=0,i=(Tsize-1)/2;i<height-(Tsize-1)/2;i++,t+=Tsize-1)  
428.     {  
429.         for(j=(Tsize-1)/2;j<width-(Tsize-1)/2;j++)  
430.         {  
431.             sum+=pow(float(*t-*mean),2);  
432.         }  
433.     }  
434.     *dev=sqrt(sum/size);  
435.     return;  
436. }  
437.   
438. /  
439. //  Linear Stretch  
440. // Temp: point to the image before stratching  
441. // widht: the width of the image  
442. // height: the height of the image  
443. // ResImg: point to the resultant image  
444. /  
445. void LinearStretch(int *Temp,int width,int height,int *mean,int *dev,BYTE *ResImg)  
446. {  
447.     BYTE *pRes;  
448.     int *t,min,max,temp,c;  
449.     int i,size=width*height;  
450.     min=*mean-3*(*dev);  
451.     max=*mean+3*(*dev);  
452.     c=255.0/(max-min)*2048;  
453.     for(pRes=ResImg,t=Temp,i=0;i<size;i++,t++,pRes++)  
454.     {  
455.         temp=((*t-min)*c)>>11;  
456.         if(temp>255) *pRes=255;  
457.         else if(temp<0) *pRes=0;  
458.         else *pRes=temp;  
459.     }  
460.   
461.     return;  
462. }  
463.   
464. /* 
465. // 
466. // GAMMA Correction 
467. // 
468. void GammaCorrection(double *OrgData,int width,int heihgt,double *ResData) 
469. { 
470. int i,size; 
471. double *pOrg,*pRes; 
472. for(i=0,pOrg=OrgDat,pRes=ResData;i<size;i++) 
473. { 
474. *pRes++=pow(*pOrg++,1.0/GAMMA); 
475. }  
476. } 
477. */  
478. //  
479. // Contrast  
480. //  
481. void Contrast(BYTE *OrgImg,int x,int y,int width,int height,int blockw,int blockh,double &wcontrast,double &mcontrast)  
482. {  
483.     BYTE *pCur=OrgImg+x*blockw+y*width*blockh;  
484.     double min=10000,max=-1;  
485.     int i,j;  
486.   
487.     for(i=0;i<blockh;i++,pCur=pCur+width-blockw)  
488.     {  
489.         for(j=0;j<blockw;j++,pCur++)  
490.         {  
491.             if(*pCur<min) min=*pCur;  
492.             if(*pCur>max) max=*pCur;  
493.         }  
494.     }  
495.     mcontrast=(max-min+1)/(max+min+1);  
496.     if(min==0) wcontrast=(max+5)/(min+5);  
497.     else wcontrast=max/min;  
498.   
499.     return;  
500. }  
501.   
502. //  
503. //  Measure of Performance  
504. //  
505. void Measure(BYTE *ResImg,int width,int height,double &emee,double &ame)  
506. {  
507.     int k1=8,k2=8,i,j,blockw,blockh;  
508.     double wcontrast,mcontrast;  
509.     blockw=width/k2;  
510.     blockh=height/k1;  
511.     emee=0;  
512.     ame=0;  
513.     for(i=0;i<k1;i++)  
514.     {  
515.         for(j=0;j<k2;j++)  
516.         {  
517.             Contrast(ResImg,i,j,width,height,blockw,blockh,wcontrast,mcontrast);  
518.             emee+=pow(wcontrast,ALPHA_c)*log(wcontrast);  
519.             ame+=pow(mcontrast,ALPHA_c)*log(mcontrast);  
520.         }  
521.     }  
522.     emee=ALPHA_c*emee/(k1*k2);  
523.     ame=ALPHA_c*ame/(k1*k2);  
524.     return ;  
525. }  
526.   
527. /  
528. // Gray Image Process  
529. /  
530. void GrayImageProcess(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
531. {  
532.     int *Data,*LogImg,*Template,mean,dev;  
533.     int Tsize;  
534.     double c;  
535.     Tsize=5;  
536.     c=20;  
537.     Template=new int[Tsize*Tsize];  
538.     Data=new int[width*height];  
539.     LogImg=new int[width*height];  
540.     LocalNonlinearStretch(OrgImg,width,height,Data);      
541.     LogarithmTransform(Data,width,height,LogImg);  
542.     //GaussianTemplate(Template,Tsize,c);  
543.     GaussianTemplate2(Template,0.5);Tsize=3;  
544.     SSR(LogImg,OrgImg,width,height,Data,Template,Tsize);  
545.     GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
546.     LinearStretch(Data,width,height,&mean,&dev,ResImg);  
547.     delete Template;  
548.     delete Data;  
549.     delete LogImg;  
550.     return;  
551. }  
552.   
553. /  
554. // Color Image Process  
555. /  
556. void ColorImageProcess(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
557. {  
558.     BYTE *Value;  
559.     int i,j,Tsize,temp;  
560.     int *Data,*LogImg,*Template,*Percent,mean,dev;  
561.     double c,emee,ame;  
562.     Tsize=5;  
563.     c=0.75;  
564.     Template=new int[Tsize*Tsize];  
565.     Data=new int[width*height];  
566.     LogImg=new int[width*height];  
567.     Percent=new int[width*height*3];  
568.     Value=new BYTE[width*height];  
569.     memset(Value,0,sizeof(BYTE)*width*height);  
570.     for(j=0;j<width*height;j++)  
571.     {    
572.         temp=0;  
573.         for(i=0;i<3;i++)  
574.         {   
575.             temp+=*(OrgImg+j*3+i);  
576.         }  
577.         for(i=0;i<3;i++)  
578.         {   
579.             *(Percent+j*3+i)=*(OrgImg+j*3+i)/double(temp)*2048;  
580.             *(Value+j)+=(*(OrgImg+j*3+i)*(*(Percent+j*3+i)))>>11;  
581.         }  
582.     }  
583.     LocalNonlinearStretch(Value,width,height,Data);  
584.     LogarithmTransform(Data,width,height,LogImg);  
585.     //  GaussianTemplate(Template,Tsize,c);  
586.     GaussianTemplate2(Template,0.5);Tsize=3;  
587.     SSR(LogImg,Value,width,height,Data,Template,Tsize);  
588.     GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
589.     LinearStretch(Data,width,height,&mean,&dev,Value);  
590.     for(j=0;j<width*height;j++)  
591.     {  
592.         for(i=0;i<3;i++)  
593.         {   
594.             temp=(*(Percent+j*3+i)*(*(Value+j))*3)>>11;  
595.             if(temp>255) *(ResImg+j*3+i)=255;  
596.             else *(ResImg+j*3+i)=temp;  
597.         }  
598.     }  
599.     printf("*****彩色图像三通道联合处理结果******************\n");  
600.     Measure(Value,width,height,emee,ame);  
601.     cout<<"EMEE:"<<emee<<endl;  
602.     delete Template;  
603.     delete Data;  
604.     delete LogImg;  
605.     delete Value;  
606.     delete Percent;  
607.     return;  
608. }  
609.   
610.   
611. //  
612. //  Color Image process 2  
613. //  三通道分别处理  
614. //  
615. void ColorImageProcess2(BYTE *OrgImg,int width,int height,BYTE *ResImg)  
616. {  
617.     BYTE *Value;  
618.     int i,j,Tsize;  
619.     int *Data,*LogImg,*Template,mean,dev;  
620.     double c,emee=0,ame=0,x1,x2;  
621.     Tsize=5;  
622.     c=0.5;  
623.     Template=new int[Tsize*Tsize];  
624.     Value=new BYTE[width*height];     
625.     Data=new int[width*height];  
626.     LogImg=new int[width*height];  
627.     GaussianTemplate(Template,Tsize,c);  
628.     //      GaussianTemplate2(Template,0.5);    Tsize=3;  
629.     for(i=0;i<3;i++)  
630.     {  
631.         for(j=0;j<width*height;j++)  
632.         {  
633.             *(Value+j)=*(OrgImg+j*3+i);  
634.         }  
635.         LocalNonlinearStretch(Value,width,height,Data);  
636.         LogarithmTransform(Data,width,height,LogImg);  
637.         SSR(LogImg,Value,width,height,Data,Template,Tsize);   
638.         GetMeanAndDeviance(Data,width,height,Tsize,&mean,&dev);  
639.         LinearStretch(Data,width,height,&mean,&dev,Value);  
640.         Measure(Value,width,height,x1,x2);  
641.         emee+=x1;  
642.         ame+=x2;  
643.         for(j=0;j<width*height;j++)  
644.         {  
645.             *(ResImg+j*3+i)=*(Value+j);  
646.         }  
647.     }  
648.     printf("*****彩色图像三通道分别处理结果******************\n");  
649.     cout<<"EMEE:"<<emee/3.0<<endl;          
650.     delete Template;  
651.     delete Value;  
652.     delete Data;  
653.     delete LogImg;  
654.     return;  
655. }  

 

 

 

 

[cpp]  view plain  copy

 print ?

1. // Retinex.cpp : Defines the entry point for the console application.  
2. //  
3.   
4. #include "stdafx.h"  
5.   
6. #include "Retinex.h"  
7. //  
8. //  main  
9. //  
10. int main()  
11. {     
12.     BYTE *OrgImg,*ResImg,*p1,*p2;  
13.     int width,height,i,j,suc,area1,area2,Tsize;  
14.     bool isRGB,ret;  
15.     clock_t t1,t2;  
16.     double emee,ame;  
17.     char ch;  
18.     int *offsetdata=new int[2];  
19.     for(i=0;i<2;i++)  
20.     {  
21.         *(offsetdata+i)=0X80808080;  
22.     }      
23.     system( "cls" );      
24.     printf("******中心/环绕Retienx算法******************\n");  
25.     printf("       1.处理灰度图像\n");  
26.     printf("       2.处理彩色图像\n");  
27.     printf("*********************************************\n");  
28.     printf("请选择(1或2): ");     
29.     do  
30.     {  
31.         cin>>ch;   
32.     }while( ch != '1' && ch != '2');  
33.   
34.     //system ( "cls" );  
35.   
36.   
37.     if ( ch == '1')  
38.         isRGB=false;  
39.     else if ( ch == '2')  
40.         isRGB=true;  
41.     // open file  
42.   
43.     string image_name;  
44.     cout<<endl<<"输入图像名：";  
45.     cin>>image_name;  
46.   
47.     Mat image_dst;  
48.   
49.     if (!isRGB)  
50.     {  
51.         image_dst = imread(image_name,0);  
52.     }  
53.     else  
54.     {  
55.         image_dst = imread(image_name,1);  
56.     }  
57.   
58.       
59.       
60.     if(image_dst.empty())  
61.     {  
62.         return -1;  
63.     } //是否加载成功  
64.   
65.     imshow(image_name,image_dst);  
66.       
67.   
68.      width = image_dst.cols;    
69.      height = image_dst.rows;    
70.     int channel = image_dst.channels();   
71.      int step = width * channel* 1;  
72.     uchar* ps = NULL;  
73.       
74.     p1 = new BYTE[width*height*channel+8];  
75.   
76.     for (int i = 0; i < height; i++)    
77.     {    
78.         ps = image_dst.ptr<uchar>(i);   
79.         int bmp_i = height -1 -i;//图片上下倒置  
80.         for (int j = 0; j <width; j++)    
81.         {    
82.             if (1 == channel)    
83.             {    
84.                 *(p1 + i*step + j) = ps[j];    
85.                   
86.             }    
87.             else if (3 == channel)    
88.             {  //opencv的顺序是bgr,bmp本来是grb？  
89.                 *(p1 + bmp_i*step + j*3) = ps[j*3];  //b  
90.                 *(p1 + bmp_i*step + j*3 + 1) = ps[j*3 + 1];  //g  
91.                 *(p1 + bmp_i*step + j*3 + 2) = ps[j*3 + 2];  //r  
92.             }    
93.         }    
94.     }    
95.   
96.     BYTE* px;  
97.       
98.   
99.     if(!isRGB)  
100.         px = Read8BitBmpFile2Img(image_name.c_str(),&width,&height);  
101.     else  
102.         px = Read24BitBmpFile2Img(image_name.c_str(),&width,&height);  
103.   
104.   
105.     int len = width*height*channel;  
106.     int huanhang = width *channel;  
107.   
108.   
109.     FILE *pFile1 = fopen("out_opencv.txt",     "w"); // 文件打开方式 如果原来有内容也会销毁  
110.     FILE *pFile2 = fopen("out_ori.txt",     "w"); // 文件打开方式 如果原来有内容也会销毁  
111.   
112.     for (int i = 1;i <=len;i++)  
113.     {  
114.         fprintf(pFile1,"%d\t",(int)*(p1+i));  
115.         if (i % (huanhang)==0)  
116.         {  
117.             fprintf(pFile1,"\n");  
118.           
119.         }  
120.   
121.     }  
122.     cout<<"下面是正常的"<<endl;  
123.   
124.     for (int i = 1;i <=len;i++)  
125.     {  
126.         fprintf(pFile2,"%d\t",(int)*(px+i));  
127.   
128.         if (i % huanhang==0)  
129.         {  
130.             fprintf(pFile2,"\n");  
131.         }  
132.   
133.     }  
134.   
135.     fclose(pFile1);  
136.     fclose(pFile2);  
137.     fflush(pFile1);  
138.     fflush(pFile2);  
139.   
140.     if (p1==NULL)  
141.     {    
142.         printf("*fopen err!\n");  
143.         delete p1;  
144.         return 0;  
145.     }  
146.     if(width%64!=0||height%64!=0)  
147.     {  
148.         cout<<"图像大小需是64的倍数"<<endl;  
149.         delete p1;  
150.         return 0;  
151.     }  
152.   
153.     area1=width*height;  
154.   
155.     if(!isRGB)  
156.     {  
157.         OrgImg=p1+8-int(p1)%8;    
158.         p2=new BYTE[width*height+8];  
159.         ResImg=p2+8-int(p2)%8;  
160.         t2=clock();  
161.         GrayImageProcess(OrgImg,width,height,ResImg);  
162.         t1=clock();  
163.         printf("*****灰度图像处理结果******************\n");  
164.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
165.         Measure(ResImg,width,height,emee,ame);  
166.         cout<<"EMEE:"<<emee<<endl;  
167.         suc=Write8BitImg2BmpFile(ResImg,width,height,"result_1.bmp");  
168.     }  
169.     else  
170.     {     
171.         OrgImg=p1;//+8-int(p1)%8;     
172.         p2=new BYTE[width*height*3+8];        
173.         ResImg=p2;//+8-int(p2)%8;  
174.         t2=clock();  
175.         ColorImageProcess(OrgImg,width,height,ResImg);  
176.         t1=clock();  
177.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
178.         t2=clock();  
179.         suc=Write24BitImg2BmpFile(ResImg,width,height,"result_1.bmp");  
180.         ColorImageProcess2(OrgImg,width,height,ResImg);  
181.         t1=clock();  
182.         cout<<"运行时间："<<t1-t2<<"ms"<<endl;  
183.         suc=Write24BitImg2BmpFile(ResImg,width,height,"result_2.bmp");  
184.     }  
185.     if(suc==-1)  
186.     {  
187.         printf("*fwrite err!\n");    
188.     }  
189.   
190.     Mat result1 = imread("result_1.bmp",1);  
191.     Mat result2 = imread("result_2.bmp",1);  
192.   
193.     imshow("result1",result1);  
194.     imshow("result2",result2);  
195.     //release mem  
196.     delete p1;  
197.     delete p2;  
198.   
199.   
200.     waitKey(0);  
201.     return 0;  
202. }