-
Notifications
You must be signed in to change notification settings - Fork 0
/
filtering.cpp
448 lines (360 loc) · 14.6 KB
/
filtering.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
/* -----------------------------------------------------------------
* File: filtering.cpp
* Created: 2015-09-22
* -----------------------------------------------------------------
*
* Image convolution and filtering
*
* ---------------------------------------------------------------*/
#include "filtering.h"
#include <cmath>
#include <cassert>
using namespace std;
Image boxBlur(const Image &im, int k, bool clamp) {
// --------- HANDOUT PS02 ------------------------------
// Convolve an image with a box filter of size k by k
// It is safe to asssume k is odd.
// return im; // change this
// --------- SOLUTION PS02 ------------------------------
// Create a new empty image
Image filtered(im.width(), im.height(), im.channels());
int sideSize = int((k-1.0f)/2.0f);
float normalizer = 1.0f/float(k*k);
float accum = 0.0f;
// For every pixel in the image
for (int z = 0; z < filtered.channels(); z++) {
for (int y = 0; y < filtered.height(); y++) {
for (int x = 0; x < filtered.width(); x++) {
// Accumulate the sum in the pixel's kxk neighborhood
accum = 0.0f;
for (int yBox = -sideSize; yBox < -sideSize + k; yBox++) {
for (int xBox = -sideSize; xBox < -sideSize + k; xBox++) {
accum += im.smartAccessor(x-xBox,y-yBox,z,clamp);
}
}
// Assign the output pixel the value from convolution (normalized)
filtered(x,y,z) = accum * normalizer;
}
}
}
return filtered;
}
Image Filter::convolve(const Image &im, bool clamp){
// --------- HANDOUT PS02 ------------------------------
// Write a convolution function for the filter class
// return im; // change this
// --------- SOLUTION PS02 ------------------------------
Image imFilter(im.width(), im.height(), im.channels());
int sideW = int((width-1.0)/2.0);
int sideH = int((height-1.0)/2.0);
float accum;
// For every pixel in the image
for (int z = 0; z < imFilter.channels(); z++) {
for (int y = 0; y < imFilter.height(); y++) {
for (int x = 0; x < imFilter.width(); x++) {
accum = 0.0;
for (int yFilter=0; yFilter<height; yFilter++) {
for (int xFilter=0; xFilter<width; xFilter++) {
// Sum the image pixel values weighted by the filter
// flipped kernel, xFilter, yFilter have different signs in filter
// and im
accum += operator()(xFilter, yFilter) *
im.smartAccessor(x-xFilter+sideW,y-yFilter+sideH,z,clamp);
}
}
// Assign the pixel the value from convolution
imFilter(x,y,z) = accum;
}
}
}
return imFilter;
}
Image boxBlur_filterClass(const Image &im, int k, bool clamp) {
// --------- HANDOUT PS02 ------------------------------
// Reimplement the box filter using the filter class.
// check that your results match those in the previous function "boxBlur"
// return im; // change this
// --------- SOLUTION PS02 ------------------------------
vector<float> fData (k*k, 1.0/(k*k) );
Filter boxFilter(fData, k, k);
Image imFilter = boxFilter.convolve(im, clamp);
return imFilter;
}
Image gradientMagnitude(const Image &im, bool clamp){
// --------- HANDOUT PS02 ------------------------------
// Uses a Sobel kernel to compute the horizontal and vertical
// components of the gradient of an image and returns the gradient magnitude.
//return im; // change this
// --------- SOLUTION PS02 ------------------------------
// Sobel filtering in x direction
Filter sobelX(vector<float>{
-1.0, 0.0, 1.0,
-2.0, 0.0, 2.0,
-1.0, 0.0, 1.0}, 3, 3);
Image imSobelX = sobelX.convolve(im, clamp);
// Sobel filtering in y direction
Filter sobelY(vector<float>{
-1.0, -2.0, -1.0,
0.0, 0.0, 0.0,
1.0, 2.0, 1.0}, 3, 3);
Image imSobelY = sobelY.convolve(im, clamp);
// squared magnitude
Image magnitude = imSobelX*imSobelX + imSobelY*imSobelY;
// take the square root
for(int i=0; i<magnitude.number_of_elements(); i++){
magnitude(i) = sqrt(magnitude(i));
}
return magnitude;
}
vector<float> gauss1DFilterValues(float sigma, float truncate){
// --------- HANDOUT PS02 ------------------------------
// Create a vector containing the normalized values in a 1D Gaussian filter
// Truncate the gaussian at truncate*sigma.
// return vector<float>();
// --------- SOLUTION PS02 ------------------------------
// compute the support of the filter
int offset = int(ceil(truncate * sigma));
int filterSize = 2*offset+1;
vector <float> fData (filterSize, 0.0f);
// compute the un-normalized value of the gaussian
float normalizer = 0.0f;
for( int i = 0; i < filterSize; i++){
fData[i] = exp( -(i - offset)*(i - offset) / (2.0f *sigma*sigma) );
normalizer += fData[i];
}
// normalize
for( int i=0; i < filterSize; i++){
fData[i] /= normalizer;
}
return fData;
}
Image gaussianBlur_horizontal(const Image &im, float sigma, float truncate, bool clamp) {
// --------- HANDOUT PS02 ------------------------------
// Gaussian blur across the rows of an image
// return im;
// --------- SOLUTION PS02 ------------------------------
// Filter in the x direction
vector<float> fData = gauss1DFilterValues(sigma, truncate);
Filter gaussX(fData, fData.size(), 1);
Image imFilter = gaussX.convolve(im, clamp);
return imFilter;
}
vector<float> gauss2DFilterValues(float sigma, float truncate) {
// --------- HANDOUT PS02 ------------------------------
// Create a vector containing the normalized values in a 2D Gaussian
// filter. Truncate the gaussian at truncate*sigma.
// return vector<float>();
// --------- SOLUTION PS02 ------------------------------
// Compute the filter size
int offset = int(ceil(truncate * sigma));
int k = 2*offset+1;
int filterSize = k*k;
vector <float> fData(filterSize, 0);
int idx = 0;
float normalizer = 0.0;
// compute the unnormalized value of the gaussian and put it in a row-major vector
for( int i=0; i<k; i++){
for( int j=0; j<k; j++){
fData[idx] = exp( -( (i-offset)*(i-offset) + (j-offset)*(j-offset) ) / (2.0 * sigma*sigma) );
normalizer += fData[idx];
idx++;
}
}
// normalize
for(int i=0; i < filterSize; i++){
fData[i] /= normalizer;
}
return fData;
}
Image gaussianBlur_2D(const Image &im, float sigma, float truncate, bool clamp) {
// --------- HANDOUT PS02 ------------------------------
// Blur an image with a full full 2D rotationally symmetric Gaussian kernel
// return im;
// --------- SOLUTION PS02 ------------------------------
// Blur using a 2D gaussian filter
vector<float> fData = gauss2DFilterValues(sigma, truncate);
int k = sqrt(fData.size());
Filter gauss(fData, k, k);
Image imFilter = gauss.convolve(im, clamp);
return imFilter;
}
Image gaussianBlur_separable(const Image &im, float sigma, float truncate, bool clamp) {
// --------- HANDOUT PS02 ------------------------------
// Use principles of seperabiltity to blur an image using 2 1D Gaussian Filters
// return im;
// --------- SOLUTION PS02 ------------------------------
// Blur using two 1D filters in the x and y directions
vector<float> fData = gauss1DFilterValues(sigma, truncate);
Filter gaussX(fData, fData.size(), 1);
Filter gaussY(fData, 1, fData.size());
Image imFilter = gaussX.convolve(im, clamp);
imFilter = gaussY.convolve(imFilter, clamp);
return imFilter;
}
Image unsharpMask(const Image &im,
float sigma,
float truncate,
float strength,
bool clamp){
// --------- HANDOUT PS02 ------------------------------
// Sharpen an image
// return im;
// --------- SOLUTION PS02 ------------------------------
// Get the low pass image
Image lowPass = gaussianBlur_separable(im, sigma, truncate, clamp);
// Subtract it from the original image to get the high pass image
Image highPass = im - lowPass;
// Increase the highPass component to sharpen
Image sharp = im + strength * highPass;
return sharp;
}
Image bilateral(const Image &im,
float sigmaRange,
float sigmaDomain,
float truncateDomain,
bool clamp){
// --------- HANDOUT PS02 ------------------------------
// Denoise an image using the bilateral filter
// return im;
// --------- SOLUTION PS02 ------------------------------
Image imFilter(im.width(), im.height(), im.channels());
// calculate the filter size
int offset = int(ceil(truncateDomain * sigmaDomain));
int sizeFilt = 2*offset + 1;
float accum,
tmp,
range_dist,
normalizer,
factorDomain,
factorRange;
// for every pixel in the image
for (int z=0; z<imFilter.channels(); z++)
for (int y=0; y<imFilter.height(); y++)
for (int x=0; x<imFilter.width(); x++)
{
// initilize normalizer and sum value to 0 for every pixel location
normalizer = 0.0f;
accum = 0.0f;
// sum over the filter's support
for (int yFilter=0; yFilter<sizeFilt; yFilter++)
for (int xFilter=0; xFilter<sizeFilt; xFilter++)
{
// calculate the distance between the 2 pixels (in range)
range_dist = 0.0f; // |R-R1|^2 + |G-G1|^2 + |B-B1|^2
for (int z1 = 0; z1 < imFilter.channels(); z1++) {
tmp = im.smartAccessor(x,y,z1,clamp); // center pixel
tmp -= im.smartAccessor(x+xFilter-offset,y+yFilter-offset,z1,clamp); // neighbor
tmp *= tmp; // square
range_dist += tmp;
}
// calculate the exponential weight from the domain and range
factorDomain = exp( - ((xFilter-offset)*(xFilter-offset) + (yFilter-offset)*(yFilter-offset) )/ (2.0 * sigmaDomain*sigmaDomain ) );
factorRange = exp( - range_dist / (2.0 * sigmaRange*sigmaRange) );
normalizer += factorDomain * factorRange;
accum += factorDomain * factorRange * im.smartAccessor(x+xFilter-offset,y+yFilter-offset,z,clamp);
}
// set pixel in filtered image to weighted sum of values in the filter region
imFilter(x,y,z) = accum/normalizer;
}
return imFilter;
}
Image bilaYUV(const Image &im, float sigmaRange, float sigmaY, float sigmaUV, float truncateDomain, bool clamp){
// --------- HANDOUT PS02 ------------------------------
// 6.865 only
// Bilateral Filter an image seperatly for
// the Y and UV components of an image
// return im;
// --------- SOLUTION PS02 ------------------------------
//convert from RGB to YUV
Image imYUV = rgb2yuv(im);
// We pass the whole imYUV to bilateral, since we want to compute the
// weight on the full YUV range
Image bilY = bilateral(imYUV, sigmaRange, sigmaY, truncateDomain, clamp);
Image bilUV = bilateral(imYUV, sigmaRange, sigmaUV, truncateDomain, clamp);
// put the Y and UV parts of the image back into one image
for(int i=0; i<im.width(); i++) {
for(int j=0; j<im.height(); j++) {
imYUV(i,j,0) = bilY(i,j,0);
imYUV(i,j,1) = bilUV(i,j,1);
imYUV(i,j,2) = bilUV(i,j,2);
}
}
// convert from YUV back to RGB
Image bilRGB = yuv2rgb(imYUV);
return bilRGB;
}
/**************************************************************
// DON'T EDIT BELOW THIS LINE //
*************************************************************/
// Create an image of 0's with a value of 1 in the middle. This function
// can be used to test that you have properly set the kernel values in your
// Filter object. Make sure to set k to be larger than the size of your kernel
Image impulseImg(int k){
// initlize a kxkx1 image of all 0's
Image impulse(k, k, 1);
// set the center pixel to have intensity 1
int center = floor(k/2);
impulse(center,center,0) = 1.0f;
return impulse;
}
// ------------- FILTER CLASS -----------------------
Filter::Filter(const vector<float> &fData, int fWidth, int fHeight)
: kernel(fData), width(fWidth), height(fHeight)
{
assert(fWidth*fHeight == (int) fData.size());
}
Filter::Filter(int fWidth, int fHeight)
: kernel(std::vector<float>(fWidth*fHeight,0)), width(fWidth), height(fHeight)
{}
Filter::~Filter()
{}
const float & Filter::operator()(int x, int y) const {
if (x < 0 || x >= width)
throw OutOfBoundsException();
if ( y < 0 || y >= height)
throw OutOfBoundsException();
return kernel[x + y*width];
}
float & Filter::operator()(int x, int y) {
if (x < 0 || x >= width)
throw OutOfBoundsException();
if ( y < 0 || y >= height)
throw OutOfBoundsException();
return kernel[x +y*width];
}
// --------- END FILTER CLASS -----------------------
// --------- HANDOUT PS07 ------------------------------
Image gradientX(const Image &im, bool clamp){
Filter sobelX(3, 3);
sobelX(0,0) = -1.0; sobelX(1,0) = 0.0; sobelX(2,0) = 1.0;
sobelX(0,1) = -2.0; sobelX(1,1) = 0.0; sobelX(2,1) = 2.0;
sobelX(0,2) = -1.0; sobelX(1,2) = 0.0; sobelX(2,2) = 1.0;
Image imSobelX = sobelX.convolve(im, clamp);
return imSobelX;
}
Image gradientY(const Image &im, bool clamp) {
// sobel filtering in y direction
Filter sobelY(3, 3);
sobelY(0,0) = -1.0; sobelY(1,0) = -2.0; sobelY(2,0) = -1.0;
sobelY(0,1) = 0.0; sobelY(1,1) = 0.0; sobelY(2,1) = 0.0;
sobelY(0,2) = 1.0; sobelY(1,2) = 2.0; sobelY(2,2) = 1.0;
Image imSobelY = sobelY.convolve(im, clamp);
return imSobelY;
}
Image maximum_filter(const Image &im, float maxiDiam) {
float mi = floor((maxiDiam) / 2);
float ma = maxiDiam - mi - 1;
Image mf(im.width(), im.height(), im.channels());
for (int c = 0; c < im.channels(); c++)
for (int i = mi; i < im.width() - ma; i++)
for (int j = mi; j < im.height() - ma; j++)
{
for (int deli = -mi; deli <= ma; deli++)
for (int delj = -mi; delj <= ma; delj++)
{
mf(i, j, c) = max(mf(i, j, c), im(i+deli, j+delj, c));
}
}
return mf;
}
// ------------------------------------------------------