C# Performance

24 03 2011

I have been working C# for nearly 8 years and never lost sleep over performance. Sure there are times when things run slowly, but the bottleneck is usually a database query or a web service call. I have never really needed to look into C# performance because it was always fast enough.

A couple weeks ago, Danny Tuppeny wrote “Why I’m Close to Giving Up on Windows Phone 7, as a User and a Developer”.  He raised some pretty valid points, but of course this was a call to arms for the “Anything But Microsoft” crowd.  A side thread about C# performance caught my eye.  Here are a couple of the comments:

  • “My favorite feature of .Net, in general, is sluggish performance, but C# is the best language to write sluggish software in by far.”
  • “ObjC performance is also usually much better than even unsafe C# code, since it is essentially C.”

As I said before, I have never had reason to consider C# “sluggish”, but I also hadn’t really looked into its performance either.  Let’s cut to the chase.

The Results

These are the results (in seconds) for a naïve nth prime finder.  11, 101, 1001, 10001, and 1000001 are the “n’s” in nth.  31, 547, 7927, 104743, and 1299721 are the actual values for the nth prime.

  Version Trial 11 101 1001 10001 100001
31 547 7927 104743 1299721
C++ 1 1 0.000 0.000 0.017 2.101 263.429
2 0.000 0.000 0.021 2.100 263.472
3 0.000 0.000 0.016 2.103 264.520
4 0.000 0.000 0.016 2.109 263.482
5 0.000 0.000 0.015 2.098 265.413
Average 0.000 0.000 0.017 2.102 264.063
C# 1 1 0.001 0.001 0.017 2.172 264.570
2 0.001 0.001 0.017 2.133 264.378
3 0.001 0.001 0.017 2.114 264.246
4 0.001 0.001 0.016 2.133 264.401
5 0.001 0.001 0.017 2.128 264.654
Average 0.001 0.001 0.017 2.136 264.450

The difference is negligible.

The Code

I chose the nth prime problem for a couple reasons.  The code could be written nearly identical in C++ and C# removing any ambiguity about language semantics.  It also has a wide range of optimizations that can be applied.  Version 1 is an extremely naïve (i.e. not optimized) approach.

C++

#include "stdafx.h"

#include <time.h>

#include <iostream>

 

const int maxPrimeIndex = 100001;

int _tmain(int argc, _TCHAR* argv[])

{

       clock_t start, finish;

       start = clock();

 

       int currentPrime = 2;

       int primeCount = 1;

 

       while (primeCount < maxPrimeIndex)

       {

              bool isPrime = false;

              int candidate = currentPrime + 1;

              while (!isPrime)

              {

                     isPrime = true;

                     for (int factor = 2; factor< candidate; factor++)

                     {

                           if (candidate % factor == 0)

                           {

                                  isPrime = false;

                                  break;

                           }

                     }

                     if (!isPrime)

                     {

                           candidate++;

                     }

              }

              currentPrime = candidate;

              primeCount++;

       }

 

       finish = clock();

       double elapsed = ((double)(finish – start)) / CLOCKS_PER_SEC;

       std::cout << elapsed << std::endl << currentPrime << std::endl;

       return 0;

}

C#

using System;

 

namespace SpeedTest

{

    class Program

    {

        const int maxPrimeIndex = 100001;

        static void Main(string[] args)

        {

            DateTime start, finish;

            start = DateTime.Now;

 

            int currentPrime = 2;

            int primeCount = 1;

 

            while (primeCount < maxPrimeIndex)

            {

                bool isPrime = false;

                int candidate = currentPrime + 1;

                while (!isPrime)

                {

                    isPrime = true;

                    for (int factor = 2; factor < candidate; factor++)

                    {

                        if (candidate % factor == 0)

                        {

                            isPrime = false;

                            break;

                        }

                    }

                    if (!isPrime)

                    {

                        candidate++;

                    }

                }

                currentPrime = candidate;

                primeCount++;

            }

 

            finish = DateTime.Now;

            TimeSpan elapsed = finish – start;

            Console.WriteLine(elapsed.ToString());

            Console.WriteLine(currentPrime);

        }

    }

}

Why?

Why are the results so similar?  Because both C++ and C# are compiled.  When comparing these programs we are really comparing their compilers.  This program is pretty straight-forward (i.e. no function calls, memory allocation, array checking, etc.) so the compiler optimizations are probably pretty similar.  I used the Microsoft C++ compiler, but I also tested it with g++ and there was negligible difference. 

