Lecture 12: Blob Analysis, Binary Image Processing, Green's Theorem, Derivative and Integral

Key Takeaways

there's quiz one which will be out today and the rules for that are slightly different from the homework problems it counts twice as much as a homework problem and it's longer not twice as long but something like that and so you have also a little bit more time to do it it's not just one week but um it's a bit of work so don't leave it too late please start on it fairly soon and this is where you're supposed to work on your own uh not not collaborate so okay and it covers uh you know whatever we did up to this point uh with a bit more emphasis on recent material and i guess only the last question has to do with the patent we're discussing so i should be able to deal with the other four right away okay a little uh sidebar here about um intellectual property and starting off with the fact that some people don't like that terminology and um you know uh how can an idea be property a little bit like uh how can a company be a person well it's you know an awkward thing in law that you have to come up with some way of formalizing these things and so in this case ideas are treated as property to some extent there are several different types of intellectual property first of all we have patents which is what we'll be talking about and they're different types of patents the so-called utility patents i.e useful and not to say the others aren't useful but their design patterns so for example jerome little lamelton's first patent was a basically a design patent for a baseball cap with a hole in it and a tube that you could blow through so this is the inverse of the beer baseball cap and there'd be a propeller and as you blow through it the propeller would turn another famous design patent is apples designed for a cell phone and their design pattern says that it should be a rectangle with rounded corners and there was a lawsuit with where they sued samsung because their phone you know with a rectangle with rounded corners and the jury awarded uh apple a billion dollars so big money involves here and you know you might say well what should a cell phone look like other than this you know this is it but it's not completely crazy because at the time they invented it you know phones were these big concrete things you held in your hand and so anyway uh not not to argue about the merits of this and and of course it's been appealed and so in the end i don't know what will happen i guess some some large amount of money will transfer uh change hands an amount of money which is astronomical to us mere mortals but probably for these companies you know it's a small change so we shouldn't feel too sorry for them so patents we mentioned last time this was the kind of social contract where you explain exactly how to do something in return for having a limited monopoly that lasts a certain number of years and the you know the rules change slightly as as time goes on uh by the way um you are supposed to explain uh how to do it and in particular if you know a good way of doing it and in your patent you only describe a lousy way of doing it that could be grounds later on for invalidating a patent and that's called best mode so there's lots and lots of terminology and in the litigation there's sort of standard things that that come up one of them is best mode so in in one machine vision patent for a wheel alignment so that one has cameras and leds and it determines the axis of rotation of the wheel and the axis of the steering that you it turns around when you turn your steering wheel and uh the pattern that was claimed to have been infringed uh did not uh disclose the best mode it disclosed the method that if you actually implement it uh it gives gives you the correct answer if you have perfect measurements but with uh realistic measurements it would have an error of you know a degree or two and that's not good enough you know if you have a car like a bmw they measure they specify those angles to 0.01 degrees so anyway so they had a problem there because they did not disclose the best mode okay utility patterns and we talked about the structure of those i'll just briefly talk about some of the others so copyright um you know you can protect artistic expression using copyright so if you write a book there's a copyright on that if you record a song there's a copyright on that if you choreograph some ballet as a copyright on that and their exceptions so for example if i want to talk about some topic i am allowed in class to present a portion of that material without violating the copyright without having to ask the author you know how much do you want for me to use your product there are also exceptions where you use extracts so for example if you're writing a news article or say you've got a blog about a movie and there's a particular part of the movie that you'd like to talk about under certain circumstances you can clip that out and use it and of course you can imagine there's a lot of fun that lawyers have with this because a lot of music is extracted from other music and put together and is that legal or is that copyright infringement copyright used to be basically for the lifetime of the author you know if you write a book you should benefit up to a certain point when you're not around and you don't benefit well of course the heirs of famous authors didn't like that so that was changed and um now we have the the sony bono uh system of copyright basically he was the person who convinced congress to change the laws and now it's uh author's life plus 75 years and they've been periodically updating it so it's basically author's life plus infinity because every time we get close to the limit they say oh well let's upgrade these poor heirs they can't live without the royalties from from this so excuse me if i'm being sarcastic so okay so that's copyright and by the way that was very important for a while uh and still is uh in software because again the rule was you couldn't copy you couldn't patent uh mathematics or or an idea abstract idea uh and so a lot of uh the courts typically held that you know that's what programs are their ideas they're abstract stuff they're not you can't hold them and weigh them and so uh the way people protected themselves was to send their programs into the copyright office and register them there and you can imagine that you know say you have some operating system like ibm 360 operating system and you send in your 750 million lines of code to the copyright office that was quite quite an exciting event and certainly if someone makes an exact copy of your program then you can sue them under under this law of course uh you know with all of these they're people who are working to get around these laws and so in the case of copyright there was this notion of a clean room so what you would do is you'd have a bunch of people that understood programs well analyze what someone else's program is doing kind of reverse engineering it but they weren't allowed to show it to another group of people they were able to they were only able to tell them what what it does and then this other group of people would write the code and supposedly because they didn't see the original code there was not a copyright violation so okay um not sure that's you know and it's not totally unreasonable in that if you look at the code line by line it's likely to be totally different you know they'll use different variable names and stuff so okay then there's a trademark so trademark is um much more restricted it's just you know if you want to call your company dunker donuts you can uh trademark that and but it has to be uh unique in the field and it has to not be you know there's a whole bunch of rules but basically you can't use a common word like you can't call your company time space because that's those are common words and so a lot of company names are slightly misspelled versions of common words the other thing is that the trademark may include particular shapes particular markings so you can distort some later and color so you know dunker donuts has two colors that they've copyrighted and i'm proud to say i have hats that have both of those colors they're very useful in hunting season which is approximately what's happening right now where i live so you can use color to protect yourself [Music] they have to be unique to the field so you may have a trademark on the name in the rubber industry uh and someone else has it in the semiconductor industry they don't conflict and there's a famous case about that where again apple apple sued the beatles why well because when the beatles started they formed their own company called apple may not know that but in any case so apple sued them and they lost because it's they're not there's no confusion of the client the customer right they're in two totally different fields one is music and the other one is computers and of course apple claimed oh well but we distribute music so it's it is the same thing well whatever anyway um then the last thing is trade secret which is you know no protection at all you just hide what you're doing you don't tell anyone exactly how you're doing it and you know classic example of that is coca-cola there's a safe down in atlanta georgia which has the formula for the ingredients in coca-cola and supposedly you know very few people know what's what's in there and uh the danger of course of i mean the good part is that it's unlimited i mean there's no like expiration date like on patents and the bad part is you know if it ever comes out there's no recourse there it is everyone now knows what to mix up to make coca-cola so um it's a risky uh thing but it's certainly much cheaper than pursuing some of the other avenues you don't have to pay lawyers for it and there is a certain legal recourse you know if somebody signs a non-disclosure agreement when joining your company and then they walk off with a formula for coca-cola and tell pepsi caller how to do it then you can sue them but you know the cat's out of the bag it's not so easy to recover from that kind of loss okay just a quick um overview of what's called the intellectual property and you know richard stallman would definitely complain that i use that terminology because he doesn't think that ideas should be treated as property just as some people don't think that companies should be treated as persons in under the law but there you are okay back to the particular patent so this particular patent is kind of real low level i wanted to start with something very simple you know finding edges very accurately finding edges very accurately so where are we going with that well the idea is that once we found edges and we're describing uh images using edges we can do more elaborate things like recognition and determining the position and determining the attitude of an object and most of uh what we'll be doing after finishing with this pattern is in in 2d where you know the world of course the 3d but there are lots of cases where 2d machine vision is incredibly successful and um you know one thing that's important to point out is that it has to be incredibly accurate you know like if if your thing works 70 of the time forget it it's got to be working 99.999 of the time so these methods are you know very carefully thought out and uh tuned to have extremely good performance and what we'll find is that in the 2d world this is possible it's a little harder once you get to 3d so once we've done position and attitude in 2d we'll progress to 3d which is a more interesting problem but you know let's get some of the basics sorted out so i'll just sort of quickly review what this patent did and so we'll start so this is and the first idea was to look at the brightness gradient and so just a quick story on that so here we've got you know brightness as a function of x and there's a point where the curvature becomes zero and changes sign and that's called an inflection point and that's what we're looking for so the edge is actually spread out over several pixel but we're trying to identify very accurately a particular point on the edge then we looked at the derivative and there we're looking at a peak we're searching for a peak and some methods use second derivative and look for zero crossing but importantly this is in the gradient direction so this is of course a 1d cross section and we're dealing with images so how do we take the cross section well we're only interested in the direction that is perpendicular to the edge and so these graphs and these ideas correspond to that okay then we talked a little bit about things that are sometimes called stencils and sometimes called computational molecules and what we're trying to do is in the discrete world we're trying to estimate derivatives and so of course um you know there's some obvious ones so there's there's a way of getting an approximation for e sub x and here's another one and these are all some of the ones that we already looked at so there are lots of ways of estimating the derivatives and how do we know which one to choose well uh their trade-offs and we saw that one way to progress is to use taylor series expansion to see what the lowest order error term is because we the higher we can push that if it's a third order error term uh it's better than if it was a second order error term and if two methods have the same order of error term then we look at the coefficient and if it's lower you know like like these two have the same uh lowest order error term and also what was very important was you know where are we trying to estimate the derivative and we decided that we get the best results at those particular points now some of those points are offset by half a pixel from the pixel grid and that's why people don't use them but that's silly because what's wrong with having some quantities on the pixel grid the image and some quantities on a grid that's offset by half as long as you know that it's offset by half you can translate any result on that grid into a result on the other group we can also analyze these in the fourier transform domain in terms of how they affect higher frequency content and we won't do that right now but that that's second set of methods we can use to decide which which of these is better now these derivative estimators can become quite complicated if you're looking for high precision so for example in in analyzing uh muscle electrical signals there's quite a bit of noise and there's quite a lot of [Music] distortion of the signal as it progresses through uh tissue and so you know when people are trying to estimate the first derivative they in one case one paper i saw use a 39 degree approximation so they use a pattern like this that's you know 39 elements long and they feel that that way they can control the trade-off between suppressing the noise and getting very accurate results so we're using very low order partly because you know it's easier and partly because one other trade-off is if you make them too long then they start to have different features interacting with each other so here we're trying to detect an edge now if you use an edge operator that's 100 pixels long and there's another edge within that 100 pixels it's going to interfere with the results so we we try to compromise on the one hand we get better results with bigger support better noise suppression but at the same time we run more into the problem of what happens when two edges get close to each other and in particular you know here's an image of the corner of a cube and yeah over here you know we can we could have potentially quite a large support but when we're back here edges get pretty close together and then a large support means that you're combining information about different edges and you won't get the results you you would like okay that's first derivative and if our model is we find the peak in the gradient direction of the first derivative then that's all we need but for some purposes we might want higher order derivatives and we'll see some examples of that later so now one way to think about this is that second derivative is just the first derivative applied twice and so we can run our numerical approximation a second time and that corresponds to convolution and the result is so that we can easily compute second order derivatives this way and let me just make sure we understand what this is so these computational molecules the way they work is that you put them down on the image grid and then you multiply the gray levels by whatever the weight is at this point so multiply that by one put it in accumulator take that one multiply it by -2 add that to the accumulator take that one multiply it by one add that to the accumulator and then finally take the result and divide by epsilon squared well often we don't care about constant factors so we might drop that but if we're actually thinking about derivatives we should we should do that so that's one thing and we're really talking about convolution right so we apply this to every place in the image slide it along if you like to produce a new image results at a bunch of points now how do we know that you know if we want a sanity check well a way of making this making a check on this is to try it on some function where you know what the answer is so right so we know that in this case uh the answer should be 2. so so let's apply this and well we have to decide where we're applying it so let's suppose that this is uh zero this is one this is minus one so we're plunking this down and so what do we get well uh in this case epsilon is obviously one so that's one times and then 1 times -1 well that's 2 no sorry 1 times yeah 1 times -1 squared is 1. then we have minus 2 times x squared with 0 and plus 1 times 1 squared and so that's 2. okay so apparently it works so you know the different ways of designing these computational molecules um certainly convolution is one and then you'd want to check them so the another thing you might want to check is well what if f of x is x uh okay well then we're going to get 1 times -1 0 1 times 1 the answer is 0 which is what it should be right what if f of x is one so you want to check it for polynomials of low order up to the order of the derivative you're trying to get so if f of x is one then of course we get one minus two plus one is zero which is again what what you what it should be because the second derivative is zero so one of the constraints on these operators if they're supposed to be derivative operators is that the weights add up to zero that makes sense okay uh well so that's just one uh one dimension so to speak what about d2 dxdy we can just do that and so right so because this is the x derivative and that's the y derivative and uh composition of differentiation corresponds to convolution so that's a we're kind of sneaking in some stuff here one of them is that you know uh we're used to dealing with linear shift invariant systems but it turns out that of course derivative operations are linear shift invariant in that if i take the derivative of a function and then take the derivative of the same function shifted what i'm going to get is the derivative shifted right if i take the derivative of the sum of two functions i'm going to get the sum of the derivatives so that's a very important thing to exploit is that um you know derivatives taking derivatives can be considered convolution and so all of the good stuff we know about that applies so let's see what this gives us now um uh what we're going to do basically is take one of these and uh flip it and superimpose it on the other one and if i superimpose it over here of course i get zero because i'm multiplying you know we assume the background is zero we assume that the only values we're showing are the non-zero values okay but then what if i move it here well then there's an overlap and then i move it over there there's an overlap and i can move it down here there's an overlap and i move it down here and again so i expect to get four values out of this convolution and so i'm going to get something like this and so my 2 by 2 stencil for estimating the mixed derivative is is just that and it makes perfect sense one thing to watch out for is that in convolution you know we flip one of the two functions and uh then then slide it across and try it in all possible places uh over here sometimes we're using computational stencils that are not flipped so we may end up with some sign reversal but you want to you can always check on a simple polynomial to see whether it's working the right way around by the way this one here we can sort of take a diagonal view of this it's a little bit like right if i project this down so there's a plus one far out on this side there's a plus one far out on that side and then there's a minus two there two minus ones in the middle and this looks awfully like a second derivative uh you know like the one we had over here just rotated 45 degrees and of course this mixed partial derivative is uh in a rotated coordinate system just d2 dx squared right so if i take so this is my original coordinate system and then i look at the world in this coordinate system then the second derivative e x prime x prime is the same as that mixed derivative so you know we kind of think of e x y as a very different animal from e x x and e y y but but it isn't oh yeah just a little bit more of this before we go on um so and we already mentioned the laplacian let's just bring that in here again so we had a second derivative operator here and so in the case of the laplacian we're adding two of them in different in orthogonal directions so we could just do right so that's just this thing plus it the same thing rotated 90 degrees and so so that's one way of writing the laplacian but that turns out to um also be a good approximation uh to the laplacian how do i know well there are lots of ways of checking one is taylor series another way is apply it to test functions see if it gets the right result and another way is fourier transform but you know those and and you'll see that this is really the same pattern rotated 45 degrees but now the separation is square root of two times as large so i end up having to change the the weight and now i mentioned that the laplacian is the lowest order linear differential operator that's rotationally symmetric well neither of these looks particularly rotationally symmetric so can we make one that's a bit more i mean on a square grid you can't but we can do better than these and one way is to combine them to take a weighted sum of these two okay so i've taken uh i don't know four times this one plus one times uh this one and and i get this operator and it's a little bit smoother in terms of rotational symmetry uh the the corner ones aren't zero uh but they also don't have the same weight as the diagonal one the the up and down ones um because they're further from the center so how do i know that's a good combination how did people get to that one well again you look at the taylor series and it turns out that the lowest order error term for this one is one larger than the lowest order error term for either of these two so it's better uh more work right if you do if you apply this to your image uh you're doing not quite twice as much work but you're doing more work because you have to take into account all of the corner ones which this this operator doesn't or conversely the up and down and west and east ones uh there so and you know if we wanted to we could do a more detailed analysis of the error terms but that's that's pretty boring so we won't do that anyway so those are kind of all of the computational molecules we're going to need and as you can see the the options the the different versions uh that are trying to compute the same sort of thing of course on a hexagonal grid this looks much better because we can we can do something like this which uh you know looks nice and rotationally symmetric yeah what's the one in the middle oh -20 so so can anyone think of how you could determine that it should be minus 20 without actually doing a lot of area algebra they have to sum to zero right why is that well because when you apply this operator to f of x equals one the laplacian is zero and so if you apply this to the function that's one everywhere obviously you'll just be adding all the weights and so they they need to add how do i know it's one over 6 well one way to get it is is to apply this to a test function like you know if f of x is x squared plus y squared where i know that the answer is 4 and then you know go from there or i can get it from this argument which is the weighted argument that i'm taking uh four times one of those and one of those and this one has a factor of 2 in there so you end up with i don't know a half plus a quarter no anyway it comes out to 6 epsilon squared okay so so it's sort of annoying that we don't have hexagonal pixels by the way there are some situations where people are very concerned about [Music] you know efficiency you know like trying to image the black hole at the center of our galaxy using using radio frequencies by you and you know each antenna you put up costs a pile of money so you want to make sure that this grid of antennas is the most efficient way of sampling the fourier transform space in that case and so it turns out that this way of sampling is four over pi as efficient as that way of sampling so there are places where people do this and and for example briefly in you know almost all chips are laid out on a rectangular pattern because that's very easy to do and check but if it comes down to packing density and you know particularly if you have something that has a very simple repeating pattern then sometimes there's an advantage so there were memory chips for a while that used the hexagonal layout but they've since disappeared because now we're stacking things vertically and right now it doesn't seem to be an efficient way to go um oh so while we're here i should mention i mentioned already that the laplacian is the lowest order linear operator differential operator that's rotationally symmetric here's a non-linear one so this operator is rotationally symmetric what do i mean by that i mean that if you rotate the coordinate system you know for example over here and you compute e sub x dash and square it and add e sub y dash and square it you'll get the same answer as if you take so so it doesn't depend on the orientation of the x axis and so this is lower order of course in laplacian because it's first order but it's non-linear so nevertheless we run into this quite a bit because remember uh roberts use these stencils which you can just think of as ways of computing the derivatives in the rotated coordinate system and then he took the square root of the sum of squares and that was his edge detector and so it's equivalent to doing this and for his purposes that took less computation and he already knew in 1965 that what you want to do is make sure that your ex and ey refer to the same point in pixel space a lesson which has since been forgotten so except here um okay so that was the front end and this has to be very efficient because it's run on every pixel and it also lends itself to special purpose hardware of course so the next step was our subpixel edge detection method so we used what is called non-maximum suppression so um this weird terminology i mean why not just say finding the maximum but there it is so what where did that terminology come from well the idea is that we apply this edge operator everywhere and in most places it has a kind of feeble response but on the edges it really kicks up and so one approach would be you know let's just threshold so if if we get a strong response from one of these molecules here then that's then we're on the edge and if not then then we're not well that involves um early decision making because once you've made that decision uh that's it you'll you you know you you'll throw away you'll throw away that edge point and never never do any computation with it again because you've decided it's below threshold or vice versa you picked it as being a threshold so a big in the patent bill silver makes a big fuss about avoiding thresholds if possible and not making decisions too early that's his sort of main motivation for not using thresholds so some previous edge detection work did work that way you apply some operator which has a strong response on the edge and then you threshold and now you get responses but they're not just right on the edge because we saw that the edge is a slow smooth transition so there'll be neighboring points that also have a strong response plus with noise there will be points in the background where noise just happens to add up and now there seems to be an edge there so that's undesirable okay so so the previous methods worked by thresholding and looking for things that had a strong edge response and here instead what we're doing is we're going to remove everything except the maximum but maximum in what sense well again just in the gradient direction and so here here we deal with the unfortunate fact that we're going to quantize the gradient directions so we're only going to have you know compass directions east northeast north northwest west etcetera eight of them and let's suppose that this is our quantized gradient direction then we step through the image and we consider those three values and the non-maximum suppression says that we will accept this as a potential edge point only if this is true right because if g minus was bigger than g zero well then that's going to be the edge point and remember that the edge is running at right angles to this lot so the the actual underlying edge is is is like this and we're looking in the gradient direction and the gradient direction is of course perpendicular to the edge okay and then we talked about how there's this asymmetry because occasionally we're going to find that g 0 is equal to g minus or g 0 is equal to g plus and we don't want to declare both of them to be edge points we want to have a tie breaker so that only one of them gets elected and which one is well it's arbitrary you can you know you could easily well have done that as long as there's a way of breaking that tie okay and then we said okay now we can plot the profile of this edge response along the gradient direction and we get a picture like this and then we can fit some curve to it for example we can fit a parabola to it and we find the peak of the parabola and that's our sub pixel edge position so so several points there one of them is you know why fit a parabola well that's arbitrary uh you know the the shape of that curve depends on the optics the image sensor the thing you're looking at and but we only get three samples of it so treating it as a smooth curve the lowest order polynomial that will work is second order so uh one option we have is to use that and not to say that that's the only option or the best one but it's you know pretty good guess uh next thing is that what what's s well s is the displacement from g0 so in here s0 would be that point and then if s is say well if we get over here then s is a half no so if we get over here then s is a half right we're halfway to that point and obviously it doesn't make sense for s to be bigger than a half because then this would have been the maximum right so and same in the other direction notice that s is not in units of pixels because in this diagonal case s equals 1 is that distance which is the square root of 2 times the pixel spacing right whereas if we happen to get the quantized direction over here then s would be the pixel spacing so um that's something to keep in mind that it all depends on the actual gradient direction okay um avoid thresholds yeah and so now um we have uh a potential edge point uh and notice we're not doing any thresholding so we get we're gonna get these points all over the show not just on on the edge we haven't done any thresholding so here you know we we mark uh this place based on uh square root of 2 s delta where delta is the pixel spacing so just to be precise so that's that's where we think the edge is so i've drawn the edge here but actually now with sub pixel interpolation we find it's it's there so if we just were to go with the peak in that curve uh on a discrete grid then we'd find that the edge runs through that point but now we know it's it's over there okay but in order to get there we will quantize the we quantize the edge direction and so we can improve things slightly right so i suppose i should draw that again and make it less messy so here's our quantized direction and here's some point we found that's uh an edge point and then suppose that the actual gradient direction is slightly different it can't be hugely different because then we would have picked a different compass direction okay so now um that's the gradient direction so the h has to be perpendicular to it so i can draw a perpendicular to the green line um [Music] that passes through this point so the edge actually is like this and when i report a potential edge point i can report any point on on that line and which one do i pick well the simplest one is just to pick the one i calculated i can however slightly improve the result by uh instead picking this one why well it's closer to the origin it's closer to the actual peak and so it's less likely to be subject to noise it's not a huge improvement but you know they decided this was a worthwhile additional step it also actually aids in a step that we will be getting to in a second where we chain together the edges okay so that's the idea so we project we project from the quantized uh gradient direction down onto the actual gradient direction and that's the plane position step what next o so then we get to the bias compensation okay so we we said that we somewhat arbitrarily picked the parabola and then as you know in the patent there's a second method which uses a rooftop a triangular shape he uses that it fits that and finds the peak of that and and that gives you another possible position for where the edge really is and we said that in certain circumstances one may be more accurate give you a more accurate answer than the other but what you can imagine doing is actually you know experimentally moving the edge by very small increments tiny sub pixel increments and seeing what your method gives you and then plotting that against what it should have been right so this is an experiment where you have the camera looking at an edge and then you move the edge or the camera by some tiny amount in increments and and you measure now ideally you know your magic method for peak finding should always give you the correct peak value but it may not so ideally you get a 45 degree line and again we're only interested in the range from minus a half to plus a half and now all of these methods should give you the correct answer in three cases no matter what you do if you act if your peak is actually a g zero then it should return g zero that position s equals zero and similarly if you were halfway between pixels it should give you the and both of these methods do they satisfy that requirement but what happens in between well uh as i mentioned that depends on you know exactly what shape the edge has so and you might end up with you know something like this or maybe something like that so typically the departures from the diagonal will be quite small and typically they'll be quite smooth and so you know you could keep a lookup table of this or something to compensate for it but it's it's not really worthwhile it's it's a relatively small correction again but in going down to the the aim is uh 1 40th of a pixel accuracy so you know if you think about it that's that's quite amazing that you can do that but and one part of it one part of it is this plane position correction and one part of it is this and so what what is done is to approximate whatever that shape is with this function and so for example for b equals 0 we get s prime equals s so that's just the diagonal line that's the ideal case for b greater than 0 then that means that s prime is uh s raised to a power greater than one and so that means it's going to bow uh this way as as the red curve does and if b is less than zero that means that s prime is s raised to a power that's slightly less than one so that's like approximating square root so that's this uh bowed upper upper curve which i'll show in green so that's green and the other one is red so what's the point the point is that this is a small correction so it doesn't really pay to try and be too fussy but you want to do it and this is a one parameter fit to the type of curve that what's the b is the one parameter and so uh you could calibrate it based on this method you know get one value of b you could calibrate it based on this method get a different value of b then you notice that if we are happen to be working east-west the spacing between pixels is much less than if we're working northeast and so it turns out you want a different value of b for for this case than you do for that case fortunately the only two cases you know they're the ones that are east west north south and then the ones in between but you you can use a different value of b for those two okay and you know again you could do something more elaborate but since it's a small correction it's only going to affect a small fraction of a pixel also um you don't want to be too clever here because this curve is going to change a little bit with circumstances you know if if your camera is slightly out of focus you'll get a slightly different result if the corners of the cardboard box you're looking at are somewhat damaged then you'll get a different edge response and so on so you don't want to be too overly clever okay now a lot of this depends on the actual edge transition and you know we we drew one just by hand and uh then uh came up with these methods for finding out where the edge actually is but uh realistically what you know what's causing these edges to be fuzzy and we already said that it's a good thing they are fuzzy otherwise we wouldn't be able to um do the sub-pixel recovery we'd be suffering from uh huge aliasing problems so it's it's a good thing they're fussy but you know why are they fuzzy well one reason is the defocus so let's just look at that uh as a special case that's of interest for other reasons so here's our lens and you know the objects up there and it's suppose that it's it's a point light source and this is the in focus plane where that distant light source star maybe is image as a point but our camera has the image plane slightly off and so the picture will be slightly out of focus so what did i call this so this is f and i guess i call this delta so [Music] when i look at the image of that star it's no longer a impulse a point it's a circle so the uniform uniform brightness and if i want to plot it as a function of x and y it would look like this and i don't know i call this a pill box i guess people don't have pill boxes anymore but it's a cylinder of constant height and if i want to describe it mathematically so what's that big r is the radius of of the little pill box and i divide by 1 over pi r squared because the same energy is being deposited into that area no matter how out of focus i am so if i'm in focus it's all that concentrated at one point if i'm out of focus that same amount of light is spread into a larger and larger circle and so i compensate for that by dividing by the pi the area of the circle and then what's this so this is the unit step function and so for r equal to zero this will be u of sum minus quantity so that will be zero so i get one minus zero and so this will just be one and then when i get out to the radius if i go past the radius if little r is bigger than big r this will be equal the step function is one right because this will be greater than zero and so i get one minus one is zero so it's it's just a fancy way of saying the same thing as that diagram uh the other thing i need is you know what is r what is big r how big is it well it's obviously going to depend on how far out of focus i am so it's going to depend on delta and it's going to depend on the size of the lens something like that right this is just similar triangles again okay um oh a f uh which is this distance here from the lens to the uh image in focus image plane okay so obviously as i go more out of focus as i move the image plane further up this radius gets bigger and the brightness gets less because i'm now putting the same energy into into a larger area and so this is uh called the point spread function psf point spread function for this system when it's out of focus and [Music] this is used a lot in understanding the effect of being out of focus you know we think of it as blurring and of course we can think of it in terms of the fourier transform in terms of removing some or suppressing higher frequency content and it's the higher frequency content that makes things look sharp and in our case blurs the edge so let's see what the effect is on on the edge well basically we're going to have a response which we can calculate geometrically by superimposing the edge and the circle so let's take a simplest case where you know there's black on one side of the edge and white on the other and um now we have a radio a circle of a certain radius and what are we looking at well we're looking at this overlap you know that's going to be what controls how bright things appear and then we're going to move this thing across the edge to see how the response varies so we take this disc and we slide it across and obviously until it touches nothing happens there's no output and also obviously once we get over here the output is constant it's one nothing more changes so there's a transition between x being minus r and x being plus r where there's some sort of change between zero and one and and that's what we're trying to calculate well unfortunately uh it's not quite that simple so we could just write the answer out by inspection but we can't get it this way okay you know there's probably a formula for a sector of a circle like that but i don't know what it is so let's do it uh using things we do know so a couple of things we know we know how to compute the area of a triangle well probably know several ways of computing the area of a triangle and we also know how to compute the area of a sector of a circle so let me draw that here so that that's the sector of the circle that's the triangle i'm looking at and it's obvious that the quantity i want is the difference right it's what's left over when i subtract that triangle from that area so this is r and i'm going to call this angle theta okay and this bran this thing here is obviously r squared minus x squared so x is the position x0 means that you're right on top of the edge that the circle is bisected by the edge and then x can get as large as plus r before it saturates or minus r before it saturates okay so the area of the sector is 2 r squared theta write this this sector here and we can check you know when theta is 2 pi um hmm okay i guess oh this is theta can only get up to pi so uh we would get uh oh pi r squared how about let's check that again so theta is the half angle of the sector and so it can only get as large as pi and if it's equal to pi that means we've covered the whole circle and we should get the area of the circle which is pi r squared right so okay so we get r squared theta and then we have to subtract from that the triangle area which is right it's the half the base times the height so the base is uh two times this quantity right because this quantity is just that section so that's the base and then the height is x and so half base times height gives me that and [Music] theta is given by this quantity and so i end up with i mean the the details aren't too exciting here i'm not going to expect you to remember this but we can plot that and see what sort of transition it gives us and what's more we can then feed it into this algorithm of this pattern to see how accurately it will determine the edge position and in particular we can plot the diagram that i just erased here which had s versus s bar right we in other words if this is how the edge is formed if this is why the edge has that smooth transition then this will allow me to calculate the error and the error ought to be pretty small but it's non-zero and if you want high accuracy you have to take it into account so another way of looking at this is to plot this diagram so so so where does that come from that that looks like a circle except it's um elliptical because it's been increased in height by two well if you think about it when i move this uh edge or the circle then i am adding or subtracting an area you know infinitesimal area that has a height equal to this quantity well in other words just you know twice the height of the circle there so actually i don't need for some purpose i don't need to do all that hairy math i can just look at that diagram and immediately write down uh that the brightness derivative is that and and oh it has a peak at zero well we expect that okay while i'm here i can look at the second derivative and it looks like that and it's x over minus it's just the derivative of that which is pretty easy but what's e well it's the integral of this thing and and we we just we just did that integral in a somewhat painful way but but it's probably just as well because i don't remember what the integral of that is and you know there it is i think okay okay so uh what can we do with this we can now feed this into the algorithm and say if this is how why the edges smoothly varying because it's out of focus then this is the relationship between the true position of the edge and the one i compute by say using the the parabola argument and therefore i can compensate for it now of course in a real imaging situation there'll be more than just the defocus of the lens so uh it probably uh doesn't pay to be too too careful about this okay now um [Music] suppose that you're a patent infringer or you're trying to infringe this pattern then there are a number of ways to go basically you look at the uh components of each claim and you see if there's one of the components that maybe you can avoid maybe you can do it in a different way or maybe you can do it in a better way just arguing whether something's better or not uh doesn't help it's it has to be different okay and now there's some some things here that aren't that uh pleasant one of them is this quantization of gradient directions and the reason it's not that great is because it introduces the awkwardness that the spacing is comes in two sizes you know pixel spacing and square root of two times pixel spacing and these these effects of defocus etc are they are on pixel grid they're in units of pixel spacing and so now we're sampling it in two different ways so we expect to get a slightly different error contribution so how can we avoid that so here are a couple of ways of so the idea is we don't want quantized gradient directions and so suppose our gradient is now if we follow the the preferred method in the patent again notice that that method doesn't show up in any of the claims it's so which is good because that means that you could easily circumvent the patent but it is the preferred method that's shown in in the specification then we would quantize this and we would work with the diagonal but let's instead say well how about this how about if i figured out what the if i knew what the value was there i know what the value is here you know how about this arrangement so this could be g0 and that's g plus and that's g minus then i would avoid the quantization of gradient direction and how do i do that well i only know the values on the actual pixel grid and so tada i interpolate so i use interpolation to go from the values on the grid to the value over here and because i'm interpolating along this line i can actually easily just use uh one 1d linear interpolation right so how does that work well if i have a function [Music] then i can say well i know the value here another value there and let's say it's a straight line in between and then the formula for this is i mean you can come up with that yourself but you can easily check it first of all it's linear in x and then at x equals a we get only this term and the b minus a cancels the b minus a and at x equals b uh this term drops out so we only have that and again the b minus a can okay so it's easy to do that in in 1d and of course you can extend that to so-called bilinear interpolation in in 2d but we don't actually even need that here okay so then we approximate the value here by interpolation and we can use more sophisticated methods of interpolation for example we can use cubic spline interpolation which we won't talk about here but it uh you know involves more points but gives you an interpolation that's a smooth cubic curve which in some cases is more accurate okay and then we perform whatever we did and we end up you know somewhere there and we we don't need the uh plane position step because we're actually on the gradient direction we're exactly where you want to be so um that you know you wonder why why didn't they do that well uh they don't say it but uh two two reasons one that i can think of one is that before we complained because sometimes we were using pixel spacing and sometimes square root of two times pixel spacing well now we can have anything in between not not just those two values and that correction graph the bias graph you know will be different for all of those values so you've got to i don't know quantize it build a lookup table no no not insurmountable problems but an extra hassle and you wonder you know is it worth it so that's one reason you might not choose this the other reason is that you did an interpolation and how accurate is that you know you've gone away from the values that you actually know for sure to something that you interpolated using a method that you chose in some arbitrary fashion you know how good is it so that that'll introduce some error as well now in many cases that error wouldn't matter it'd be so small but you know if you're going for a 40th of a pixel accuracy that even that small error can be significant okay so that's one one method uh one alternative and um one of its problems is uh changing size of this step uh from a pixel to square root of two well what if we get rid of that and we just have a fixed size so here's again our pixel grid we know the values at all of the intersections of these lines and this is g0 and now we say well let's just draw a circle and again we'll pick the we'll pick some gradient direction and now we use that value and that value and that value and those are equally spaced they're always the same distance from each other no matter what the gradient direction so i've kind of now dealt with two issues the one is there's no quantization of the gradient direction and there's none of this change of the you know what does s equals one mean uh before that could take on two values and over here it could take on a range of values well here it's always the same it's it's the the pixel spacing okay so of course i don't know the value here so i have to use interpolation which is now a 2d interpolation and so i can use bilinear or i could for example use bicubic and again you know why didn't they do this well it's extra work particularly with a bicubic we need to take into account more pixels than we have here you know go out to a five by five grid not just a three by three grid bilinear is not terribly accurate so we may not want to use that so anyway so you can see that even though a lot of things are described in the pattern there is also some hidden stuff going on they probably experimented with all of this and picked something that was very accurate and yet simple and that's how we got to that method now the last thing i want to talk about today is uh multiscale they do uh point out in this patent that you might want to do this at multiple scales they don't make a big fuss of it here but um um [Music] we already mentioned that that's something you want to do there are some edges which are very sharp and they will be best found at the highest resolution and then there's some other edges that are kind of blurry either because they're out of focus or because the object doesn't have a sharp transition you know if the object has a smoothly turning curvature then there'll be a transition in brightness from whatever the brightness is for one planar surface to the brightness for a different planar surface but they're all these in between values so the transition won't be on one pixel it'll be spread out over many pixels um see do i want to do this now or we haven't talked about chordic so um i'm just looking at time we don't have much time so when i why don't i do a codec so this was the method to go from e x e y to e zero e theta right where e zero e0 so it's just a cartesian to polar coordinate transformation and you know square roots and 810 take some computing you'd have to probably use a lookup table to speed that up so so uh what is their preferred uh implementation of this so the idea is that we rotate the coordinate system so we have we have we have some gradient e x e y and if we knew that angle then we could just rotate it down and then the length of it e x would be e0 and e y would be zero right so so we're aiming for where e zero is just e x prime and e theta e and e y prime is zero when you write it this way around that makes more sense okay now we don't know theta but we can rotate using some test angles and by an iterative method uh you know keep on improving keep on coming closer to to this situation where the y component of the gradient is reduced to zero so the superscripts in parenthesis again refer to the iteration number rather than being powers and of course we know in 2d a rotation matrix is just very simple okay and so one thing we can do is have a sequence of thetas that we try and when when we're finished with this iteration when we decide to stop the answer for the angle is just the sum of all of those increments we make so you can imagine various strategies like you can try an angle and suppose that this goes the wrong way well then maybe take half that angle and try that or you could try turning through that angle positive and through the angle negative and compare the two results see which one works better and then keep on reducing the size of the angle until you converge until or at least until you're you know happy enough with the results so so each time reduce each step reduce the magnitude of ey and and it will increase the other one right because as we rotate closer to the x-axis the projection of this onto the x-axis will will increase okay so that's the the iteration and how do we pick the thetas well we could do this theta zero is pi over two theta one is pi over four theta two is pi over eight you know and so on that'd be one sort of obvious approach so we try and turn through these different angles and see if it and we accept we accept the change if it satisfies this condition that magnitude that e y is reduced it doesn't have to be positive so it could overcompensate as long as the magnitude is reduced and you can always flip it to make it positive again if you want to um so you can do that um but it's expensive for two reasons one is you got the cosines and the sines but you could build a table obviously you could just store the cosines and sines of these angles the other one is you need each time you need four multiplications and two additions it doesn't sound like much but remember this happens several times because of the iterations and it happens at every pixel so it's expensive so how to avoid multiplication well what we can do is pick the angles very carefully suppose that we pick the angles so that there are inverse powers of 2. then that matrix there the rotation matrix just becomes this and the matrix part of this is very easy to compute you know multiplication by one doesn't cost us anything and this is just a shift and shifts are extremely cheap and then we have an addition so it reduces the whole thing down to two additions from uh four multiplications and two additions okay uh well then we have and we do this repeatedly as and you can see that the angle we're turning through gets smaller and smaller what what is that angle let's see we can compute it from that formula after a while it basically uh gets halved initially because of the non-linear nature of trigonometric functions it's not half but eventually it becomes r and then uh when you're done you you end up over here you end up with a product of all of these cosine theta i's and so what do you do with that well maybe you don't care because it's just a constant multiplier but suppose you do care then you can pre-compute it and actually it's 1.16 very quickly it converges very fast so you can just use 1.16 okay so that's the basic idea of chordic that we rotate but through special angles that have this property that were the you know the tangent of theta i is one over uh two to the i um and um well it takes a bit longer to explain it in the paper but that's the basic idea and so uh the iterative process is you um you change through that angle and if it improves the answer you keep it and then you keep track of whether it got negative or not and the first thing you have to do is get it to the first octant you know the whole idea here is that you're working in this regime but that's obviously trivial you just look at the signs of x and y and whether y is greater than x and you can reduce it to the first auction okay so next time we'll talk about multi-scale and we'll talk a little bit about sampling and aliasing because because that's part of the multi-scale story so again please start early on the quiz it's more work than the typical homework problem