- Email: [email protected]
Case study
15 Case study: Median filtering 15.1 Introduction to the case studies
The case studies in Chapters 15 to 17 discuss the development of three real programs. Each illustrates multiple facets of the design processes, with different emphases. Each of the problems being solved is a subproblem of scientific and engineering R&D, and each represents the sort of task an individual scientist or engineer might face in practice. Because the programs are not long, the requirements for specification, modularization and detailed design are modest compared to large software systems. But the principles and practices are constant, as we will see.
15.2 Introduction to this chapter
In this chapter we will develop an w-point median filter for signal processing. Median filtering is widely applied in signal conditioning (immediately after collecting measurements fi-om a sensor) and later manipulation of signals for enhancement, coding or analysis. The particular application that we'll focus on is in video signal processing, and we'll see how this presents some important constraints on the design. On the other hand, we'll always keep in mind that we'd like to develop software that is as widely usefiil as possible. First, we'll look at the context: why you might want to filter signals with a non-linear operator like a median rather than a conventional filter. Then we turn to the problem: the particular video analysis application of median filters that motivated this case study. Next, the steps of the design process will illustrate: Rapid prototyping Rules Make the first version durabie not functionai, and get it running early. • The search for relevant prior art Rule 11 Exploit existing solutions. • Use of a doubly-linked list data structure • Implementation and testing strategy Rule 14 Test to break. • Finishing: polishing the final version.
Case study: Median filtering 311 Several versions of the code will be presented. Beware: the second one contains a subtle bug that does not get resolved until the statistical testing stage. The final version implements a robust and efficient median filter class that solves the problem.
15.3 Background
The median of a set of numbers is the middle one when they are listed in ascending or descending order. So 3 is the median of {5,4,3,2,1}, of {1,2,3,5,8}, and of {3,1,4,1,5,9,2} (which, like most real input sets, isn't in order yet). Those examples all involved an odd number of numbers; when you have an even number, the median is the average of the two middle items in a sorted hst. So the median of {2, 0, 0, 4} is 1. On the face of it, finding the median of m = 2A: + 1 nimibers means sorting the set, then choosing the ^h, while finding the median ofm = 2k numbers means sorting the set, then averaging the A: — 1 and kth values (assuming the sorted list is indexedfi*om0). There are faster ways to calculate a median than doing afixUsort, but we need not consider them here, because our task is to design a median filter. In filtering we have an input signal or sequence 5 = (xj, X2,..., JC„) of numbersfi*oma limited range, and apply to them a sliding window fiinctionf of iQngth m so that the output j ^ is given by yt =f(Xj\J = i - o,... ,i-\- m-
o)
where/could be any fiinction of the samples in the window and o is an offset. Often filters are centred on the output index, and may even be symmetric or antisymmetric. In these cases, the window length is odd, and, as above, we can write m = 2A: + 1 and yi=f(Xj\j
= i - K ... , « + ^)
Some important filter fimctions are linear, and in this case the filtering operation is known as convolution. But non-linear filters can be used too, and one of the most popular is the median filter: yi = median(xj\j = i — k, ... ,i + k) Of course, median filtering could use an even window length, but in all practical cases, an odd-length window is appropriate. We therefore confine our discussion to this case.
15.4 Why use median filtering?
Signals are often filtered to remove noise. This can be done with a linear filter, such as a simple average of the samples in the window, but that often causes unwanted side effects. For example, suppose you had a signal with values {0, 0, 0, 0, 0, 100, 100, 100, 100, 100} that was corrupted by random noise to become {3, 15, —7, 2, —10, 94, 102, 100,112,91}. You could apply an averaging filter with window size 3 to obtain {?, 3.667, 3.333, - 5 , 30.667, 62, 96.667, 104.667, 101, ?} where the '?'s indicate that the filter overlaps the end of the signal. This certainly smoothes out some of the noise, but it also blurs the strong step that is in the original. If, instead, you use a length 3 median filter, the result is {?, 3,2, -7,2,94,100,102,100, ?} which preserves the step discontinuity much better. As a result median filters are commonly used for signal enhancement. If applied to an image, for example,
312 Software Design for Engineers and Scientists a. median filter will remove graininess and point noise but preserve the sharp edges between objects. Linear filters are fairly easy to implement in software. Averaging is particularly easy and efficient: simply keep a running sum of the m samples in the window; as it is slid to the next position, subtract the sample just leaving the window and add the sample just entering; then divide the sum by the window size. A median filter is more complicated - perhaps an ordered list can be maintained with items removed and added as in heaps (see page 214) - but before considering alternatives let us turn to a particular application and look at its demands.
15.5 The application
One important module in a particular kind of video analysis is a temporal median filter. Each output videofi*ameis the median ofm consecutive input fi-ames. That is, each individual output pixel is the median of the m values the pixel at that position had in m consecutive input frames. In afi*ameofrXc pixels (r for rows, c for columns), there will he rX c independent, simultaneously running median filters. By employing such a temporal median filter, a system can distinguish short-term temporal 'noise' (such as an object whizzing across the foreground of a picture over a couple offi*amesand occluding the scene behind)fi*oma sequence of fast cuts between different scenes. The problem that faces us is how to implement a median filter that can do this accurately and efficiently.
15.6 Approacliing the problem
Several features of the temporal median filter problem are immediately evident: • The obvious decomposition of the total problem is to develop a single median filter and then replicate it many times (equal to the number of pixels in a picture - maybe quarter of a million). Objectbased coding is therefore ideal: we will define an abstract data type that encapsulates median filter fimctionality, then instantiate a 2D array of objects of that type. • There is a simple naive algorithm for the median filter: read the m values in the current window into an array; sort the array; return the middle element. Such an implementation could be prototyped very rapidly. • Median filtering is a well-known operation so presumably there are establishe4 tested and true algorithms. Bearing these three in mind, our best approach is to begin coding a class and a test program with the naive algorithm; use this to evolve our public class interface and identify important design and test issues, then turn to prior art and do the 'real' implementation based on a study of the literature. Why is this approach best? Why is postponing the search for a good algorithm until after the first implementation good policy? Partly the reason is the primacy of first version as fi:amework: we want to concentrate on getting the overall structure right first. But a second and probably more important reason is to do with testing. The key problem faced by this design is algorithmic. The naive method will work, but it won't be good enough because it will
Case study: Median filtering 313 be slow and clunky. Yet a high-efficiency method is likely to be tricky to code correctly. Because we are fairly sure that our naive initial method will have little in common with the final version, it will make an excellent reference version for statistical black box testing. At this stage we can already see that once we have basicfimctionality,we will want to generate many test cases automatically. Our initial naive implementation will provide a set of results for these. We can't regard the naive implementation as the gold standard - it may have bugs but it provides a cross check.
15.7 R a p i d prototyping
We can use the framework generator on page 190 to create medianfilter.h, medianfilter.cpp, testmedianfilter.cpp and a Makefile. Then opening medianfilter.h and medianfilter.cpp in side-by-side windows, we begin by adding expected data items and fimctions to medianfilter.h. / / Header for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / John Robinson #ifndef medianfilter_H #define medianfilter_H class medianfilter { double *buf; int length; int bufptr; / / Index into circular buffer public: medianfilter (const int len); ^medianfilter (); double addnext (const double newval); / / Returns value displaced from buffer double getmedianO const; void print() const;// Debug function to show current state of / / buffer };
#endif As usual the choice of data structures affects the choice of algorithms and vice versa. We have to store the members in the set somewhere, so our first decision is to have a data item buf that will point to an array of length length. But how will items be added and removed? Assuming we want our data type to function as a filter, we should keep the current window in the object and each time we slide the window by one position of the input, we load the new item in and displace the item that's just left the window. So we need to store the window contents in a circular buffer. Therefore we need bufptr, an index into the buf array showing where the next item will be placed. We have also defined the public interface for the median filter: a constructor in which the length of the window is specified, a member fimction addnext() which adds a new value, causing the displacement of the oldest value (which it returns, which may be helpfiil in
314 Software Design for Engineers and Scientists debugging), and a member function getmedian() which simply returns the current median. Finally there is a print() function to show the current state of the buffer, which, again, will be useful in testing. What seems like a big decision has also been made - that the samples will be doubles. Indeed, the type of the samples is an important issue: since we're dealing with video coding, our samples will actually all be integers between 0 and 255, so could be represented as ints, unsigned ints or even unsigned chars. At this stage, we're fine choosing a type (double) and developing with that. Later we can make the code into a template version to handle diverse types, or perhaps have a single typedef declaration that can be changed depending on the application. The next thing we do is think about consts. Already the class declaration includes some consts, where values are passed into functions. But we should also decide at this stage which member functions leave the object unchanged and preface their declaration with const. printQ should just output current state information and should therefore be const. In our naive implementation, we would expect to code the sorting as part of getmedian() but that need not alter the object itself. Therefore, we also provisionally declare getmedianQ as const. We now move on to write our initial implementation file medianfilter.cpp: / / Implementation for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / John Robinson #include "medianfilter.h" #include #include using namespace std; / / Private function for qsort. WiU remove when / / we improve median finding static int comparedoubles(const void *a, const void *b) { const double *da((const double *)a); const double *db((const double *)b); const double diff(*da-*db); if (diff < 0) return-1; else if (diff > 0 ) return 1; return 0; } medianfilter: :medianfilter(const int len) { length = len; buf = new double [length]; for (int i = 0; i < length; i++) buf[i] = 0; bufptr = 0; } medianfilter: :~medianfilter() { }
Case study: Median filtering 315 double medianfilter::addnext(const double newval) { double retval = buf[bufptr]; buf[bufptr] = newval; if (++bufptr == length) bufptr = 0; return retval; } double medianfilter::getmedian() const { / / Initial implementation: do sort each time. Must improve! double tempbuf[length]; for (int i = 0; i < length; i++) tempbuf[i] = buf[i]; qsort(tempbuf, length, sizeof(double), comparedoubles); if (length & 1) / / length is odd return tempbuf[length/2]; return (tempbuf[length/2] + tempbuf[length/2-l])/2; } void medianfilter::print() const { for (int i = 0; i < length; i++) cout « "buf[" « i « "] = " « buf[i] « endl; cout « "Current buffer index is " « bufptr « endl; } The code here is straightforward. Only two things to note. First, our window is a circular buffer, and therefore we have the test on bufptr in addnextO which cycles it back to the start. Second, we have used qsort() as discussed in Chapter 4 to sort a copy of the buffer. (If we were using the standard template library, we would have accomplished the same thing with the function sort() on a vector.) Obviously, as the code itself says, this is a clunky implementation, but it is easy to code quickly, and thus it forms our first attempt. We now develop the test framework for our new class. Initially we want to check that the circular buffer is working correctly, the true medians are being calculated, and the implementation does support arbitrary window sizes. We therefore make our initial test program an interactive loop where we enter values by hand and at each stage check that both the output and the internal state are correct: / / Test program for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 #include "medianfilter.h" #include using namespace std; int main(int argc, char *argv[]) { int filterlengh; int inputlength; double val; cout « "Enter filter length:"; cin » filterlengh; cout « "Enter input length (number of iterations):"; cin »inputlength;
316 Software Design for Engineers and Scientists medianfilter testmedianfilter(filterlengh); cout « "Buffer state:\n"; testmedianfilter.printO; for (int i = 0; i < inputlength; i++) { cout « "Next input value:"; cin » val; val = testmedianfilter.addnext(val); cout « "Value displaced was " « val « endl; val = testmedianfilter.getmedianO; cout « "Current median is " « val « endl; cout « "Buffer state:\n"; testmedianfilter.printO; } return 0; }
The test program uses all the debugging features provided by the class itself to monitor what is going on. When we compile and test the program it seems to work - except for negative filter lengths. Adding numbers one by one, we keep a running hand calculation of the median, and the program agrees. We therefore have an initialfi*ameworkwhich can serve as reference comparison for our final version.
15.8 Exploit existing solutions
Our search for existing solutions begins with a search for 'median filtering' on Citeseer, Yahoo! and Google (see Chapter 11). Variations like 'median filter' + algorithm also are worth trying. Very quickly wefinda pointer to http://iris.usc.eduA^ision-Notes/bibliography/twod266.html which is a whole list of pubUcations on medianfiltering.Digging deeper draws a blank on online resources but some good leads on papers. One that sounds particularly promising is: M. Juhola, J. Katajainen and T. Raita, 'Comparison of algorithms for standard median filtering', IEEE Transactions on Signal Processing, Vol. 39, No. 1, January 1991, pp. 204-208. It's a little old, but most of the more recent stuff seems to be about hardware, and it has two attractive features: (1) it is a comparison paper, so it will be concise in describing algorithms and will give advice about how to choose an appropriate one, (2) it's only five pages long. So we head for the library and get a copy. With permission of the authors and the IEEE, that copy is reproduced in the Appendix, page 398. This paper seems like just what we want. After one read-through, it looks likely that we will want to concentrate on method D or E, since these are clearly the fastest. See Figure 4 on page 401 and Table 1 on page 402. Our next step is therefore to rough out an implementation based on one of these. When I originally attempted the design, I chose method E, because it's fastest overall, and its restriction to integers doesn't matter for the end application of video processing. I began by ensuring I understood the method, sketching out a histogram and running some numbers through the algorithm by hand, 'playing computer'. To be sure, I wrote the code for the key central section of the algorithm that involves stepping a counter up (or down) through the
Case study: Median filtering 317 number of samples that are equal to the current median as samples are added to each side. This was important because the authors' description of the algorithm didn't explicitly say what happens when a new item has the same value as the current median, and I needed to check the assumptions I was making about that. Next I rewrote the header file and implementation just to use integers fi*om lowest to highest. This introduced a lot of extra range checks so I decided to restrict the range to [0,255], since these would be the values in the video fi-ames. And then ... at last... I realized what was wrong with method E. Although the histogram method is fast, each median filter requires a histogram of bins (i.e. an array of counters) in the range of input values. Even restricting that range to [0, 255], and the window size to 255 (so that each histogram bin can be an unsigned char), means that each filter object must maintain 256 bytes of histogram. If we have a quarter million median filters running at once (for a typicalsized picture), that's 64 Mbytes of memory for what is just one component of the system. Sixty-four Mbytes is perhaps not prohibitive, but it is still a huge overhead on what method D requires. Given that D can cope with unrestricted ranges, it surely is the better choice. How long did I waste on method E? About an hour. Could I have avoided that? Certainly - right there on page 207 of the paper (page 402 of this book), the authors state that the space requirement of method E is 0 ( t / + m) where U'ls the number of integers in the range and m is the window size. All their other tested methods required space 0(m). I missed that, or perhaps I didn't register it until I'd started programming. But the cost was not too great, because developing incrementally exposed the problem more quickly than if I'd been occupied trying to write the whole of method E in one go. Method D is a doubly linked list method, and as we know, linked list manipulation can be tricky. We always have to keep in mind what happens at the ends of the lists. So here is a carefully considered implementation that treats the ends specially. First, the header file which includes the whole class definition for a linked list node called bufitem. // // // //
Header for medianfilter class, created Thu Aug 28 12:42:56 2003 Revised Aug 28 15:15 John Robinson
#ifndef medianfilter_H #de£ine medianfilter_H class bufitem { / / All items private. Just used by friends double value; bufitem ^smaller;// Using pointers rather than indices / / means bufitems can be collected into other DSs than just / / arrays, bufitem *greater; bufitem 0 { value = 0; smaller = 0;
318 Software Design for Engineers and Scientists greater = 0; } '-bufitemO { }
void set(const double v) { value = v; } double get() const { return value; } void initpointers(bufitem *const predecessor,bufitem *const successor) { smaller = predecessor; greater = successor; } void putafter(bufitem *predecessor) { / / Break existing link from predecessor to its successor / / and insert this one bufitem *successor = predecessor->greater; predecessor->greater = this; if (successor) successor->smaller = this; smaller = predecessor; greater = successor; }
void putbefore(bufitem *successor) { / / Can't just use putafter(successor->before) because / / this may be 0 on list bufitem *predecessor = successor->smaller; successor->smaller = this; if (predecessor) predecessor->greater = this; smaller = predecessor; greater = successor; } void removefromlistsO { if (smaller) smaller->greater = greater; if (greater) greater->smaller = smaller; } bufitem *next() const { return greater; } bufitem *prev() const { return smaller; } friend class medianfilter; };
class medianfilter { bufitem *buf; int length;
Case study: Median filtering 319 int nextinsertion; / / Index into circular buffer bufitem *medianloc; / / Not index but pointer public: medianfilter (const int len); ~medianfilter (); double addnext(const double newval); / / Returns value displaced from buffer double getmedianO const; void print0 const; / / Debug function to show current //state of buffer };
#endif Here is the implementation file. Note that most of the fimctionality is in addnextO as implied by the paper. Although care has been taken with pointer manipulation, this version does have a bug. See if you can track how the data structures and algorithms on page 401 have been implemented in the code, and perhaps spot the bug. // // // //
Implementation for medianfilter class, created Thu Aug 28 12:42:56 2003 Revised Aug 28 15:15 John Robinson
#include "medianfilter.h" #include #include using namespace std; medianfilter: :medianiilter(const int len) { length = len; buf = new bufitem[length]; / / Have to treat the items at each end differently to the others / / because they are at the ends of the linked fists buf[0].set(0); buf[0].initpointers(0, &buf[l]); for (int i = 1; i < length-1; i++) { buf[i].set(0); buf[i] .initpointers(&buf[i-1] ,&buf[i+1]); } buf[length-l].set(0); buf[length-1] .initpointers(&buf[length-2] ,0); nextinsertion = 0; medianloc = &buf[length/2];// Actually doesn't matter / / because initial buf values are all equal to 0 } medianfilter: :~medianfilter() { }
double medianfilter::addnext(const double newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length)
320 Software Design for Engineers and Scientists nextinsertion = 0; double oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf[thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) {// May need to move this node upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; } if (p->next()->get() >= newval) { buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); break; } p = p->next(); } if(!p->next()){ buf[thisinsertion] .remove£romlists(); buf[thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prevO; while(p->prev()) { if ((p == medicLnloc)&&!movedmedian) { medianloc = medianloc->prev(); movedmedian = 1; } if (p->prev()->getO <= newval) { buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putbefore(p);
Case study: Median filtering 321 break; } p = p->prev(); } if(!p->prev()){ buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); } } return oldval; } double mediaiifilter::getmedian() const { / / Initial implementation: do sort each time. Must improve! return medianloc->get(); } void medianfilter::print() const { for (int i = 0; i < length; i++) { cout « "buf[" « i « "] = " « buf[i].get(); cout « " has address " « (void *) &buf[i]; cout « ",prev pointer" « (void *) buf[i].prev(); cout « ", next pointer " « (void *) buf[i].next() « endl; } cout « "Current buffer index is " « nextinsertion « endl; }
After compiling and quick hand testing, this version seems to work. But here is where we need to do rigorous, extensive testing. Because the code is a little tricky to understand, white box testing can't be guaranteed to uncover situations that have not been properly catered for in the code. Black box testing by hand is also susceptible to missed cases: we could enumerate possibilities like add-value-greater-than-currentmedian-when-current-median-equals-value-to-remove-from-thewindow-but-more-than-one-other-value-still-in-window-equalscurrent-median, but it's going to be tedious to devise test sets that cover all combinations. This is a case when statistical testing is justified: we can generate long streams of random numbers, filter with both our initial naive method and our 'provisionally final' method and look for differences. Generating long random sequences is easy and we can even use the stable data type discussed in Chapter 17 to help analyse outputs. When we find a difference, we can check by hand which of the two systems is wrong, then work step by step using printQ (or other debugging tools) to find the problem. When we do such testing on the above implementation we discover its problem: addnextQ does not properly reset the pointer to the median when the new item becomes the median. This can be fixed as demonstrated in the version of addnext() listed below. The location of the fix is flagged by comments. double medianfilter::addnext(const double newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length)
322 Software Design for Engineers and Scientists nextinsertion = 0; double oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf[thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) { / / May need to move this node / / upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { / / THE NEXT 11 LINES HX THE BUG BY REORDERING //TESTS / / A CORRESPONDING 1 l-LINE FDC IS FLAGGED //LOWER DOWN if (p->next()->get() >= newval) { buf [thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->next(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; } p = p->next(); } if(!p->next()){ buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prev(); while(p->prev()) {
Case study: Median filtering 323 / / THE NEXT 11 LINES FIX THE BUG BY REORDERING //TESTS / / A CORRESPONDING 11-LINE FIX IS FLAGGED //HIGHER UP if (p->prev()->get() <= newval) { buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->prev(); break; } if ((p == medianloc)&&!movedinedian) { medianloc = medianloc->prev(); movedmedian = 1; } p = p->prev(); } if(!p->prev()){ buf [thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); } } return oldval; }
15.9 Finishing
Now extensive random number tests succeed for the implementation. We therefore consider what should be done to finish the program. The helper class bufitem is only used inside medianfilter. It could be an embedded class, but at the very least it can be moved out of the header file and into the implementation. How shall we make our median filter work on other types than doubles (we're particularly interested in ints for video processing)? One possibility is to use templates, but because export is not supported by many compilers, that would mean relocating a large part of the implementation to the header file. Instead we use a typedef in the header file with a comment to remind us that this has to be changed and things recompiled for different types. We also need to fix the problem identified in the very first version: a negative argument to the constructor causes errors. There are a number of ways of dealing with this, but the best option is to signal a problem but proceed anyway with a default value replacing what the user asked for We need a default value for the constructor too, so that arrays of medianfilters can be instantiated. Having made all those changes, here is the final version. First the header file: / / Header for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 //Revised Aug 28 //Version 0.1 / / John Robinson #i&idef medianfilter_H #define medianfilter_H
324 Software Design for Engineers and Scientists typedef double itemtype; / / Change and recompile for other / / input sample types class medianfilter { class bufitem *buf; int length; int nextinsertion; / / Index into circular buffer class bufitem *medianloc; / / Not index but pointer public: medianfilter (const int len =15); -medianfilter (); itemtype addnext(const itemtype newval); / / Returns value displaced from buffer itemtype getmedianQ const; void printQ const; / / Debug function to show current state / / o f buffer };
#endif Note the addition of the class keyword before bufitem's first mention, because now we have taken the declaration of bufitem out of the header. Here is the final implementation file: / / Implementation for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / Revised Aug 28 ... Corrected same day //Version 0.1 / / John Robinson #include "medianfilter.h" #include #include using namespace std; / / bufitem now defined in implementation file because not part / / of interface class bufitem { / / All items private. Just used by friend medianfilter itemtype value; bufitem ^smaller;// Using pointers rather than indices means / / bufitems can be collected into other DSs / / than just arrays, bufitem *greater; bufitem 0 { value = 0; smaller = 0; greater = 0; } ~bufitem() { }
void set(const itemtype v) { value = v;
Case study: Median filtering 325 } itemtype get() const { return value; } void initpointers(bufitem *const predecessor.bufitem *const successor) { smaller = predecessor; gfreater = successor; } void putafter(bufitem *predecessor) { / / Break existing link from predecessor to its successor / / and insert this one bufitem *successor = predecessor->greater; predecessor->greater = this; if (successor) successor->smaller = this; smaller = predecessor; greater = successor; } void putbefore(bufitem *successor) { / / Can't just use putafter(successor->before) because / / this may be 0 on list bufitem ^predecessor = successor->smaller; successor->smaller = this; if (predecessor) predecessor->greater = this; smaller = predecessor; greater = successor; } void removefromlistsO { if (smaller) smaller->greater = greater; if (greater) greater->smaller = smaller; } bufitem *next() const { return greater; } bufitem *prev() const { return smaller; } friend class medianfilter; }; medianfilter: :medianfilter(const int len) { if (len < 0) { / / Can add an upper limit later if necessary cout « "Specified length" « len « " is out of boundsXn"; cout « "Using length = 15 insteadXn"; length = 15; } else
326 Software Design for Engineers and Scientists length = len; buf = new biifitem[length]; / / Have to treat the items at each end differently to the others / / because they are at the ends of the Unked lists buf[0].set(0); buf[0].initpointers(0, &buf[l]); for (int i = 1; i < length-1; i++) { buf[i].set(0); buf [i] .initpointers(&buf [i-1 ] ,&buf[i+1 ]); } buf[length-l].set(0); buf[length-1 ] .initpointers(&buf [length-2] ,0); nextinsertion = 0; medianloc = &buf{length/2]; / / Actually doesn't matter because initial / / buf values are all equal to 0 } medianfilter::~medianfilter() { } itemtype medianfilter::addnext(const itemtype newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length) nextinsertion = 0; itemtype oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf [thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) {// May need to move this node upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { if (p->next()->get() >= newval) { buf[thisinsertion] .remove£romlists(); buf[thisinsertion] .putafter(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->next(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; }
Case study: Median filtering 327 p = p->next(); } if (!p->next()) { buf [thisinsertion] .remove£romlists(); buf [thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prev(); while(p->prev()) { if (p->prev()->get() <= newval) { buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->prev(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->prev(); movedmedian = 1; } p = p->prev(); } if(!p->prev()){ buf[thisinsertion] .removefromlistsQ; buf [thisinsertion] .putbefore(p); } } return oldval; }
itemtype medianfilter::getmedian() const { return medianloc->get(); } void medianfilter::print() const { for (int i = 0; i < length; i++) { cout « "buf[" « i « "] = " « buf[i].get(); cout « " has address " « (void *) &buf[i]; cout « ",prev pointer " « (void *) buf[i].prev(); cout « ", next pointer" « (void *) buf[i].next() « endl; } cout « "Current buffer index is " « nextinsertion « endl; }
328 Software Design for Engineers and Scientists Is this the end of the story? Not quite. addnextQ can be rewritten to avoid repetition. The two cases (newval < oldval and newval > oldval) can be combined by adding extra functions to bufitem that switch on a parameter to say whether the search is going upwards or downwards. Making the change certainly tightens up the code, but the clarity of shortness and saying something only once is traded for more variables being maintained and repeatedly checked. The upshot is an execution time increase of about 25%, so such a change is probably best avoided. Finally we install the module in the system where it will solve the problem. The median filter developed here was integrated in a video processing system, where it achieved an order of magnitude speedup over a previous median filter that sorted the window each time.
The case studies in Chapters 15 to 17 discuss the development of three real programs. Each illustrates multiple facets of the design processes, with different emphases. Each of the problems being solved is a subproblem of scientific and engineering R&D, and each represents the sort of task an individual scientist or engineer might face in practice. Because the programs are not long, the requirements for specification, modularization and detailed design are modest compared to large software systems. But the principles and practices are constant, as we will see.
15.2 Introduction to this chapter
In this chapter we will develop an w-point median filter for signal processing. Median filtering is widely applied in signal conditioning (immediately after collecting measurements fi-om a sensor) and later manipulation of signals for enhancement, coding or analysis. The particular application that we'll focus on is in video signal processing, and we'll see how this presents some important constraints on the design. On the other hand, we'll always keep in mind that we'd like to develop software that is as widely usefiil as possible. First, we'll look at the context: why you might want to filter signals with a non-linear operator like a median rather than a conventional filter. Then we turn to the problem: the particular video analysis application of median filters that motivated this case study. Next, the steps of the design process will illustrate: Rapid prototyping Rules Make the first version durabie not functionai, and get it running early. • The search for relevant prior art Rule 11 Exploit existing solutions. • Use of a doubly-linked list data structure • Implementation and testing strategy Rule 14 Test to break. • Finishing: polishing the final version.
Case study: Median filtering 311 Several versions of the code will be presented. Beware: the second one contains a subtle bug that does not get resolved until the statistical testing stage. The final version implements a robust and efficient median filter class that solves the problem.
15.3 Background
The median of a set of numbers is the middle one when they are listed in ascending or descending order. So 3 is the median of {5,4,3,2,1}, of {1,2,3,5,8}, and of {3,1,4,1,5,9,2} (which, like most real input sets, isn't in order yet). Those examples all involved an odd number of numbers; when you have an even number, the median is the average of the two middle items in a sorted hst. So the median of {2, 0, 0, 4} is 1. On the face of it, finding the median of m = 2A: + 1 nimibers means sorting the set, then choosing the ^h, while finding the median ofm = 2k numbers means sorting the set, then averaging the A: — 1 and kth values (assuming the sorted list is indexedfi*om0). There are faster ways to calculate a median than doing afixUsort, but we need not consider them here, because our task is to design a median filter. In filtering we have an input signal or sequence 5 = (xj, X2,..., JC„) of numbersfi*oma limited range, and apply to them a sliding window fiinctionf of iQngth m so that the output j ^ is given by yt =f(Xj\J = i - o,... ,i-\- m-
o)
where/could be any fiinction of the samples in the window and o is an offset. Often filters are centred on the output index, and may even be symmetric or antisymmetric. In these cases, the window length is odd, and, as above, we can write m = 2A: + 1 and yi=f(Xj\j
= i - K ... , « + ^)
Some important filter fimctions are linear, and in this case the filtering operation is known as convolution. But non-linear filters can be used too, and one of the most popular is the median filter: yi = median(xj\j = i — k, ... ,i + k) Of course, median filtering could use an even window length, but in all practical cases, an odd-length window is appropriate. We therefore confine our discussion to this case.
15.4 Why use median filtering?
Signals are often filtered to remove noise. This can be done with a linear filter, such as a simple average of the samples in the window, but that often causes unwanted side effects. For example, suppose you had a signal with values {0, 0, 0, 0, 0, 100, 100, 100, 100, 100} that was corrupted by random noise to become {3, 15, —7, 2, —10, 94, 102, 100,112,91}. You could apply an averaging filter with window size 3 to obtain {?, 3.667, 3.333, - 5 , 30.667, 62, 96.667, 104.667, 101, ?} where the '?'s indicate that the filter overlaps the end of the signal. This certainly smoothes out some of the noise, but it also blurs the strong step that is in the original. If, instead, you use a length 3 median filter, the result is {?, 3,2, -7,2,94,100,102,100, ?} which preserves the step discontinuity much better. As a result median filters are commonly used for signal enhancement. If applied to an image, for example,
312 Software Design for Engineers and Scientists a. median filter will remove graininess and point noise but preserve the sharp edges between objects. Linear filters are fairly easy to implement in software. Averaging is particularly easy and efficient: simply keep a running sum of the m samples in the window; as it is slid to the next position, subtract the sample just leaving the window and add the sample just entering; then divide the sum by the window size. A median filter is more complicated - perhaps an ordered list can be maintained with items removed and added as in heaps (see page 214) - but before considering alternatives let us turn to a particular application and look at its demands.
15.5 The application
One important module in a particular kind of video analysis is a temporal median filter. Each output videofi*ameis the median ofm consecutive input fi-ames. That is, each individual output pixel is the median of the m values the pixel at that position had in m consecutive input frames. In afi*ameofrXc pixels (r for rows, c for columns), there will he rX c independent, simultaneously running median filters. By employing such a temporal median filter, a system can distinguish short-term temporal 'noise' (such as an object whizzing across the foreground of a picture over a couple offi*amesand occluding the scene behind)fi*oma sequence of fast cuts between different scenes. The problem that faces us is how to implement a median filter that can do this accurately and efficiently.
15.6 Approacliing the problem
Several features of the temporal median filter problem are immediately evident: • The obvious decomposition of the total problem is to develop a single median filter and then replicate it many times (equal to the number of pixels in a picture - maybe quarter of a million). Objectbased coding is therefore ideal: we will define an abstract data type that encapsulates median filter fimctionality, then instantiate a 2D array of objects of that type. • There is a simple naive algorithm for the median filter: read the m values in the current window into an array; sort the array; return the middle element. Such an implementation could be prototyped very rapidly. • Median filtering is a well-known operation so presumably there are establishe4 tested and true algorithms. Bearing these three in mind, our best approach is to begin coding a class and a test program with the naive algorithm; use this to evolve our public class interface and identify important design and test issues, then turn to prior art and do the 'real' implementation based on a study of the literature. Why is this approach best? Why is postponing the search for a good algorithm until after the first implementation good policy? Partly the reason is the primacy of first version as fi:amework: we want to concentrate on getting the overall structure right first. But a second and probably more important reason is to do with testing. The key problem faced by this design is algorithmic. The naive method will work, but it won't be good enough because it will
Case study: Median filtering 313 be slow and clunky. Yet a high-efficiency method is likely to be tricky to code correctly. Because we are fairly sure that our naive initial method will have little in common with the final version, it will make an excellent reference version for statistical black box testing. At this stage we can already see that once we have basicfimctionality,we will want to generate many test cases automatically. Our initial naive implementation will provide a set of results for these. We can't regard the naive implementation as the gold standard - it may have bugs but it provides a cross check.
15.7 R a p i d prototyping
We can use the framework generator on page 190 to create medianfilter.h, medianfilter.cpp, testmedianfilter.cpp and a Makefile. Then opening medianfilter.h and medianfilter.cpp in side-by-side windows, we begin by adding expected data items and fimctions to medianfilter.h. / / Header for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / John Robinson #ifndef medianfilter_H #define medianfilter_H class medianfilter { double *buf; int length; int bufptr; / / Index into circular buffer public: medianfilter (const int len); ^medianfilter (); double addnext (const double newval); / / Returns value displaced from buffer double getmedianO const; void print() const;// Debug function to show current state of / / buffer };
#endif As usual the choice of data structures affects the choice of algorithms and vice versa. We have to store the members in the set somewhere, so our first decision is to have a data item buf that will point to an array of length length. But how will items be added and removed? Assuming we want our data type to function as a filter, we should keep the current window in the object and each time we slide the window by one position of the input, we load the new item in and displace the item that's just left the window. So we need to store the window contents in a circular buffer. Therefore we need bufptr, an index into the buf array showing where the next item will be placed. We have also defined the public interface for the median filter: a constructor in which the length of the window is specified, a member fimction addnext() which adds a new value, causing the displacement of the oldest value (which it returns, which may be helpfiil in
314 Software Design for Engineers and Scientists debugging), and a member function getmedian() which simply returns the current median. Finally there is a print() function to show the current state of the buffer, which, again, will be useful in testing. What seems like a big decision has also been made - that the samples will be doubles. Indeed, the type of the samples is an important issue: since we're dealing with video coding, our samples will actually all be integers between 0 and 255, so could be represented as ints, unsigned ints or even unsigned chars. At this stage, we're fine choosing a type (double) and developing with that. Later we can make the code into a template version to handle diverse types, or perhaps have a single typedef declaration that can be changed depending on the application. The next thing we do is think about consts. Already the class declaration includes some consts, where values are passed into functions. But we should also decide at this stage which member functions leave the object unchanged and preface their declaration with const. printQ should just output current state information and should therefore be const. In our naive implementation, we would expect to code the sorting as part of getmedian() but that need not alter the object itself. Therefore, we also provisionally declare getmedianQ as const. We now move on to write our initial implementation file medianfilter.cpp: / / Implementation for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / John Robinson #include "medianfilter.h" #include
Case study: Median filtering 315 double medianfilter::addnext(const double newval) { double retval = buf[bufptr]; buf[bufptr] = newval; if (++bufptr == length) bufptr = 0; return retval; } double medianfilter::getmedian() const { / / Initial implementation: do sort each time. Must improve! double tempbuf[length]; for (int i = 0; i < length; i++) tempbuf[i] = buf[i]; qsort(tempbuf, length, sizeof(double), comparedoubles); if (length & 1) / / length is odd return tempbuf[length/2]; return (tempbuf[length/2] + tempbuf[length/2-l])/2; } void medianfilter::print() const { for (int i = 0; i < length; i++) cout « "buf[" « i « "] = " « buf[i] « endl; cout « "Current buffer index is " « bufptr « endl; } The code here is straightforward. Only two things to note. First, our window is a circular buffer, and therefore we have the test on bufptr in addnextO which cycles it back to the start. Second, we have used qsort() as discussed in Chapter 4 to sort a copy of the buffer. (If we were using the standard template library, we would have accomplished the same thing with the function sort() on a vector.) Obviously, as the code itself says, this is a clunky implementation, but it is easy to code quickly, and thus it forms our first attempt. We now develop the test framework for our new class. Initially we want to check that the circular buffer is working correctly, the true medians are being calculated, and the implementation does support arbitrary window sizes. We therefore make our initial test program an interactive loop where we enter values by hand and at each stage check that both the output and the internal state are correct: / / Test program for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 #include "medianfilter.h" #include
316 Software Design for Engineers and Scientists medianfilter testmedianfilter(filterlengh); cout « "Buffer state:\n"; testmedianfilter.printO; for (int i = 0; i < inputlength; i++) { cout « "Next input value:"; cin » val; val = testmedianfilter.addnext(val); cout « "Value displaced was " « val « endl; val = testmedianfilter.getmedianO; cout « "Current median is " « val « endl; cout « "Buffer state:\n"; testmedianfilter.printO; } return 0; }
The test program uses all the debugging features provided by the class itself to monitor what is going on. When we compile and test the program it seems to work - except for negative filter lengths. Adding numbers one by one, we keep a running hand calculation of the median, and the program agrees. We therefore have an initialfi*ameworkwhich can serve as reference comparison for our final version.
15.8 Exploit existing solutions
Our search for existing solutions begins with a search for 'median filtering' on Citeseer, Yahoo! and Google (see Chapter 11). Variations like 'median filter' + algorithm also are worth trying. Very quickly wefinda pointer to http://iris.usc.eduA^ision-Notes/bibliography/twod266.html which is a whole list of pubUcations on medianfiltering.Digging deeper draws a blank on online resources but some good leads on papers. One that sounds particularly promising is: M. Juhola, J. Katajainen and T. Raita, 'Comparison of algorithms for standard median filtering', IEEE Transactions on Signal Processing, Vol. 39, No. 1, January 1991, pp. 204-208. It's a little old, but most of the more recent stuff seems to be about hardware, and it has two attractive features: (1) it is a comparison paper, so it will be concise in describing algorithms and will give advice about how to choose an appropriate one, (2) it's only five pages long. So we head for the library and get a copy. With permission of the authors and the IEEE, that copy is reproduced in the Appendix, page 398. This paper seems like just what we want. After one read-through, it looks likely that we will want to concentrate on method D or E, since these are clearly the fastest. See Figure 4 on page 401 and Table 1 on page 402. Our next step is therefore to rough out an implementation based on one of these. When I originally attempted the design, I chose method E, because it's fastest overall, and its restriction to integers doesn't matter for the end application of video processing. I began by ensuring I understood the method, sketching out a histogram and running some numbers through the algorithm by hand, 'playing computer'. To be sure, I wrote the code for the key central section of the algorithm that involves stepping a counter up (or down) through the
Case study: Median filtering 317 number of samples that are equal to the current median as samples are added to each side. This was important because the authors' description of the algorithm didn't explicitly say what happens when a new item has the same value as the current median, and I needed to check the assumptions I was making about that. Next I rewrote the header file and implementation just to use integers fi*om lowest to highest. This introduced a lot of extra range checks so I decided to restrict the range to [0,255], since these would be the values in the video fi-ames. And then ... at last... I realized what was wrong with method E. Although the histogram method is fast, each median filter requires a histogram of bins (i.e. an array of counters) in the range of input values. Even restricting that range to [0, 255], and the window size to 255 (so that each histogram bin can be an unsigned char), means that each filter object must maintain 256 bytes of histogram. If we have a quarter million median filters running at once (for a typicalsized picture), that's 64 Mbytes of memory for what is just one component of the system. Sixty-four Mbytes is perhaps not prohibitive, but it is still a huge overhead on what method D requires. Given that D can cope with unrestricted ranges, it surely is the better choice. How long did I waste on method E? About an hour. Could I have avoided that? Certainly - right there on page 207 of the paper (page 402 of this book), the authors state that the space requirement of method E is 0 ( t / + m) where U'ls the number of integers in the range and m is the window size. All their other tested methods required space 0(m). I missed that, or perhaps I didn't register it until I'd started programming. But the cost was not too great, because developing incrementally exposed the problem more quickly than if I'd been occupied trying to write the whole of method E in one go. Method D is a doubly linked list method, and as we know, linked list manipulation can be tricky. We always have to keep in mind what happens at the ends of the lists. So here is a carefully considered implementation that treats the ends specially. First, the header file which includes the whole class definition for a linked list node called bufitem. // // // //
Header for medianfilter class, created Thu Aug 28 12:42:56 2003 Revised Aug 28 15:15 John Robinson
#ifndef medianfilter_H #de£ine medianfilter_H class bufitem { / / All items private. Just used by friends double value; bufitem ^smaller;// Using pointers rather than indices / / means bufitems can be collected into other DSs than just / / arrays, bufitem *greater; bufitem 0 { value = 0; smaller = 0;
318 Software Design for Engineers and Scientists greater = 0; } '-bufitemO { }
void set(const double v) { value = v; } double get() const { return value; } void initpointers(bufitem *const predecessor,bufitem *const successor) { smaller = predecessor; greater = successor; } void putafter(bufitem *predecessor) { / / Break existing link from predecessor to its successor / / and insert this one bufitem *successor = predecessor->greater; predecessor->greater = this; if (successor) successor->smaller = this; smaller = predecessor; greater = successor; }
void putbefore(bufitem *successor) { / / Can't just use putafter(successor->before) because / / this may be 0 on list bufitem *predecessor = successor->smaller; successor->smaller = this; if (predecessor) predecessor->greater = this; smaller = predecessor; greater = successor; } void removefromlistsO { if (smaller) smaller->greater = greater; if (greater) greater->smaller = smaller; } bufitem *next() const { return greater; } bufitem *prev() const { return smaller; } friend class medianfilter; };
class medianfilter { bufitem *buf; int length;
Case study: Median filtering 319 int nextinsertion; / / Index into circular buffer bufitem *medianloc; / / Not index but pointer public: medianfilter (const int len); ~medianfilter (); double addnext(const double newval); / / Returns value displaced from buffer double getmedianO const; void print0 const; / / Debug function to show current //state of buffer };
#endif Here is the implementation file. Note that most of the fimctionality is in addnextO as implied by the paper. Although care has been taken with pointer manipulation, this version does have a bug. See if you can track how the data structures and algorithms on page 401 have been implemented in the code, and perhaps spot the bug. // // // //
Implementation for medianfilter class, created Thu Aug 28 12:42:56 2003 Revised Aug 28 15:15 John Robinson
#include "medianfilter.h" #include
double medianfilter::addnext(const double newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length)
320 Software Design for Engineers and Scientists nextinsertion = 0; double oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf[thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) {// May need to move this node upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; } if (p->next()->get() >= newval) { buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); break; } p = p->next(); } if(!p->next()){ buf[thisinsertion] .remove£romlists(); buf[thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prevO; while(p->prev()) { if ((p == medicLnloc)&&!movedmedian) { medianloc = medianloc->prev(); movedmedian = 1; } if (p->prev()->getO <= newval) { buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putbefore(p);
Case study: Median filtering 321 break; } p = p->prev(); } if(!p->prev()){ buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); } } return oldval; } double mediaiifilter::getmedian() const { / / Initial implementation: do sort each time. Must improve! return medianloc->get(); } void medianfilter::print() const { for (int i = 0; i < length; i++) { cout « "buf[" « i « "] = " « buf[i].get(); cout « " has address " « (void *) &buf[i]; cout « ",prev pointer" « (void *) buf[i].prev(); cout « ", next pointer " « (void *) buf[i].next() « endl; } cout « "Current buffer index is " « nextinsertion « endl; }
After compiling and quick hand testing, this version seems to work. But here is where we need to do rigorous, extensive testing. Because the code is a little tricky to understand, white box testing can't be guaranteed to uncover situations that have not been properly catered for in the code. Black box testing by hand is also susceptible to missed cases: we could enumerate possibilities like add-value-greater-than-currentmedian-when-current-median-equals-value-to-remove-from-thewindow-but-more-than-one-other-value-still-in-window-equalscurrent-median, but it's going to be tedious to devise test sets that cover all combinations. This is a case when statistical testing is justified: we can generate long streams of random numbers, filter with both our initial naive method and our 'provisionally final' method and look for differences. Generating long random sequences is easy and we can even use the stable data type discussed in Chapter 17 to help analyse outputs. When we find a difference, we can check by hand which of the two systems is wrong, then work step by step using printQ (or other debugging tools) to find the problem. When we do such testing on the above implementation we discover its problem: addnextQ does not properly reset the pointer to the median when the new item becomes the median. This can be fixed as demonstrated in the version of addnext() listed below. The location of the fix is flagged by comments. double medianfilter::addnext(const double newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length)
322 Software Design for Engineers and Scientists nextinsertion = 0; double oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf[thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) { / / May need to move this node / / upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { / / THE NEXT 11 LINES HX THE BUG BY REORDERING //TESTS / / A CORRESPONDING 1 l-LINE FDC IS FLAGGED //LOWER DOWN if (p->next()->get() >= newval) { buf [thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->next(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; } p = p->next(); } if(!p->next()){ buf[thisinsertion] .removefromlistsQ; buf[thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prev(); while(p->prev()) {
Case study: Median filtering 323 / / THE NEXT 11 LINES FIX THE BUG BY REORDERING //TESTS / / A CORRESPONDING 11-LINE FIX IS FLAGGED //HIGHER UP if (p->prev()->get() <= newval) { buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->prev(); break; } if ((p == medianloc)&&!movedinedian) { medianloc = medianloc->prev(); movedmedian = 1; } p = p->prev(); } if(!p->prev()){ buf [thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); } } return oldval; }
15.9 Finishing
Now extensive random number tests succeed for the implementation. We therefore consider what should be done to finish the program. The helper class bufitem is only used inside medianfilter. It could be an embedded class, but at the very least it can be moved out of the header file and into the implementation. How shall we make our median filter work on other types than doubles (we're particularly interested in ints for video processing)? One possibility is to use templates, but because export is not supported by many compilers, that would mean relocating a large part of the implementation to the header file. Instead we use a typedef in the header file with a comment to remind us that this has to be changed and things recompiled for different types. We also need to fix the problem identified in the very first version: a negative argument to the constructor causes errors. There are a number of ways of dealing with this, but the best option is to signal a problem but proceed anyway with a default value replacing what the user asked for We need a default value for the constructor too, so that arrays of medianfilters can be instantiated. Having made all those changes, here is the final version. First the header file: / / Header for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 //Revised Aug 28 //Version 0.1 / / John Robinson #i&idef medianfilter_H #define medianfilter_H
324 Software Design for Engineers and Scientists typedef double itemtype; / / Change and recompile for other / / input sample types class medianfilter { class bufitem *buf; int length; int nextinsertion; / / Index into circular buffer class bufitem *medianloc; / / Not index but pointer public: medianfilter (const int len =15); -medianfilter (); itemtype addnext(const itemtype newval); / / Returns value displaced from buffer itemtype getmedianQ const; void printQ const; / / Debug function to show current state / / o f buffer };
#endif Note the addition of the class keyword before bufitem's first mention, because now we have taken the declaration of bufitem out of the header. Here is the final implementation file: / / Implementation for medianfilter class, created Thu Aug 28 / / 12:42:56 2003 / / Revised Aug 28 ... Corrected same day //Version 0.1 / / John Robinson #include "medianfilter.h" #include
void set(const itemtype v) { value = v;
Case study: Median filtering 325 } itemtype get() const { return value; } void initpointers(bufitem *const predecessor.bufitem *const successor) { smaller = predecessor; gfreater = successor; } void putafter(bufitem *predecessor) { / / Break existing link from predecessor to its successor / / and insert this one bufitem *successor = predecessor->greater; predecessor->greater = this; if (successor) successor->smaller = this; smaller = predecessor; greater = successor; } void putbefore(bufitem *successor) { / / Can't just use putafter(successor->before) because / / this may be 0 on list bufitem ^predecessor = successor->smaller; successor->smaller = this; if (predecessor) predecessor->greater = this; smaller = predecessor; greater = successor; } void removefromlistsO { if (smaller) smaller->greater = greater; if (greater) greater->smaller = smaller; } bufitem *next() const { return greater; } bufitem *prev() const { return smaller; } friend class medianfilter; }; medianfilter: :medianfilter(const int len) { if (len < 0) { / / Can add an upper limit later if necessary cout « "Specified length" « len « " is out of boundsXn"; cout « "Using length = 15 insteadXn"; length = 15; } else
326 Software Design for Engineers and Scientists length = len; buf = new biifitem[length]; / / Have to treat the items at each end differently to the others / / because they are at the ends of the Unked lists buf[0].set(0); buf[0].initpointers(0, &buf[l]); for (int i = 1; i < length-1; i++) { buf[i].set(0); buf [i] .initpointers(&buf [i-1 ] ,&buf[i+1 ]); } buf[length-l].set(0); buf[length-1 ] .initpointers(&buf [length-2] ,0); nextinsertion = 0; medianloc = &buf{length/2]; / / Actually doesn't matter because initial / / buf values are all equal to 0 } medianfilter::~medianfilter() { } itemtype medianfilter::addnext(const itemtype newval) { int thisinsertion = nextinsertion++; if (nextinsertion == length) nextinsertion = 0; itemtype oldval = buf[thisinsertion].get(); if (newval == oldval) / / Everything fine. Nothing to do return oldval; bufitem *p = &buf [thisinsertion]; p->set(newval); int movedmedian = 0; if (newval > oldval) {// May need to move this node upwards if (!p->next() | | p->next()->get() >= newval) / / Again, don't need to change return oldval; / / Otherwise have to move upwards if (p == medianloc) { medianloc = p->next(); movedmedian = 1; } p = p->next(); while(p->next()) { if (p->next()->get() >= newval) { buf[thisinsertion] .remove£romlists(); buf[thisinsertion] .putafter(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->next(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->next(); movedmedian = 1; }
Case study: Median filtering 327 p = p->next(); } if (!p->next()) { buf [thisinsertion] .remove£romlists(); buf [thisinsertion] .putafter (p); } } else{ if (!p->prev() | | p->prev()->get() <= newval) / / Again, don't need to change return oldval; / / Otherwise have to move downwards if (p == medianloc) { medianloc = p->prev(); movedmedian = 1; } p = p->prev(); while(p->prev()) { if (p->prev()->get() <= newval) { buf[thisinsertion] .removefromlists(); buf[thisinsertion] .putbefore(p); if ((p == medianloc)&&!movedmedian) medianloc = medianloc->prev(); break; } if ((p == medianloc)&&!movedmedian) { medianloc = medianloc->prev(); movedmedian = 1; } p = p->prev(); } if(!p->prev()){ buf[thisinsertion] .removefromlistsQ; buf [thisinsertion] .putbefore(p); } } return oldval; }
itemtype medianfilter::getmedian() const { return medianloc->get(); } void medianfilter::print() const { for (int i = 0; i < length; i++) { cout « "buf[" « i « "] = " « buf[i].get(); cout « " has address " « (void *) &buf[i]; cout « ",prev pointer " « (void *) buf[i].prev(); cout « ", next pointer" « (void *) buf[i].next() « endl; } cout « "Current buffer index is " « nextinsertion « endl; }
328 Software Design for Engineers and Scientists Is this the end of the story? Not quite. addnextQ can be rewritten to avoid repetition. The two cases (newval < oldval and newval > oldval) can be combined by adding extra functions to bufitem that switch on a parameter to say whether the search is going upwards or downwards. Making the change certainly tightens up the code, but the clarity of shortness and saying something only once is traded for more variables being maintained and repeatedly checked. The upshot is an execution time increase of about 25%, so such a change is probably best avoided. Finally we install the module in the system where it will solve the problem. The median filter developed here was integrated in a video processing system, where it achieved an order of magnitude speedup over a previous median filter that sorted the window each time.