blur.cpp

Go to the documentation of this file.
00001 //==============================================
00002 //  copyright            : (C) 2003-2005 by Will Stokes
00003 //==============================================
00004 //  This program is free software; you can redistribute it
00005 //  and/or modify it under the terms of the GNU General
00006 //  Public License as published by the Free Software
00007 //  Foundation; either version 2 of the License, or
00008 //  (at your option) any later version.
00009 //==============================================
00010 
00011 //Systemwide includes
00012 #include <qimage.h>
00013 #include <qstring.h>
00014 #include <math.h>
00015 
00016 //Projectwide includes
00017 #include "blur.h"
00018 
00019 //----------------------------------------------
00020 // Inputs:
00021 // -------
00022 // QImage& - image to be blurred
00023 // Float sigma - amount of blurring desired
00024 //
00025 // Description:
00026 // ------------
00027 // This method blurs a color image using the approach published
00028 // by Ian T. Young and Lucas J. van Vliet, "Recursive implementation
00029 // of the Gaussian filter", 1995 Signal Processing.
00030 //
00031 // The first solution one might take to blurring an image
00032 // is to convolve a gaussian filter with an image. This is expensive since
00033 // as the gaussian filter matrix grows, processing becomes quite expensive
00034 // ( O(M^2 * N^2) where M is the width/height of the gaussian matrix and N is the 
00035 // width/height of the image. )
00036 //
00037 // An alternative is to approximate a gaussian by repeatedly convolving a kernel 
00038 // such as uniform filter. While faster, this is only an approximation of a 
00039 // gaussian, and the weights one uses for each iteration must be hand picked in
00040 // order to follow the Well's equation:
00041 // (W_1^2 + W_2^2 + W_3^2 + W_4^2 = 12*sigma^2 + 4)
00042 //
00043 // This process is tedious and limits us in applying intermediate blur amounts.
00044 //
00045 // In theory Young and Vliet's technique (which is O(4N^2) = O(N^2)) is not only
00046 // more accurate but faster as well. The only penalty we pay is the extra storage
00047 // of a row or column while that row or column is being processed. For a 
00048 // 2048x2048 image my implementation will require an extra 16kb for this extra 
00049 // storage space, a drop in the bucket compared to the memory required to hold
00050 // the entire image in float memory while being processed.
00051 //
00052 // Before being blurred color values are converted from integer [0-255] to float
00053 // [0.0-1.0] space. This is necessary to avoid compound rounding errors. In order
00054 // to minimize space required to store the float buffer, the image is processed 
00055 // in three passes for each of the red, green, and blue color channels.
00056 //
00057 // Since the code is broken up nicely, well commented, and is a direct 
00058 // implementation of Young and Vliet's techinque I point you to their paper for a
00059 // better understanding of why/how this techinque actually works. Vliet has a 
00060 // copy on his personal web site, last seen at:
00061 // http://www.ph.tn.tudelft.nl/~lucas/
00062 //----------------------------------------------
00063 
00064 void computeCoeffs( float sigma );
00065 void fillBuffer( QImage &image, int channel );
00066 void blurBuffer();
00067 
00068 void blurRow( int row );
00069 void blurColumn( int column );
00070 
00071 void blurRegionsInRow( int y );
00072 void blurRegionsInCol( int x );
00073 
00074 void resetImageData( QImage &image, int channel, bool blurEdges);
00075 
00076 float edgeValue(int x, int y);
00077 
00078 float q, b0, b1, b2, b3, B;
00079 int width, height;
00080 float* buffer;
00081 float* rowBuffer;
00082 float* colBuffer;
00083 
00084 float* regionRowBuffer;
00085 float* regionColBuffer;
00086 
00087 QImage* edgeImage;
00088 int* regionMap;
00089 int regionCount;
00090 QPoint displayOffset;
00091 QSize fullRes;
00092 
00093 //==============================================
00094 void blurImage( QImage &image, float sigma )
00095 {
00096   //supply dummy data for edges, notably NULL for the edge image pointer.
00097   //other values have no effect
00098   blurImage( image, sigma, QPoint(0,0), image.size(), NULL, NULL, 0, false );
00099 }
00100 //==============================================
00101 void blurImage( QImage &image, float sigma,
00102                 QPoint offset, QSize fullImageRes,
00103                 QImage* edges, int* regions, int numRegions,
00104                 bool targetEdges)
00105 {
00106   edgeImage = edges;
00107   regionMap = regions;
00108   regionCount = numRegions;
00109   displayOffset = offset;
00110   fullRes = fullImageRes;
00111   
00112   //compute blurring coeffecients
00113   computeCoeffs(sigma);
00114   
00115   //store image dimensions
00116   width = image.width();
00117   height = image.height();
00118   
00119   //Construct float buffer that is the size of the image/
00120   //In order to conserve memory process image three times, once for
00121   //each color channel.
00122   buffer = new float[ width * height ];
00123 
00124   rowBuffer = new float[width];
00125   colBuffer = new float[height];
00126   
00127   regionRowBuffer = new float[width * numRegions];
00128   regionColBuffer = new float[height * numRegions];
00129   
00130   //iterate over each color channel
00131   int channel;
00132   for( channel = 0; channel <=2; channel++)
00133   {
00134     //copy color data into float buffer
00135     fillBuffer( image, channel );
00136     
00137     //blur buffer data
00138     blurBuffer();
00139     
00140     //reset image data used blurred buffer
00141     resetImageData(image, channel, targetEdges);
00142   }
00143   
00144   //delete buffer
00145   delete[] buffer;
00146   delete[] rowBuffer;
00147   delete[] colBuffer;
00148 }
00149 //==============================================
00150 void computeCoeffs( float sigma )
00151 {
00152   //compute q as a function of sigma
00153   if( sigma >= 2.5f )
00154   {
00155     q = 0.98711f*sigma - 0.96330f;  
00156   }
00157   else
00158   {
00159     q = 3.97156f - 4.14554f * sqrt( 1.0f - 0.26891f*sigma );
00160   }
00161 
00162   //compute b0, b1, b2, and b3
00163   b0 = 1.57825f + (2.44413f*q) + (1.4281f * q*q ) + (0.422205f * q*q*q );
00164   b1 = (2.44413f * q) + (2.85619f * q*q) + (1.26661 * q*q*q );
00165   b2 = -((1.4281 * q*q) + (1.26661 * q*q*q));
00166   b3 = 0.422205 * q*q*q;
00167   
00168   //compute B
00169   B = 1.0f - ((b1 + b2 + b3) / b0); 
00170 }
00171 //==============================================
00172 void fillBuffer( QImage &image, int channel )
00173 {
00174   //precompute 1/255
00175   float multiplier = 1.0f / 255.0f;
00176   
00177   //iterate over each selected scanline 
00178   int x, y;
00179   QRgb* rgb;
00180   uchar* scanLine;
00181   for( y=0; y<image.height(); y++)
00182   {   
00183     //iterate over each pixel in scanline
00184     scanLine = image.scanLine(y);
00185     for( x=0; x<image.width(); x++)
00186     {
00187       //get handle on rgb value in image
00188       rgb = ((QRgb*)scanLine+x);
00189 
00190       //compute index where float value is stored in buffer
00191       int index = x + y*image.width();
00192       
00193       //convert and store correct channel in buffer
00194       if( channel == 0 )
00195         buffer[index] = multiplier * qRed( *rgb );
00196       else if( channel == 1 )
00197         buffer[index] = multiplier * qGreen( *rgb );
00198       else
00199         buffer[index] = multiplier * qBlue( *rgb );
00200     } //x
00201   } //y
00202 }
00203 //==============================================
00204 void blurBuffer()
00205 {
00206   //blur rows, then columns
00207   int index;
00208   
00209   if(regionMap == NULL || edgeImage == NULL )
00210   {
00211     for(index=0; index < height; index++)
00212     { blurRow( index ); }
00213     
00214     for(index=0; index< width; index++)
00215     { blurColumn( index ); }
00216   }
00217   else
00218   {
00219     for(index=0; index < height; index++)
00220     { blurRegionsInRow( index ); }
00221     
00222     for(index=0; index< width; index++)
00223     { blurRegionsInCol( index ); }
00224   }
00225 }
00226 //==============================================
00227 int regionIndex(int x, int y)
00228 {
00229   int edgeX = ((edgeImage->width()-1) * (x+displayOffset.x())) / (fullRes.width()-1);
00230   int edgeY = ((edgeImage->height()-1) * (y+displayOffset.y())) / (fullRes.height()-1);
00231   return edgeY*edgeImage->width() + edgeX;
00232 }
00233 //==============================================
00234 float edgeValue(int x, int y)
00235 {
00236   //compute floating point x and y coordinates for edge image
00237   float edgeX = ((edgeImage->width()-1.0f) * (x+displayOffset.x())) / (fullRes.width()-1);
00238   float edgeY = ((edgeImage->height()-1.0f) * (y+displayOffset.y())) / (fullRes.height()-1);
00239   
00240   //compute 4 int values of coordinates
00241   int x1 = (int)edgeX;
00242   int y1 = (int)edgeY;
00243   int x2, y2;
00244   if( edgeX > x1 )
00245     x2 = x1+1;
00246   else
00247     x2 = x1;
00248   if( edgeY > y1 )
00249     y2 = y1+1;
00250   else
00251     y2 = y1;
00252   
00253   //compute the four indices
00254   int index1, index2, index3, index4;
00255   index1 = x1 + y1*edgeImage->width();
00256   index2 = x2 + y1*edgeImage->width();
00257   index3 = x1 + y2*edgeImage->width();
00258   index4 = x2 + y2*edgeImage->width();
00259   
00260   //find edge quantity for each corner
00261   float v1, v2, v3, v4; 
00262   uchar* scanline = edgeImage->scanLine( y1 );
00263   QRgb* rgb = ((QRgb*)scanline+x1);
00264   v1 = ((float) qRed( *rgb )) / 255.0f;
00265   rgb = ((QRgb*)scanline+x2);
00266   v2 = ((float) qRed( *rgb )) / 255.0f;
00267 
00268   scanline = edgeImage->scanLine( y2 );
00269   rgb = ((QRgb*)scanline+x1);
00270   v3 = ((float) qRed( *rgb )) / 255.0f;
00271   rgb = ((QRgb*)scanline+x2);
00272   v4 = ((float) qRed( *rgb )) / 255.0f;
00273   
00274   //blur combine left-right
00275   v1 = (edgeX-x1)*v2 + (1 - edgeX + x1)*v1;
00276   v3 = (edgeX-x1)*v4 + (1 - edgeX + x1)*v3;
00277   
00278   //combine top-bottom
00279   v1 = (edgeY-y1)*v3 + (1 - edgeY + y1)*v1;
00280   
00281   //return result
00282   return v1;
00283 }
00284 //==============================================
00285 void blurRow( int row )
00286 {
00287   int i;
00288   int rtw = row*width;
00289   
00290   //forward
00291   rowBuffer[0] = buffer[ 0 + rtw ];
00292   for(i=1; i<width; i++)
00293   {
00294      rowBuffer[i] = B*buffer[ i + rtw ] +
00295      ( b1*rowBuffer[ QMAX(i-1, 0) ] +
00296        b2 * rowBuffer[ QMAX(i-2, 0) ] +
00297        b3 * rowBuffer[ QMAX(i-3, 0) ]) / b0;    
00298     
00299   }
00300   
00301   //reverse
00302   for(i=width-1; i>=0; i--)
00303   {
00304      buffer[ i + rtw ] = B*rowBuffer[ i ] +
00305      ( b1 * buffer[ QMIN(i+1, width-1) + rtw ] +
00306        b2 * buffer[ QMIN(i+2, width-1) + rtw ] +
00307        b3 * buffer[ QMIN(i+3, width-1) + rtw ]) / b0;    
00308   }
00309 }
00310 //==============================================
00311 void blurRegionsInRow( int y )
00312 {
00313   //---------------------------------
00314   //populate region row buffer. a row has been allocated for 
00315   //each region. Pixels between regions
00316   //take the closest pixel value in that row from that region
00317   int yTimesWidth = y*width;
00318   int regionTimesWidth;
00319   int region,x,x2;
00320   
00321   //for each region
00322   for(region=0; region<regionCount; region++)
00323   {
00324     regionTimesWidth = region*width;
00325     int lastX = -1;
00326     for(x=0; x<width; x++)
00327     {
00328       //if pixel belongs to this region then update lastX index and copy value over
00329       //if lastX is mroe than one pixel away than fill inbetween region
00330       if( region == regionMap[regionIndex(x, y)] )
00331       {
00332         //fill empty region preceeding this region blob
00333         if( lastX < x-1)
00334         {
00335           //no preceeding region, spread left!
00336           if(lastX == -1)
00337           {
00338             for(x2=0; x2<x; x2++) { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
00339           }
00340           //else spread from both left and right of empty stretch
00341           else
00342           {
00343             int xMid = lastX + ((x-1) - lastX)/2;
00344             
00345             for(x2=lastX+1; x2<=xMid; x2++)
00346             { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }
00347 
00348             for(x2=xMid+1; x2<x; x2++)
00349             { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
00350           }
00351         }
00352         
00353         regionRowBuffer[x + regionTimesWidth] = buffer[x + yTimesWidth];
00354         lastX = x;
00355       } 
00356     } //x
00357     
00358     //if last stretch is empty, fill right
00359     if( region != regionMap[regionIndex(width-1, y)] )
00360     {
00361       for(x2=lastX+1; x2<width; x2++)
00362       { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }
00363     }
00364     
00365   } //region
00366   //---------------------------------
00367   //blur the region row buffers
00368 
00369   //for each region
00370   for(region=0; region<regionCount; region++)
00371   {
00372     regionTimesWidth = region*width;
00373     
00374     //forward
00375     rowBuffer[0] = regionRowBuffer[ 0 + regionTimesWidth ];
00376     for(x=1; x<width; x++)
00377     {
00378       rowBuffer[x] = B*regionRowBuffer[ x + regionTimesWidth ] +
00379       ( b1*rowBuffer[ QMAX(x-1, 0) ] +
00380         b2 * rowBuffer[ QMAX(x-2, 0) ] +
00381         b3 * rowBuffer[ QMAX(x-3, 0) ]) / b0;    
00382     }
00383     
00384     //reverse
00385     for(x=width-1; x>=0; x--)
00386     {
00387       regionRowBuffer[ x + regionTimesWidth ] = B*rowBuffer[ x ] +
00388       ( b1 * regionRowBuffer[ QMIN(x+1, width-1) + regionTimesWidth ] +
00389         b2 * regionRowBuffer[ QMIN(x+2, width-1) + regionTimesWidth ] +
00390         b3 * regionRowBuffer[ QMIN(x+3, width-1) + regionTimesWidth ]) / b0;    
00391     }
00392   }
00393   //---------------------------------
00394   //copy data from the region row buffers back to the
00395   //buffer. for each pixel we choose the correct region
00396   //row buffer basedon the original regionidentity of hte pixel
00397   for(x=0; x<width; x++)
00398   {
00399     int ri = regionIndex(x,y);
00400     int region = regionMap[ri];
00401     float bufferVal = regionRowBuffer[ x + region*width ];
00402     buffer[x + yTimesWidth] = bufferVal;
00403     
00404 //    buffer[x + yTimesWidth] = regionRowBuffer[x + regionMap[regionIndex(x,y)]*width]; 
00405   }
00406   //---------------------------------
00407 }
00408 //==============================================
00409 void blurColumn( int column )
00410 {
00411   int i;
00412 
00413   //forward
00414   colBuffer[0] = buffer[ column + 0*width ];
00415   for(i=1; i<height; i++)
00416   {
00417      colBuffer[i] = B*buffer[ column + i*width ] +
00418      ( b1 * colBuffer[ QMAX(i-1, 0) ] +
00419        b2 * colBuffer[ QMAX(i-2, 0) ] +
00420        b3 * colBuffer[ QMAX(i-3, 0) ]) / b0;    
00421   }
00422   
00423   //reverse
00424   for(i=height-1; i>=0; i--)
00425   {
00426      buffer[ column + i*width ] = B*colBuffer[ i ] +
00427      ( b1 * buffer[ column + QMIN(i+1, height-1)*width ] +
00428        b2 * buffer[ column + QMIN(i+2, height-1)*width ] +
00429        b3 * buffer[ column + QMIN(i+3, height-1)*width ]) / b0;        
00430   }
00431   
00432 }
00433 //==============================================
00434 void blurRegionsInCol( int x )
00435 {
00436   //---------------------------------
00437   //populate region col buffer. a col has been allocated for 
00438   //each region. Pixels between regions
00439   //take the closest pixel value in that col from that region
00440 //  int yTimesWidth = y*width;
00441   int regionTimesHeight;
00442   int region,y,y2;
00443   
00444   //for each region
00445   for(region=0; region<regionCount; region++)
00446   {
00447     regionTimesHeight = region*height;
00448     int lastY = -1;
00449     for(y=0; y<height; y++)
00450     {
00451       //if pixel belongs to this region then update lastY index and copy value over
00452       //if lastY is more than one pixel away than fill inbetween region
00453       if( region == regionMap[regionIndex(x, y)] )
00454       {
00455         //fill empty region preceeding this region blob
00456         if( lastY < y-1)
00457         {
00458           //no preceeding region, spread left!
00459           if(lastY == -1)
00460           {
00461             for(y2=0; y2<y; y2++) { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
00462           }
00463           //else spread from both left and right of empty stretch
00464           else
00465           {
00466             int yMid = lastY + ((y-1) - lastY)/2;
00467             
00468             for(y2=lastY+1; y2<=yMid; y2++)
00469             { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
00470             
00471             for(y2=yMid+1; y2<y; y2++)
00472             { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
00473           }
00474         }
00475         
00476         regionColBuffer[y + regionTimesHeight] = buffer[x + y*width];
00477         lastY = y;
00478       } 
00479     } //y
00480     
00481     //if last stretch is empty, fill right
00482     if( region != regionMap[regionIndex(x, height-1)] )
00483     {
00484       for(y2=lastY+1; y2<height; y2++)
00485       { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
00486     }
00487     
00488   } //region
00489   //---------------------------------
00490   //blur the region col buffers
00491   
00492   //for each region
00493   for(region=0; region<regionCount; region++)
00494   {
00495     regionTimesHeight = region*height;
00496     
00497     //forward
00498     colBuffer[0] = regionColBuffer[ 0 + regionTimesHeight ];
00499     for(y=1; y<height; y++)
00500     {
00501       colBuffer[y] = B*regionColBuffer[ y + regionTimesHeight ] +
00502       ( b1 * colBuffer[ QMAX(y-1, 0) ] +
00503         b2 * colBuffer[ QMAX(y-2, 0) ] +
00504         b3 * colBuffer[ QMAX(y-3, 0) ]) / b0;    
00505     }
00506     
00507     //reverse
00508     for(y=height-1; y>=0; y--)
00509     {
00510       regionColBuffer[ y + regionTimesHeight ] = B*colBuffer[ y ] +
00511       ( b1 * regionColBuffer[ QMIN(y+1, height-1) + regionTimesHeight ] +
00512         b2 * regionColBuffer[ QMIN(y+2, height-1) + regionTimesHeight ] +
00513         b3 * regionColBuffer[ QMIN(y+3, height-1) + regionTimesHeight ]) / b0;    
00514     }
00515   }
00516   //---------------------------------
00517   //copy data from the region row buffers back to the
00518   //buffer. for each pixel we choose the correct region
00519   //row buffer basedon the original regionidentity of hte pixel
00520   for(y=0; y<height; y++)
00521   {
00522     buffer[x + y*width] = regionColBuffer[y + regionMap[regionIndex(x,y)]*height]; 
00523   }
00524   //---------------------------------
00525 }
00526 //==============================================
00527 void resetImageData( QImage &image, int channel, bool blurEdges)
00528 {
00529   //iterate over each selected scanline 
00530   int x, y;
00531   QRgb *rgb;
00532   uchar* imageScanline = NULL;
00533   for( y=0; y<image.height(); y++)
00534   {   
00535     imageScanline = image.scanLine(y);
00536     for( x=0; x<image.width(); x++)
00537     {
00538       //get handle on rgb value in image
00539       rgb = ((QRgb*)imageScanline+x);
00540       
00541       //compute index where float value is stored in buffer
00542       int index = x + y*image.width();
00543 
00544       //convert blured value to 0-255 range
00545       int blurredColor = QMAX( QMIN( ((int) (255*buffer[index])), 255 ), 0 );
00546       
00547       //blur the entire thing!
00548       float alpha;
00549       if( edgeImage == NULL)
00550         alpha = 1.0f;
00551       else
00552       {
00553         alpha = edgeValue( x, y ); 
00554         if(!blurEdges)
00555           alpha = 1.0f - alpha;
00556       }
00557 
00558       //convert and store correct channel in buffer
00559       if( channel == 0 )
00560         *rgb = qRgb( (int) (alpha*blurredColor + (1-alpha)*qRed(*rgb)),
00561                      qGreen(*rgb), qBlue(*rgb) );
00562       else if( channel == 1 )
00563         *rgb = qRgb( qRed(*rgb), 
00564                      (int) (alpha*blurredColor + (1-alpha)*qGreen(*rgb)),
00565                      qBlue(*rgb) );
00566       else
00567         *rgb = qRgb( qRed(*rgb), qGreen(*rgb), 
00568                      (int) (alpha*blurredColor + (1-alpha)*qBlue(*rgb)) );
00569     } //x
00570   } //y
00571 }
00572 //==============================================

Generated on Thu Jan 3 10:52:45 2008 for AlbumShaper by  doxygen 1.5.4