If you really need to see the C++ compiler “beat” the C# compiler change candidate from an int to a long.  The C++ code will run a little more than twice as fast.  My guess is that the % operator is much less efficient on longs in C#, but that’s just a guess.

C# is typically JIT compiled, but it doesn’t have to be.  You can use a tool called ngen to compile the image before running it.  You should have a really good reason before doing this though because it can cause headaches when managing updates and the results in many cases will not be dramatic.

Optimization

As I said before, this code was intentionally not optimized.  We are going to apply a very simple but very effective optimization by only checking possible factors up to (and including) the square root of the candidate.  Here are the results:

  Version Trial 11 101 1001 10001 100001
31 547 7927 104743 1299721
C++ 2 1 0.000 0.000 0.000 0.016 0.483
2 0.000 0.000 0.000 0.017 0.479
3 0.000 0.000 0.000 0.015 0.472
4 0.000 0.000 0.000 0.015 0.468
5 0.000 0.000 0.000 0.014 0.481
Average 0.000 0.000 0.000 0.015 0.477
C# 2 1 0.001 0.001 0.002 0.016 0.484
2 0.001 0.001 0.001 0.018 0.474
3 0.001 0.001 0.001 0.016 0.470
4 0.001 0.001 0.002 0.018 0.473
5 0.001 0.001 0.002 0.019 0.470
Average 0.001 0.001 0.002 0.017 0.474

With one simple optimization we have drastically increased our efficiency as n increases.  This raises another question: how fast is fast enough?  For the 100001st it was definitely worth our while.  For the 10001st we save a couple of seconds.  For the 1001st and below it took us more time to write this simple optimization than we’ll ever save.  Context is important.

Conclusion

The point is that people that make blanket statements about performance often have no idea what they are talking about.  Good programmers take the time to understand bottlenecks.  Their gut reaction isn’t “throw more hardware at it” or “should have written it in C”.  A great compiler isn’t going to save your application from a bad coder.

If you disagree, let me know.  Maybe you’re curious how JavaScript or Python compares?  I’d be happy to oblige.  You can find me on Twitter (@azzlsoft) or email (rich@azzlsoft.com).





Phone Vision 18 – Erosion and Dilation

23 03 2011

Our exploration of set theory last time was a precursor to discussing dilation and erosion.  Simply put, dilation makes things bigger and erosion makes things smaller.  Seems obvious right?  Let’s see if the mathematicians can muck it up.

Translation

This is not necessarily a fundamental set operation, but it is critical for what we will be doing. Translation can be described as shifting the pixels in the image. We can move it in both the x and y directions and we will represent the amount of movement using a point. To translate the set then, we simply add that point to every element of the set. For instance, using the point (20, 20) the puzzle is moved down 20 pixels and right 20 pixels like this:

image

The light blue area marks the original location.

Dilation

Dilation is an extension of translation from a point to a set of points. If you union the sets created by each translation you end up with dilation.  In the following diagram, the black pixels represent the original puzzle piece.  The light blue halo represents several translations of the original.  Each translation is translucent so I have tried to make the image large enough so you can see the overlap.  Notice how the hole is missing?

image

The set of points that were used to create the translated sets is called the structuring element.  A common structuring element is

{ (-1,-1), (0,-1), (1,-1),(-1,0), (0,0), (1,0),(-1,1), (0,1), (1,1) }

This will expand the image by one pixel in every direction.  If we think about this as a binary image, it might be represented by this:

image

Kinda reminds me of one of our convolution kernels.  Structuring elements are often called “probes”.  This one is centered at (0,0).

Using the PixelSet we created last time, the code is pretty straightforward.

public void Dilate(PixelSet structuringElem)

{

    PixelSet dilatedSet = new PixelSet();

    foreach(int translateY in structuringElem._pixels.Keys)

    {

        foreach(int translateX in structuringElem._pixels[translateY])

        {

            PixelSet translatedSet = new PixelSet();

            foreach (int y in _pixels.Keys)

            {

                foreach(int x in _pixels[y])

                {

                    translatedSet.Add(x + translateX, y + translateY);

                }

            }

            dilatedSet.Union(translatedSet);

        }

    }

    _pixels = new Dictionary<int,List<int>>(dilatedSet._pixels);

}

Here is the puzzle piece dilated by one pixel in all directions.  The hole is gone.

image

Erosion

