AlbumShaper  1.0a3
blur.cpp
Go to the documentation of this file.
1 //==============================================
2 // copyright : (C) 2003-2005 by Will Stokes
3 //==============================================
4 // This program is free software; you can redistribute it
5 // and/or modify it under the terms of the GNU General
6 // Public License as published by the Free Software
7 // Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9 //==============================================
10 
11 //Systemwide includes
12 #include <qimage.h>
13 #include <qstring.h>
14 #include <math.h>
15 
16 //Projectwide includes
17 #include "blur.h"
18 
19 //----------------------------------------------
20 // Inputs:
21 // -------
22 // QImage& - image to be blurred
23 // Float sigma - amount of blurring desired
24 //
25 // Description:
26 // ------------
27 // This method blurs a color image using the approach published
28 // by Ian T. Young and Lucas J. van Vliet, "Recursive implementation
29 // of the Gaussian filter", 1995 Signal Processing.
30 //
31 // The first solution one might take to blurring an image
32 // is to convolve a gaussian filter with an image. This is expensive since
33 // as the gaussian filter matrix grows, processing becomes quite expensive
34 // ( O(M^2 * N^2) where M is the width/height of the gaussian matrix and N is the
35 // width/height of the image. )
36 //
37 // An alternative is to approximate a gaussian by repeatedly convolving a kernel
38 // such as uniform filter. While faster, this is only an approximation of a
39 // gaussian, and the weights one uses for each iteration must be hand picked in
40 // order to follow the Well's equation:
41 // (W_1^2 + W_2^2 + W_3^2 + W_4^2 = 12*sigma^2 + 4)
42 //
43 // This process is tedious and limits us in applying intermediate blur amounts.
44 //
45 // In theory Young and Vliet's technique (which is O(4N^2) = O(N^2)) is not only
46 // more accurate but faster as well. The only penalty we pay is the extra storage
47 // of a row or column while that row or column is being processed. For a
48 // 2048x2048 image my implementation will require an extra 16kb for this extra
49 // storage space, a drop in the bucket compared to the memory required to hold
50 // the entire image in float memory while being processed.
51 //
52 // Before being blurred color values are converted from integer [0-255] to float
53 // [0.0-1.0] space. This is necessary to avoid compound rounding errors. In order
54 // to minimize space required to store the float buffer, the image is processed
55 // in three passes for each of the red, green, and blue color channels.
56 //
57 // Since the code is broken up nicely, well commented, and is a direct
58 // implementation of Young and Vliet's techinque I point you to their paper for a
59 // better understanding of why/how this techinque actually works. Vliet has a
60 // copy on his personal web site, last seen at:
61 // http://www.ph.tn.tudelft.nl/~lucas/
62 //----------------------------------------------
63 
64 void computeCoeffs( float sigma );
65 void fillBuffer( QImage &image, int channel );
66 void blurBuffer();
67 
68 void blurRow( int row );
69 void blurColumn( int column );
70 
71 void blurRegionsInRow( int y );
72 void blurRegionsInCol( int x );
73 
74 void resetImageData( QImage &image, int channel, bool blurEdges);
75 
76 float edgeValue(int x, int y);
77 
78 float q, b0, b1, b2, b3, B;
80 float* buffer;
81 float* rowBuffer;
82 float* colBuffer;
83 
86 
87 QImage* edgeImage;
88 int* regionMap;
91 QSize fullRes;
92 
93 //==============================================
94 void blurImage( QImage &image, float sigma )
95 {
96  //supply dummy data for edges, notably NULL for the edge image pointer.
97  //other values have no effect
98  blurImage( image, sigma, QPoint(0,0), image.size(), NULL, NULL, 0, false );
99 }
100 //==============================================
101 void blurImage( QImage &image, float sigma,
102  QPoint offset, QSize fullImageRes,
103  QImage* edges, int* regions, int numRegions,
104  bool targetEdges)
105 {
106  edgeImage = edges;
107  regionMap = regions;
108  regionCount = numRegions;
109  displayOffset = offset;
110  fullRes = fullImageRes;
111 
112  //compute blurring coeffecients
113  computeCoeffs(sigma);
114 
115  //store image dimensions
116  width = image.width();
117  height = image.height();
118 
119  //Construct float buffer that is the size of the image/
120  //In order to conserve memory process image three times, once for
121  //each color channel.
122  buffer = new float[ width * height ];
123 
124  rowBuffer = new float[width];
125  colBuffer = new float[height];
126 
127  regionRowBuffer = new float[width * numRegions];
128  regionColBuffer = new float[height * numRegions];
129 
130  //iterate over each color channel
131  int channel;
132  for( channel = 0; channel <=2; channel++)
133  {
134  //copy color data into float buffer
135  fillBuffer( image, channel );
136 
137  //blur buffer data
138  blurBuffer();
139 
140  //reset image data used blurred buffer
141  resetImageData(image, channel, targetEdges);
142  }
143 
144  //delete buffer
145  delete[] buffer;
146  delete[] rowBuffer;
147  delete[] colBuffer;
148 }
149 //==============================================
150 void computeCoeffs( float sigma )
151 {
152  //compute q as a function of sigma
153  if( sigma >= 2.5f )
154  {
155  q = 0.98711f*sigma - 0.96330f;
156  }
157  else
158  {
159  q = 3.97156f - 4.14554f * sqrt( 1.0f - 0.26891f*sigma );
160  }
161 
162  //compute b0, b1, b2, and b3
163  b0 = 1.57825f + (2.44413f*q) + (1.4281f * q*q ) + (0.422205f * q*q*q );
164  b1 = (2.44413f * q) + (2.85619f * q*q) + (1.26661 * q*q*q );
165  b2 = -((1.4281 * q*q) + (1.26661 * q*q*q));
166  b3 = 0.422205 * q*q*q;
167 
168  //compute B
169  B = 1.0f - ((b1 + b2 + b3) / b0);
170 }
171 //==============================================
172 void fillBuffer( QImage &image, int channel )
173 {
174  //precompute 1/255
175  float multiplier = 1.0f / 255.0f;
176 
177  //iterate over each selected scanline
178  int x, y;
179  QRgb* rgb;
180  uchar* scanLine;
181  for( y=0; y<image.height(); y++)
182  {
183  //iterate over each pixel in scanline
184  scanLine = image.scanLine(y);
185  for( x=0; x<image.width(); x++)
186  {
187  //get handle on rgb value in image
188  rgb = ((QRgb*)scanLine+x);
189 
190  //compute index where float value is stored in buffer
191  int index = x + y*image.width();
192 
193  //convert and store correct channel in buffer
194  if( channel == 0 )
195  buffer[index] = multiplier * qRed( *rgb );
196  else if( channel == 1 )
197  buffer[index] = multiplier * qGreen( *rgb );
198  else
199  buffer[index] = multiplier * qBlue( *rgb );
200  } //x
201  } //y
202 }
203 //==============================================
205 {
206  //blur rows, then columns
207  int index;
208 
209  if(regionMap == NULL || edgeImage == NULL )
210  {
211  for(index=0; index < height; index++)
212  { blurRow( index ); }
213 
214  for(index=0; index< width; index++)
215  { blurColumn( index ); }
216  }
217  else
218  {
219  for(index=0; index < height; index++)
220  { blurRegionsInRow( index ); }
221 
222  for(index=0; index< width; index++)
223  { blurRegionsInCol( index ); }
224  }
225 }
226 //==============================================
227 int regionIndex(int x, int y)
228 {
229  int edgeX = ((edgeImage->width()-1) * (x+displayOffset.x())) / (fullRes.width()-1);
230  int edgeY = ((edgeImage->height()-1) * (y+displayOffset.y())) / (fullRes.height()-1);
231  return edgeY*edgeImage->width() + edgeX;
232 }
233 //==============================================
234 float edgeValue(int x, int y)
235 {
236  //compute floating point x and y coordinates for edge image
237  float edgeX = ((edgeImage->width()-1.0f) * (x+displayOffset.x())) / (fullRes.width()-1);
238  float edgeY = ((edgeImage->height()-1.0f) * (y+displayOffset.y())) / (fullRes.height()-1);
239 
240  //compute 4 int values of coordinates
241  int x1 = (int)edgeX;
242  int y1 = (int)edgeY;
243  int x2, y2;
244  if( edgeX > x1 )
245  x2 = x1+1;
246  else
247  x2 = x1;
248  if( edgeY > y1 )
249  y2 = y1+1;
250  else
251  y2 = y1;
252 
253  //compute the four indices
254  int index1, index2, index3, index4;
255  index1 = x1 + y1*edgeImage->width();
256  index2 = x2 + y1*edgeImage->width();
257  index3 = x1 + y2*edgeImage->width();
258  index4 = x2 + y2*edgeImage->width();
259 
260  //find edge quantity for each corner
261  float v1, v2, v3, v4;
262  uchar* scanline = edgeImage->scanLine( y1 );
263  QRgb* rgb = ((QRgb*)scanline+x1);
264  v1 = ((float) qRed( *rgb )) / 255.0f;
265  rgb = ((QRgb*)scanline+x2);
266  v2 = ((float) qRed( *rgb )) / 255.0f;
267 
268  scanline = edgeImage->scanLine( y2 );
269  rgb = ((QRgb*)scanline+x1);
270  v3 = ((float) qRed( *rgb )) / 255.0f;
271  rgb = ((QRgb*)scanline+x2);
272  v4 = ((float) qRed( *rgb )) / 255.0f;
273 
274  //blur combine left-right
275  v1 = (edgeX-x1)*v2 + (1 - edgeX + x1)*v1;
276  v3 = (edgeX-x1)*v4 + (1 - edgeX + x1)*v3;
277 
278  //combine top-bottom
279  v1 = (edgeY-y1)*v3 + (1 - edgeY + y1)*v1;
280 
281  //return result
282  return v1;
283 }
284 //==============================================
285 void blurRow( int row )
286 {
287  int i;
288  int rtw = row*width;
289 
290  //forward
291  rowBuffer[0] = buffer[ 0 + rtw ];
292  for(i=1; i<width; i++)
293  {
294  rowBuffer[i] = B*buffer[ i + rtw ] +
295  ( b1*rowBuffer[ QMAX(i-1, 0) ] +
296  b2 * rowBuffer[ QMAX(i-2, 0) ] +
297  b3 * rowBuffer[ QMAX(i-3, 0) ]) / b0;
298 
299  }
300 
301  //reverse
302  for(i=width-1; i>=0; i--)
303  {
304  buffer[ i + rtw ] = B*rowBuffer[ i ] +
305  ( b1 * buffer[ QMIN(i+1, width-1) + rtw ] +
306  b2 * buffer[ QMIN(i+2, width-1) + rtw ] +
307  b3 * buffer[ QMIN(i+3, width-1) + rtw ]) / b0;
308  }
309 }
310 //==============================================
311 void blurRegionsInRow( int y )
312 {
313  //---------------------------------
314  //populate region row buffer. a row has been allocated for
315  //each region. Pixels between regions
316  //take the closest pixel value in that row from that region
317  int yTimesWidth = y*width;
318  int regionTimesWidth;
319  int region,x,x2;
320 
321  //for each region
322  for(region=0; region<regionCount; region++)
323  {
324  regionTimesWidth = region*width;
325  int lastX = -1;
326  for(x=0; x<width; x++)
327  {
328  //if pixel belongs to this region then update lastX index and copy value over
329  //if lastX is mroe than one pixel away than fill inbetween region
330  if( region == regionMap[regionIndex(x, y)] )
331  {
332  //fill empty region preceeding this region blob
333  if( lastX < x-1)
334  {
335  //no preceeding region, spread left!
336  if(lastX == -1)
337  {
338  for(x2=0; x2<x; x2++) { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
339  }
340  //else spread from both left and right of empty stretch
341  else
342  {
343  int xMid = lastX + ((x-1) - lastX)/2;
344 
345  for(x2=lastX+1; x2<=xMid; x2++)
346  { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }
347 
348  for(x2=xMid+1; x2<x; x2++)
349  { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
350  }
351  }
352 
353  regionRowBuffer[x + regionTimesWidth] = buffer[x + yTimesWidth];
354  lastX = x;
355  }
356  } //x
357 
358  //if last stretch is empty, fill right
359  if( region != regionMap[regionIndex(width-1, y)] )
360  {
361  for(x2=lastX+1; x2<width; x2++)
362  { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }
363  }
364 
365  } //region
366  //---------------------------------
367  //blur the region row buffers
368 
369  //for each region
370  for(region=0; region<regionCount; region++)
371  {
372  regionTimesWidth = region*width;
373 
374  //forward
375  rowBuffer[0] = regionRowBuffer[ 0 + regionTimesWidth ];
376  for(x=1; x<width; x++)
377  {
378  rowBuffer[x] = B*regionRowBuffer[ x + regionTimesWidth ] +
379  ( b1*rowBuffer[ QMAX(x-1, 0) ] +
380  b2 * rowBuffer[ QMAX(x-2, 0) ] +
381  b3 * rowBuffer[ QMAX(x-3, 0) ]) / b0;
382  }
383 
384  //reverse
385  for(x=width-1; x>=0; x--)
386  {
387  regionRowBuffer[ x + regionTimesWidth ] = B*rowBuffer[ x ] +
388  ( b1 * regionRowBuffer[ QMIN(x+1, width-1) + regionTimesWidth ] +
389  b2 * regionRowBuffer[ QMIN(x+2, width-1) + regionTimesWidth ] +
390  b3 * regionRowBuffer[ QMIN(x+3, width-1) + regionTimesWidth ]) / b0;
391  }
392  }
393  //---------------------------------
394  //copy data from the region row buffers back to the
395  //buffer. for each pixel we choose the correct region
396  //row buffer basedon the original regionidentity of hte pixel
397  for(x=0; x<width; x++)
398  {
399  int ri = regionIndex(x,y);
400  int region = regionMap[ri];
401  float bufferVal = regionRowBuffer[ x + region*width ];
402  buffer[x + yTimesWidth] = bufferVal;
403 
404 // buffer[x + yTimesWidth] = regionRowBuffer[x + regionMap[regionIndex(x,y)]*width];
405  }
406  //---------------------------------
407 }
408 //==============================================
409 void blurColumn( int column )
410 {
411  int i;
412 
413  //forward
414  colBuffer[0] = buffer[ column + 0*width ];
415  for(i=1; i<height; i++)
416  {
417  colBuffer[i] = B*buffer[ column + i*width ] +
418  ( b1 * colBuffer[ QMAX(i-1, 0) ] +
419  b2 * colBuffer[ QMAX(i-2, 0) ] +
420  b3 * colBuffer[ QMAX(i-3, 0) ]) / b0;
421  }
422 
423  //reverse
424  for(i=height-1; i>=0; i--)
425  {
426  buffer[ column + i*width ] = B*colBuffer[ i ] +
427  ( b1 * buffer[ column + QMIN(i+1, height-1)*width ] +
428  b2 * buffer[ column + QMIN(i+2, height-1)*width ] +
429  b3 * buffer[ column + QMIN(i+3, height-1)*width ]) / b0;
430  }
431 
432 }
433 //==============================================
434 void blurRegionsInCol( int x )
435 {
436  //---------------------------------
437  //populate region col buffer. a col has been allocated for
438  //each region. Pixels between regions
439  //take the closest pixel value in that col from that region
440 // int yTimesWidth = y*width;
441  int regionTimesHeight;
442  int region,y,y2;
443 
444  //for each region
445  for(region=0; region<regionCount; region++)
446  {
447  regionTimesHeight = region*height;
448  int lastY = -1;
449  for(y=0; y<height; y++)
450  {
451  //if pixel belongs to this region then update lastY index and copy value over
452  //if lastY is more than one pixel away than fill inbetween region
453  if( region == regionMap[regionIndex(x, y)] )
454  {
455  //fill empty region preceeding this region blob
456  if( lastY < y-1)
457  {
458  //no preceeding region, spread left!
459  if(lastY == -1)
460  {
461  for(y2=0; y2<y; y2++) { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
462  }
463  //else spread from both left and right of empty stretch
464  else
465  {
466  int yMid = lastY + ((y-1) - lastY)/2;
467 
468  for(y2=lastY+1; y2<=yMid; y2++)
469  { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
470 
471  for(y2=yMid+1; y2<y; y2++)
472  { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
473  }
474  }
475 
476  regionColBuffer[y + regionTimesHeight] = buffer[x + y*width];
477  lastY = y;
478  }
479  } //y
480 
481  //if last stretch is empty, fill right
482  if( region != regionMap[regionIndex(x, height-1)] )
483  {
484  for(y2=lastY+1; y2<height; y2++)
485  { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
486  }
487 
488  } //region
489  //---------------------------------
490  //blur the region col buffers
491 
492  //for each region
493  for(region=0; region<regionCount; region++)
494  {
495  regionTimesHeight = region*height;
496 
497  //forward
498  colBuffer[0] = regionColBuffer[ 0 + regionTimesHeight ];
499  for(y=1; y<height; y++)
500  {
501  colBuffer[y] = B*regionColBuffer[ y + regionTimesHeight ] +
502  ( b1 * colBuffer[ QMAX(y-1, 0) ] +
503  b2 * colBuffer[ QMAX(y-2, 0) ] +
504  b3 * colBuffer[ QMAX(y-3, 0) ]) / b0;
505  }
506 
507  //reverse
508  for(y=height-1; y>=0; y--)
509  {
510  regionColBuffer[ y + regionTimesHeight ] = B*colBuffer[ y ] +
511  ( b1 * regionColBuffer[ QMIN(y+1, height-1) + regionTimesHeight ] +
512  b2 * regionColBuffer[ QMIN(y+2, height-1) + regionTimesHeight ] +
513  b3 * regionColBuffer[ QMIN(y+3, height-1) + regionTimesHeight ]) / b0;
514  }
515  }
516  //---------------------------------
517  //copy data from the region row buffers back to the
518  //buffer. for each pixel we choose the correct region
519  //row buffer basedon the original regionidentity of hte pixel
520  for(y=0; y<height; y++)
521  {
523  }
524  //---------------------------------
525 }
526 //==============================================
527 void resetImageData( QImage &image, int channel, bool blurEdges)
528 {
529  //iterate over each selected scanline
530  int x, y;
531  QRgb *rgb;
532  uchar* imageScanline = NULL;
533  for( y=0; y<image.height(); y++)
534  {
535  imageScanline = image.scanLine(y);
536  for( x=0; x<image.width(); x++)
537  {
538  //get handle on rgb value in image
539  rgb = ((QRgb*)imageScanline+x);
540 
541  //compute index where float value is stored in buffer
542  int index = x + y*image.width();
543 
544  //convert blured value to 0-255 range
545  int blurredColor = QMAX( QMIN( ((int) (255*buffer[index])), 255 ), 0 );
546 
547  //blur the entire thing!
548  float alpha;
549  if( edgeImage == NULL)
550  alpha = 1.0f;
551  else
552  {
553  alpha = edgeValue( x, y );
554  if(!blurEdges)
555  alpha = 1.0f - alpha;
556  }
557 
558  //convert and store correct channel in buffer
559  if( channel == 0 )
560  *rgb = qRgb( (int) (alpha*blurredColor + (1-alpha)*qRed(*rgb)),
561  qGreen(*rgb), qBlue(*rgb) );
562  else if( channel == 1 )
563  *rgb = qRgb( qRed(*rgb),
564  (int) (alpha*blurredColor + (1-alpha)*qGreen(*rgb)),
565  qBlue(*rgb) );
566  else
567  *rgb = qRgb( qRed(*rgb), qGreen(*rgb),
568  (int) (alpha*blurredColor + (1-alpha)*qBlue(*rgb)) );
569  } //x
570  } //y
571 }
572 //==============================================
float * regionColBuffer
Definition: blur.cpp:85
QPoint displayOffset
Definition: blur.cpp:90
void blurColumn(int column)
Definition: blur.cpp:409
float * rowBuffer
Definition: blur.cpp:81
void blurImage(QImage &image, float sigma)
Definition: blur.cpp:94
QImage * edgeImage
Definition: blur.cpp:87
int regionCount
Definition: blur.cpp:89
float B
Definition: blur.cpp:78
float * colBuffer
Definition: blur.cpp:82
float edgeValue(int x, int y)
Definition: blur.cpp:234
void blurRegionsInCol(int x)
Definition: blur.cpp:434
float b3
Definition: blur.cpp:78
float q
Definition: blur.cpp:78
int regionIndex(int x, int y)
Definition: blur.cpp:227
void blurRow(int row)
Definition: blur.cpp:285
float b0
Definition: blur.cpp:78
void computeCoeffs(float sigma)
Definition: blur.cpp:150
int width
Definition: blur.cpp:79
void fillBuffer(QImage &image, int channel)
Definition: blur.cpp:172
void resetImageData(QImage &image, int channel, bool blurEdges)
Definition: blur.cpp:527
float * buffer
Definition: blur.cpp:80
void blurBuffer()
Definition: blur.cpp:204
int * regionMap
Definition: blur.cpp:88
float b1
Definition: blur.cpp:78
float b2
Definition: blur.cpp:78
void blurRegionsInRow(int y)
Definition: blur.cpp:311
QSize fullRes
Definition: blur.cpp:91
float * regionRowBuffer
Definition: blur.cpp:84
int height
Definition: blur.cpp:79