00001 #ifndef IMG_FILTER_H_
00002 #define IMG_FILTER_H_
00003
00004
00005
00006
00007 #include <math.h>
00008 #include <vector>
00009 #include <algorithm>
00010
00011 namespace img {
00012
00013 template<int Channels, typename ScalarType, bool Safe>
00014 inline void normalize(Image<Channels,ScalarType,Safe> &image)
00015 {
00016 assert(image.isValid());
00017 assert(image.attributes.hasRange(0,255));
00018 if(Safe){
00019 if(!image.isValid()) throw ImageException("Invalid image");
00020 if(!image.attributes.hasRange(0,255)) throw ImageException("Invalid range attribute");
00021 }
00022 ScalarType max=maxValue(image);
00023 ScalarType min=minValue(image);
00024 ScalarType scale=255.0f/(max-min);
00025
00026 ScalarType* array = image.dataValues();
00027 int length =image.dataValuesSize();
00028
00029 for(int offset=0;offset<length;offset++)
00030 array[offset] = (array[offset]-min)*scale;
00031
00032 }
00033
00034 template<int Channels, typename ScalarType, bool Safe>
00035 inline Image<Channels,ScalarType,Safe> getNormalized(const Image<Channels,ScalarType,Safe> &image)
00036 {
00037 Image<Channels,ScalarType,Safe> i(image);
00038 normalize(i);
00039 return i;
00040 }
00041
00042 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00043 inline void convolution(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,const DestScalarType *matrix,int matrix_width,int matrix_height)
00044 {
00045 assert(source.isValid());
00046 assert(matrix != NULL);
00047 assert((matrix_width > 0) && ((matrix_width)%2 == 1));
00048 assert((matrix_height > 0) && ((matrix_height)%2 == 1));
00049 if(SrcSafe || DestSafe){
00050 if(!source.isValid()) throw ImageException("Invalid image");
00051 if( matrix == NULL) throw ImageException("NULL convolution matrix");
00052 if( !((matrix_width > 0) && ((matrix_width)%2 == 1)) ) throw ImageException("Matrix width must be a positive odd number");
00053 if( !((matrix_height > 0) && ((matrix_height)%2 == 1)) ) throw ImageException("Matrix height must be a positive odd number");
00054 }
00055 destination.setZero(source.width(),source.height());
00056 destination.attributes=source.attributes;
00057
00058 int x_radius=(matrix_width-1)/2;
00059 int y_radius=(matrix_height-1)/2;
00060
00061
00062 for(int channel=0; channel < Channels ;channel++){
00063
00064
00065 for (int y = 0; y < source.height(); y++){
00066
00067
00068 for (int x = 0; x < source.width(); x++){
00069 DestScalarType sum=0.0f;
00070 int offset=0;
00071
00072 for(int my = y-y_radius; my <= y+y_radius; my++)
00073 for(int mx = x-x_radius; mx <= x+x_radius; mx++)
00074 sum += matrix[offset++] * static_cast<DestScalarType>(source.getValueAsClamped(mx,my,channel));
00075
00076 destination.setValue(x,y,channel,sum);
00077
00078 }
00079 }
00080 }
00081 }
00082
00083 template<int Channels, typename ScalarType, bool Safe>
00084 inline Image<Channels,ScalarType,Safe> getConvolved(const Image<Channels,ScalarType,Safe> &image,const ScalarType *matrix,int matrix_width,int matrix_height)
00085 {
00086 Image<Channels,ScalarType,Safe> i;
00087 convolution(image,i,matrix,matrix_width,matrix_height);
00088 return i;
00089 }
00090
00091 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00092 inline void boxFilter(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,int radius)
00093 {
00094 assert(radius > 0);
00095 if(SrcSafe || DestSafe){
00096 if(radius <= 0) throw ImageException("Nonpositive radius");
00097 }
00098 int matrix_side=2*radius+1;
00099 int matrix_size=matrix_side*matrix_side;
00100 DestScalarType* matrix=new DestScalarType[matrix_size];
00101 DestScalarType val=1.0f/matrix_size;
00102 for(int i=0; i<matrix_size; i++)
00103 matrix[i]=val;
00104 convolution(source,destination,matrix,matrix_side,matrix_side);
00105 delete [] matrix;
00106 }
00107
00108 template<int Channels, typename ScalarType, bool Safe>
00109 inline Image<Channels,ScalarType,Safe> getBoxFiltered(const Image<Channels,ScalarType,Safe> &image,int radius)
00110 {
00111 Image<Channels,ScalarType,Safe> i;
00112 boxFilter(image,i,radius);
00113 return i;
00114 }
00115
00116 template<typename ScalarType>
00117 inline void _gaussian(const int &radius,const ScalarType &sigma,ScalarType * &matrix,int &matrix_side)
00118 {
00119 matrix_side=2*radius+1;
00120 int matrix_size=matrix_side*matrix_side;
00121 matrix=new ScalarType[matrix_size];
00122
00123
00124 int offset=0;
00125 const ScalarType c = -0.5f / (sigma*sigma);
00126
00127 for(int y = -radius;y <= radius; y++)
00128 for(int x = -radius;x <= radius; x++)
00129 matrix[offset++] = exp( (x*x+y*y)*c );
00130
00131
00132 ScalarType sum=0.0f;
00133 for(int i=0;i<matrix_size;i++)
00134 sum+=matrix[i];
00135 for(int i=0;i<matrix_size;i++)
00136 matrix[i]/=sum;
00137 }
00138
00139 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00140 inline void GaussianSmooth(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,const int radius)
00141 {
00142 assert(radius > 0.0f);
00143 if(SrcSafe || DestSafe){
00144 if(radius <= 0.0f) throw ImageException("Nonpositive radius");
00145 }
00146 const DestScalarType sigma = radius/3.0f;
00147 DestScalarType *matrix=NULL;
00148 int matrix_side=0;
00149 _gaussian(radius,sigma,matrix,matrix_side);
00150
00151 convolution(source,destination,matrix,matrix_side,matrix_side);
00152 delete [] matrix;
00153 }
00154
00155 template<int Channels, typename ScalarType, bool Safe>
00156 inline Image<Channels,ScalarType,Safe> getGaussianSmoothed(const Image<Channels,ScalarType,Safe> &image,const int radius)
00157 {
00158 Image<Channels,ScalarType,Safe> i;
00159 GaussianSmooth(image,i,radius);
00160 return i;
00161 }
00162
00163 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00164 inline void LaplacianFilter(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination)
00165 {
00166 DestScalarType *matrix=new DestScalarType[9];
00167 matrix[0]= 0.0f; matrix[1]= 1.0f; matrix[2]= 0.0f;
00168 matrix[3]= 1.0f; matrix[4]=-4.0f; matrix[5]= 1.0f;
00169 matrix[6]= 0.0f, matrix[7]= 1.0f; matrix[8]= 0.0f;
00170
00171 convolution(source,destination,matrix,3,3);
00172 delete [] matrix;
00173 }
00174
00175 template<int Channels, typename ScalarType, bool Safe>
00176 inline Image<Channels,ScalarType,Safe> getLaplacianFiltered(const Image<Channels,ScalarType,Safe> &image)
00177 {
00178 Image<Channels,ScalarType,Safe> i;
00179 LaplacianFilter(image,i);
00180 return i;
00181 }
00182
00183 template<typename ScalarType>
00184 inline void _laplacian_of_gaussian(const int &radius,const ScalarType &sigma,ScalarType * &matrix,int &matrix_side)
00185 {
00186 matrix_side=2*radius+1;
00187 int matrix_size=matrix_side*matrix_side;
00188 matrix=new ScalarType[matrix_size];
00189
00190
00191 int offset=0;
00192
00193 const ScalarType c1 = -0.5f/(sigma*sigma);
00194 ScalarType c;
00195
00196 for(int y = -radius;y <= radius; y++)
00197 for(int x = -radius;x <= radius; x++){
00198 c = (x*x+y*y) * c1;
00199 matrix[offset++] = (1+c) * exp(c);
00200 }
00201
00202
00203 ScalarType sum=0.0f;
00204 for(int i=0;i<matrix_size;i++)
00205 sum+=matrix[i];
00206 for(int i=0;i<matrix_size;i++)
00207 matrix[i]/=sum;
00208 }
00209
00210 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00211 inline void LoGFilter(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,int radius)
00212 {
00213 assert(radius > 0.0f);
00214 if(SrcSafe || DestSafe){
00215 if(radius <= 0.0f) throw ImageException("Nonpositive radius");
00216 }
00217 DestScalarType *matrix=NULL;
00218 int matrix_side=0;
00219 DestScalarType sigma = radius/3.0f;
00220 _laplacian_of_gaussian(radius,sigma,matrix,matrix_side);
00221
00222 convolution(source,destination,matrix,matrix_side,matrix_side);
00223 delete [] matrix;
00224 }
00225
00226 template<int Channels, typename ScalarType, bool Safe>
00227 inline Image<Channels,ScalarType,Safe> getLoGFiltered(const Image<Channels,ScalarType,Safe> &image,int radius)
00228 {
00229 Image<Channels,ScalarType,Safe> i;
00230 LoGFilter(image,i,radius);
00231 return i;
00232 }
00233
00234 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00235 inline void DoGFilter(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,int radius1,int radius2)
00236 {
00237 assert(radius1 > 0.0f);
00238 assert(radius2 > 0.0f);
00239 assert(radius2 > radius1);
00240 if(SrcSafe || DestSafe){
00241 if(radius1 <= 0.0f) throw ImageException("Nonpositive radius1");
00242 if(radius2 <= 0.0f) throw ImageException("Nonpositive radius2");
00243 if(radius2 <= radius1) throw ImageException("radius2 is less than radius1");
00244 }
00245
00246 int matrix_side=0;
00247
00248 DestScalarType *m1=NULL,*m2=NULL;
00249 const DestScalarType sigma1=radius1/3.0f;
00250 const DestScalarType sigma2=radius2/3.0f;
00251
00252
00253 _gaussian(radius2,sigma1,m1,matrix_side);
00254
00255 _gaussian(radius2,sigma2,m2,matrix_side);
00256
00257 int matrix_size=matrix_side*matrix_side;
00258
00259
00260
00261 for(int i=0;i<matrix_size;i++)
00262 m1[i] -= m2[i];
00263
00264
00265 DestScalarType sum=0.0f;
00266 for(int i=0;i<matrix_size;i++)
00267 sum+=m1[i];
00268 for(int i=0;i<matrix_size;i++)
00269 m1[i]/=sum;
00270
00271 convolution(source,destination,m1,matrix_side,matrix_side);
00272 delete [] m1;
00273 delete [] m2;
00274 }
00275
00276 template<int Channels, typename ScalarType, bool Safe>
00277 inline Image<Channels,ScalarType,Safe> getDoGFiltered(const Image<Channels,ScalarType,Safe> &image,int radius1,int radius2)
00278 {
00279 Image<Channels,ScalarType,Safe> i;
00280 DoGFilter(image,i,radius1,radius2);
00281 return i;
00282 }
00283
00284 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00285 inline void UnsharpMask(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,int radius,float factor)
00286 {
00287 assert(radius > 0);
00288 assert(factor > 0.0f);
00289 if(SrcSafe || DestSafe){
00290 if(radius <= 0.0f) throw ImageException("Nonpositive radius");
00291 if(factor <= 0.0f) throw ImageException("Nonpositive factor");
00292 }
00293
00294 GaussianSmooth(source,destination,radius);
00295
00296 DestScalarType* source_array = source.dataValues();
00297 DestScalarType* destination_array = destination.dataValues();
00298 int length = source.dataValuesSize();
00299
00300
00301 for(int offset=0;offset<length;offset++)
00302 destination_array[offset] = source_array[offset]+factor*(source_array[offset]-destination_array[offset]);
00303
00304 }
00305
00306 template<int Channels, typename ScalarType, bool Safe>
00307 inline Image<Channels,ScalarType,Safe> getUnsharpMasked(const Image<Channels,ScalarType,Safe> &image,int radius,float factor)
00308 {
00309 Image<Channels,ScalarType,Safe> i;
00310 UnsharpMask(image,i,radius,factor);
00311 return i;
00312 }
00313
00314 template<int Channels, typename SrcScalarType, bool SrcSafe,typename DestScalarType, bool DestSafe>
00315 inline void medianFilter(const Image<Channels,SrcScalarType,SrcSafe> &source,Image<Channels,DestScalarType,DestSafe> &destination,int radius)
00316 {
00317 assert(source.isValid());
00318 assert(radius > 0);
00319 if(SrcSafe || DestSafe){
00320 if(!source.isValid()) throw ImageException("Invalid image");
00321 if(radius <= 0) throw ImageException("Nonpositive radius");
00322 }
00323 destination.setZero(source.width(),source.height());
00324 destination.attributes=source.attributes;
00325
00326
00327 for(int channel=0; channel < Channels;channel++){
00328
00329
00330 for (int y = 0; y < source.height(); y++){
00331
00332
00333 for (int x = 0; x < source.width(); x++){
00334
00335
00336 std::vector<DestScalarType> v;
00337 for(int my = y-radius; my <= y+radius; my++)
00338 for(int mx = x-radius; mx <= x+radius; mx++)
00339 if (source.isInside(mx,my))
00340 v.push_back(static_cast<DestScalarType>(source.getValue(mx,my,channel)));
00341
00342
00343 int s=v.size();
00344 assert(s>0);
00345 nth_element (v.begin(), v.begin()+(s/2), v.end());
00346 DestScalarType median = *(v.begin()+(s/2));
00347 if((s%2)==0) {
00348 nth_element (v.begin(), v.begin()+(s/2)+1, v.end());
00349 median = ( *(v.begin()+(s/2)+1) + median ) /2.0f;
00350 }
00351
00352 destination.setValue(x,y,channel,median);
00353
00354 }
00355 }
00356 }
00357 }
00358
00359 template<int Channels, typename ScalarType, bool Safe>
00360 inline Image<Channels,ScalarType,Safe> getMedianFiltered(const Image<Channels,ScalarType,Safe> &image,int radius)
00361 {
00362 Image<Channels,ScalarType,Safe> i;
00363 medianFilter(image,i,radius);
00364 return i;
00365 }
00366
00367 template<int Channels, typename ScalarType1, bool Safe1,typename ScalarType2, bool Safe2>
00368 inline void channels_mean(const img::Image<Channels,ScalarType1,Safe1> &channels_image, img::Image<1,ScalarType2,Safe2> &mean_image)
00369 {
00370 assert(channels_image.isValid());
00371 if(Safe1 || Safe2){
00372 if(!channels_image.isValid()) throw img::ImageException("channels_image rgb image");
00373 }
00374 mean_image.setZero(channels_image.width(),channels_image.height());
00375 mean_image.attributes=channels_image.attributes;
00376
00377 for (int y_coord = 0; y_coord < channels_image.height(); ++y_coord)
00378 for (int x_coord = 0; x_coord < channels_image.width(); ++x_coord){
00379 ScalarType2 sum = ScalarType2(0.0);
00380 for (int channel=0; channel<Channels; ++channel)
00381 sum += static_cast<ScalarType2>(channels_image.getValue(x_coord, y_coord, channel));
00382 mean_image.setValue(x_coord, y_coord, 0, sum / ScalarType2(Channels) );
00383 }
00384 }
00385
00386 }
00387
00388 #endif