Erosion is similar to dilation except that our translation is subtraction instead of addition and we are finding the intersection instead of the union.  Here is an image that has been eroded by 10 pixels in all directions.  The light blue area is the original puzzle.

 

image

 

public void Erode(PixelSet structuringElem)

{

    PixelSet erodedSet = new PixelSet(this);

    foreach (int translateY in structuringElem._pixels.Keys)

    {

        foreach (int translateX in structuringElem._pixels[translateY])

        {

            PixelSet translatedSet = new PixelSet();

            foreach (int y in _pixels.Keys)

            {

                foreach (int x in _pixels[y])

                {

                    translatedSet.Add(x – translateX, y – translateY);

                }

            }

            erodedSet.Intersect(translatedSet);

        }

    }

    _pixels = new Dictionary<int,List<int>>(erodedSet._pixels);

}

 

Here is the puzzle eroded by one pixel in all directions.  That hole is a little bigger.

image

Summary

Like I said at the beginning, dilation makes things bigger and erosion makes things smaller.  Interestingly, erosion is simply the dilation of the background and dilation is the erosion of the background.

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2018%20-%20Erosion%20and%20Dilation.zip

Up Next: Opening and Closing





Phone Vision 17–Sets

16 03 2011

Last time we were able to get a pretty decent binary image.  The question is, what do we do with it?  Before we can answer that, we need to make sure you’ve got a basic understanding of set theory.

Sets

In a binary image we can think of the objects of interest as a set of points.  Here is our binary image from last time:

image

I am going to focus on just the top puzzle piece so let’s extract it.  (The code for extracting this is not important for this lesson, but it is included in the code as GetPuzzleImage).

image

We are going to be working with this image for a while.  In this case, the black pixels will represent members of our set.

Universe

The universe is the set of all possibilities.  In the case of image processing this is typically the set of all of the pixels.  I have seen this referred to as E2 (for 2-dimensional Euclidean space).  In our case, the universe is…

image

As expected, the set of every possibility is all black.  This is important some operations rely on this set.

Empty Set

The empty set is, well, empty.

image

Since the boundary for our image is not a specified size, this is just an empty image. If the universe was a specified size then the empty set would often be represented by a completely white image.

Complement

The complement is everything in the universe that is not in the set.  The complement of our puzzle piece:

image

The puzzle piece is now a hole in the black image.

Union (⋃)

The union is the set of all elements in either of two sets.  If you have been playing along, this is very similar to the | (binary OR) operation mentioned in lesson three. Below we union the puzzle set with the set of its complement.  The output is the universe.

image

image

image

Intersection (⋂)

Finally, the intersection is the set of all shared members of two sets.  Below we intersect the puzzle set with its complement.  This is similar to the & (binary AND) operations mentioned way back in lesson 2.  Since the complement of a set  doesn’t share any elements with the original set the output is going to be the empty set.

imageimageimage

Pixel Set

In order to demonstrate these concepts (and a few in the future) I created a class called PixelSet that behaves similar to a mathematical set. 

Caution: This set is not optimized for use in actual applications.  You’ve been warned.

Elements

The elements of the set are represented by a Dictionary.

Dictionary<int, List<int>> _pixels = new Dictionary<int, List<int>>();

The key is the y-value for the pixel and the List<int> represents all x-values in that row.

 

Add

 

public void Add(int x, int y)

{

    List<int> columns = new List<int>();

    if (!_pixels.ContainsKey(y))

    {

        _pixels.Add(y, columns);

    }

    else

    {

        columns = _pixels[y];

    }

 

    if (!columns.Contains(x))

    {

        columns.Add(x);

    }

}

 

The add function is a little tricky because if the row (y-value) doesn’t have a representative in the set yet then we need to create a new list to contain them. If it does, then we need to use the existing list.

 

Remove

 

public void Remove(int x, int y)

{

    if (_pixels.ContainsKey(y))

    {

        if (_pixels[y].Contains(x))

        {

            _pixels[y].Remove(x);

            if (_pixels[y].Count == 0)

            {

                _pixels.Remove(y);

            }

        }

    }

}

 

Remove is included for completeness, but I’m not calling it from any location in the code.

 

Union

 

public void Union(PixelSet otherSet)

{

    foreach (int key in otherSet._pixels.Keys)

    {

        foreach (int col in otherSet._pixels[key])

        {

            Add(col, key);

        }

    }

}

 

A possible optimization here is to loop through the smallest set.  You would need to maintain a count, but since new members are added through Add that should be easy.

 

Intersect

 

public void Intersect(PixelSet otherSet)

{

    PixelSet intersection = new PixelSet();

    foreach (int key in otherSet._pixels.Keys)

    {

        if (_pixels.ContainsKey(key))

        {

            foreach (int col in otherSet._pixels[key])

            {

                if (_pixels[key].Contains(col))

                {

                    intersection.Add(col, key);

                }

            }

        }

    }

    _pixels = new Dictionary<int, List<int>>(intersection._pixels);

}

 

Once again you could change this to loop through the smallest set to speed this up.

 

Summary
 

I expect this will have been a review for many readers, but it’s important to have a grasp of this information before we head into the next lesson regarding erosion and dilation.  I’d also like to reiterate that this is not optimized for image processing.  My HD7 takes almost 8 times longer to execute some of these tight loops.

 

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2017%20-%20Set%20Operations.zip

Up Next: Erosion and Dilation





Phone Vision 15 – Median Filters

1 03 2011

Filters, filters, filters… Tired of them yet?  Me too.  After this we’ll leave them alone for a while.  I promise.  Until then, we have work to do.

What is a median filter?  If you understand the average filter (or mean filter) then the median filter should give you no troubles at all – assuming you know what the median is.  We’ll start there.

What is the Median?

The median is the middle number of a sample.  To find it, you sort your list and pick the middle item.  Yes, it really is that simple.

Here is a random set of 9 numbers:

image

Sort them:

image

Pick the middle value:

image

Like I said, it’s simple.

Why use median instead of mean?

The median is particularly effective against certain types of noise – specifically “salt and pepper” noise (black and white specks).  Let’s assume this is the neighborhood we’re dealing with:

image

Performing the mean gives:

image

image

While the median returns:

image

image

If we expect that the 0 values are noise then the median is probably a much closer estimate.

Code

private WriteableBitmap MedianFilter(WriteableBitmap grayscale, int radius)

{

    // we are still going to create a new image

    // because we don’t want to modify the

    // old image as we are processing it

    WriteableBitmap filtered =

        new WriteableBitmap(

            grayscale.PixelWidth,

            grayscale.PixelHeight);

 

    // boiler plate code for our

    // histogram stuff

    int[] histogram = new int[256];

    int maxIntensity = 0;

 

    // the math is still easier if we create two loops

    // instead of one

    for (int y = 0; y < grayscale.PixelHeight; y++)

    {

        for (int x = 0; x < grayscale.PixelWidth; x++)

        {

            //here’s the pixel we’re centered on

            int pixel = x + y * grayscale.PixelWidth;

            byte intensity = (byte)grayscale.Pixels[pixel];

 

            // if we are on an edge we are going to leave it

            // as the original intensity.  you will see the

            // edges increasingly unsmoothed as the window

            // size increases.  here we are using the radius

            // to determine our bounds

            if (y <= radius – 1 ||

                x <= radius – 1 ||

                y >= grayscale.PixelHeight – radius ||

                x >= grayscale.PixelWidth – radius)

            {

 

                histogram[intensity]++;

                if (histogram[intensity] > maxIntensity)

                {

                    maxIntensity = histogram[intensity];

                }

                continue;

            }

 

            /////////////////////////////////////////////////////////

            // IMPORTANT PART ///////////////////////////////////////

            /////////////////////////////////////////////////////////

            // this list is the key

            // it contains all of the neighboring pixels

            List<byte> localIntensities = new List<byte>();

            for (int yoffset = -radius; yoffset <= radius; yoffset++)

            {

                for (int xoffset = -radius;

                    xoffset <= radius;

                    xoffset++)

                {

                    localIntensities.Add((

                        (byte)grayscale.Pixels[(x + xoffset)

                        + (y + yoffset) * grayscale.PixelWidth]));

                }

            }

            //sort the intensities

            localIntensities.Sort();

            //pick the middle value

            int medianLocalIntensity =

             localIntensities[(int)(localIntensities.Count/2.0+.5)];

            /////////////////////////////////////////////////////////

            // END IMPORTANT PART ///////////////////////////////////

            /////////////////////////////////////////////////////////

                   

            // and now just set the color

            filtered.Pixels[pixel] = (255 << 24)

                | (byte)medianLocalIntensity << 16

                | (byte)medianLocalIntensity << 8

                | (byte)medianLocalIntensity;

 

            histogram[(byte)medianLocalIntensity]++;

            if (histogram[(byte)medianLocalIntensity] > maxIntensity)

            {

                maxIntensity = histogram[(byte)medianLocalIntensity];

            }

        }

    }

 

    PlotHistogram(histogram, maxIntensity);

    return filtered;

}

 

Results

Taking a simple gray image I added salt and pepper noise to the image then performed a median filter and an average filter to it.  The results are stunning.

imageimage

                flat gray image                              10% salt and pepper noise

imageimage

           after median filtering                      after mean filtering

There are a few specks after applying the median filter, but the noise is removed pretty well.  The average filter performed dismally to say the list.

Summary

If you expect salt and pepper noise then the median filter is great tool to have in your toolbox.  If you want you can explore max, min, and mode filters.  I don’t think we’ll cover them here unless we have a specific application for it.

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2015%20-%20Median%20Filter.zip

(code includes salt and pepper noise generation)

Up Next: Binary Images

Feel free to leave a comment or follow @azzlsoft on Twitter if you have any questions.





Phone Vision 14 – Sobel Operators

25 02 2011

Last time we explored how we can use the rate of change between adjacent pixels to enhance a boundary between two features in an image.  Today we are going to take that one step further and isolate edges using the Sobel operators.  Hold on tight…

The Sobel Operators

A search for the Sobel operators will quickly turn up these two filters:

image

image

 

 

 

 

 

 

 

I created a rather colorful graphic to demonstrate some of the concepts.

image

 

imageimageimage

What do you notice?  Here’s what I notice:

  • All of the non-edges turn black. 
  • Diagonal edges show up in both filters. 
  • The first filter misses horizontal edges
  • The second filter misses vertical edges.

The union of these two images will give you a more complete view of the edges in the system.  If that’s all you wanted to know then you’re done and you can stop here.

I’m curious how they work so I’m going to dig deeper.

Gradient

When you are driving down a mountain you might see a sign (like the one to the right) indicating a steep grade.  The grade in the context of a mathematical function is no different.  We are trying to find the slope of the function at a given point.  This is called the gradient.

Last time we looked at how fast a function was changing by subtracting adjacent intensities.    If you think about it, it’s easy to imagine areas with a lot of change having a steep slope, right?  So have we already found the gradient?  Yes and no.  Hold on, this is going to get rough.

Subtracting the adjacent intensities, we determined the slope from one pixel to the next, but we didn’t determine the slope at the pixel.  Because we are dealing with discrete values (i.e. there is no ‘in-between’ pixels) the best we can hope for is an approximation of the slope at a pixel.  I am hoping that an example will clear this up a bit. 

Imagine the sin function represents some hilly terrain you are driving over.

image

It should be fairly apparent that at the top of the hill our grade is flat (slope = 0).  Let’s imagine we have a GPS recording our position every so often so we get an approximation of the terrain.

image

Using the same technique as the last lesson, let’s try to approximate the slope at the top of the hill.  Last time we said the change was f(x+1) – f(x).

image

Well, that’s not very good.  It turns out that if we adjust our formula ever so slightly we get much better results.  What if we use the position right before we got to the top of the hill instead of the position at the top?

image

That is a much better approximation,but what now?

Gradient Edge Detection

Let’s convert our new difference function to our filter representation. 

Remember,

  • f(x+1) = next position (right after the top)
  • f(x) = current position (at the top)
  • f(x-1) = previous position (right before the top)

so our gradient function looks like:

-1*f(x-1) + 0* f(x) + 1*f(x+1)

or

image

Following the same process in the y direction we end up with

-1*f(y-1) + 0* f(y) +1*f(y+1)

or

image

If we apply these we end up with:

imageimageimage

It’s a much weaker signal, but I can still spot the edges. 

If these filters make sense to you then you should be able to answer the following questions:

  • Why are the horizontal lines still missing in the first image?
  • And the vertical lines still missing in the second image? 
  • Why do the diagonals show up?

Separable Filters

You might be thinking, “um, I guess this is nice, but what does it have to do with the Sobel?”  Everything.

We don’t have the tools yet to completely analyze this, but the Sobel operator is a separable filter.  This means that it is a combination of two 1D filters.  Sobel is a smooth filter multiplied with the gradient filter we just built.

A couple of notes:

  • [1,2,1] is a 1D weighted smoothing filter
  • This is matrix multiplication.  If you need a refresher, check out the Kahn Academy video on the subject.  I love that place.

 

image

and

image

Whoa!  Mind blown.

Summary

We are starting to add some useful tools to our toolbox, but we are dancing around some really complex concepts.  Over the next few lessons we will continue working with a few more simple blur, sharpen, and edge detection routines before we jump head first into frequency analysis.

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2014%20-%20Sobel%20Operator.zip

Up Next: Median Filter

Feel free to leave a comment or follow @azzlsoft on Twitter if you have any questions.





Phone Vision – 13 Sharpening Filters

21 02 2011

Smoothing an image was pretty intuitive.  When we are smoothing an image we are actually trying to suppress details in the image.  When we sharpen an image we are trying to do the opposite.  So, how do we do that?

Rates of Change

When we are sharpening an image we are trying to enhance the differences.  What are the differences?  Well, let’s take a look at a 1D strip of intensities:

image

If we take a cursory look at these pixels we’ll see that the intensities gradually decrease from left to right before jumping back up dramatically then dropping back down.  The dramatic ‘jumps’ are what will ultimately be enhanced when we sharpen the image.  Here they are highlighted.

image

We want to know how fast the intensities are changing.  If ∆ represents the rate of change between adjacent pixels then our rate of change looks like:

image

We just subtract the value of the next pixel from the value of the current pixel.  A little more formally, this looks like

∆(x) = f(x+1) – f(x)

Under normal conditions, some natural variation will occur pixel by pixel.  We want to ignore that natural variation and only concern ourselves when the change is drastic.  For that we will have to perform the subtraction one more time.

image

Now we can see how fast each intensity is changing.  Once again, let’s formalize this.

2(x) =∆(x+1) – ∆(x)

2(x) = (f(x+1) – f(x)) – (f(x) – f(x-1)

2(x) = f(x-1) – 2 * f(x) + f(x+1)

Implementing the 1D Case

Then the filter for our 1D case  looks like

image

If we apply this to our image we end up with this from above:

image

If we subtract this from our original image something magical happens:

image

We have sharpened the image!

In order to achieve that result, we subtracted ∆2  from the original image f.  Our new formula looks like:

sharp = f – ∆2

so in the x direction,

sharp(x) = f(x) – ∆2(x)

sharp(x) = f(x) – (f(x-1) – 2 * f(x) + f(x+1))

sharp(x) = f(x) + -f(x-1) + 2 * f(x) – f(x+1)

sharp(x) = –f(x-1) + 3*f(x) – f(x+1)

Hence, our final 1D filter:

image

Expanding to 2D

If all images were just a single strip of pixels we’d be done.  Unfortunately they aren’t, but accounting this is a pretty simple adjustment.  The first thing to notice is that the 1D case works in any direction.  If we wanted to apply the filter vertically instead of horizontally it would look like this:

image

So, how do we expand this to two dimensions?  It’s pretty easy, actually:

image

Now you have a sharpening filter that will work in two dimensions.  If you want to include the diagonals, there is no reason you can’t.  Using the exact same method the final sharpening filter becomes:

image

And we’re done!

Summary

With some luck, you will have a solid grasp of how the sharpening filter works.  If you are interested in researching this topic more, we walked through creating the Laplacian operator.  If you are comfortable with calculus, the difference functions we discussed were actually just the first derivative (∆) and the second derivative (∆2).

Having discussed the filter code at length, the code is not posted here.  The code from previous lessons will work fantastically.

Up Next: Sobel Operators





Phone Vision 12–Smoothing Filters

15 02 2011

What does it mean to smooth an image?  When I imagine a smooth surface it doesn’t have a lot of bumps or pits.  If we were to try to smooth out a rough piece of wood we would probably sand down the peaks and fill in holes.  If you apply that analogy to an image it sounds a lot like averaging a neighborhood.

Average Filter

If you remember from last time we described a neighborhood filter using a number of weights as follows:

image

A simple average would allow all of the pixels in the neighborhood to contribute equally to the average.  In that case would represent the filter as:

image

Using our code from last time, the filter looks like

int[,] filter = new int[,]

{

    { 1, 1, 1 },

    { 1, 1, 1 },

    { 1, 1, 1 }

};

And here are the results when we run it:

image

image

The original is on the left and the smoothed version is on the right.  Notice how the large objects are relatively unchanged, but the small objects have lost a lot of information.  In addition to removing noise, if we are trying to segment large objects, we can use smoothing to blur out the small objects.

Modified Filter Code

If we run this filter enough times we will eventually get every pixel to be some sort of gray.  It would be more efficient just to increase the size of the window.  Let’s modify the filtering code to accept a square filter of arbitrary size.

private WriteableBitmap Filter(WriteableBitmap grayscale, int[,] filter)

{

    // we are still going to create a new image

    // because we don’t want to modify the

    // old image as we are processing it

    WriteableBitmap filtered =

        new WriteableBitmap(

            grayscale.PixelWidth,

            grayscale.PixelHeight);

 

    // boiler plate code for our

    // histogram stuff

    int[] histogram = new int[256];

    int maxIntensity = 0;

 

    //here’s where the magic starts

    //assume the filter is square

    int length = (int)Math.Sqrt(filter.Length);

    int filterMagnitude = 0;

    for (int x = 0; x < length; x++)

    {

        for (int y = 0; y < length; y++)

        {

            filterMagnitude += filter[x, y];

        }

    }

 

    // this will be used for our loops

    // it’s important that our size is odd

    // so that our current pixel can be the

    // center point.

    int radius = length / 2;

 

    // the math is still easier if we create two loops

    // instead of one

    for (int y = 0; y < grayscale.PixelHeight; y++)

    {

        for (int x = 0; x < grayscale.PixelWidth; x++)

        {

            //here’s the pixel we’re centered on

            int pixel = x + y * grayscale.PixelWidth;

            byte intensity = (byte)grayscale.Pixels[pixel];

 

            // if we are on an edge we are going to leave it

            // as the original intensity.  you will see the

            // edges increasingly unsmoothed as the window

            // size increases.  here we are using the radius

            // to determine our bounds

            if (y <= radius – 1 ||

                x <= radius – 1 ||

                y >= grayscale.PixelHeight – radius ||

                x >= grayscale.PixelWidth – radius)

            {

                //maintain the original / non-smoothed

                filtered.Pixels[pixel] = (255 << 24)

                    | (byte)intensity << 16

                    | (byte)intensity << 8

                    | (byte)intensity;

 

                histogram[intensity]++;

                if (histogram[intensity] > maxIntensity)

                {

                    maxIntensity = histogram[intensity];

                }

                continue;

            }

 

            int newIntensity = 0;

            //going from -radius to radius makes the

     //math easier here too

            for (int yoffset = -radius; yoffset <= radius; yoffset++)

            {

                for (int xoffset = -radius;

                    xoffset <= radius;

                    xoffset++)

                {

                    // we loop through each pixel in the neighborhood

                    // and apply the filter. by ‘apply the filter’

                    // I mean multiply it by the appropriate weight

                    newIntensity +=

                        ((byte)grayscale.Pixels

                            [(x + xoffset)

                            + (y + yoffset) * grayscale.PixelWidth])

                        * filter[(yoffset + radius),

                        (xoffset + radius)];

                }

            }

 

            // scale the new intensity back

            newIntensity /= filterMagnitude;

            newIntensity =

                Math.Max(Math.Min(newIntensity, 255), 0);

 

            // and now just set the color

            filtered.Pixels[pixel] = (255 << 24)

                | (byte)newIntensity << 16

                | (byte)newIntensity << 8

                | (byte)newIntensity;

 

            histogram[(byte)newIntensity]++;

            if (histogram[(byte)newIntensity] > maxIntensity)

            {

                maxIntensity = histogram[(byte)newIntensity];

            }

        }

    }

 

    PlotHistogram(histogram, maxIntensity);

    return filtered;

}

 

Now we can really smooth the image.

A Large Filter

A 21×21 filter is about as large as I can do without having the mask wrap to the next line.

int[,] filter = new int[,]

{

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

    { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },

};

 

image

We coded it so that the edges are unfiltered.  To work around this issue you can perform your average with the pixels you have available.  Notice here that the large objects are still relatively unchanged, while the smaller objects have lost a substantial amount of information.

Summary

This technique is pretty easy.  More important however, is that the masks can be changed easily and the effects can be seen instantly.  Once again I recommend that you play around with the masks on your own.  As we’ll see soon they can be used for things far more powerful than blurring an image.

Download Code

http://cid-88e82fb27d609ced.office.live.com/embedicon.aspx/Blog%20Files/PhoneVision/PhoneVision%2012%20-%20Smoothing%20Filters.zip

Up next: Sharpening Filters