{"id":1116,"date":"2014-08-12T12:37:23","date_gmt":"2014-08-12T16:37:23","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=1116"},"modified":"2019-11-01T11:17:16","modified_gmt":"2019-11-01T15:17:16","slug":"it-aint-easy-seeing-green-unless-you-have-matlab","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2014\/08\/12\/it-aint-easy-seeing-green-unless-you-have-matlab\/","title":{"rendered":"It Ain&#8217;t Easy Seeing Green (Unless You Have MATLAB)"},"content":{"rendered":"\r\n<div class=\"content\"><!--introduction--><p><i>I'd like to welcome guest blogger and ace MATLAB training content developer Matt Tearle for today's post. Thanks, Matt!<\/i><\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#e34c44e4-a8d1-44a0-ae6c-8c429e8a835e\">An Unusual Assignment<\/a><\/li><li><a href=\"#4ce9e85d-4007-4d77-83c0-aea317d34933\">MATLAB to the Rescue<\/a><\/li><li><a href=\"#35f1533b-a07a-4f05-9015-072538644796\">Real Life is Complicated<\/a><\/li><li><a href=\"#7734ebe3-6a5e-480f-8876-f68c6d6f4c72\">Developing a Green-Screen Algorithm<\/a><\/li><li><a href=\"#397c67b9-3168-4f87-8be8-ba56cfdf18ab\">Going Further. Because I Can.<\/a><\/li><li><a href=\"#3b7802d0-7639-46f7-94c1-dcc892387750\">The Moral of the Story<\/a><\/li><li><a href=\"#10a49a95-463a-45d1-ac32-dd4cad138b26\">Post Script: Fortune Favors the Brave (and\/or Crazy)<\/a><\/li><\/ul><\/div><h4>An Unusual Assignment<a name=\"e34c44e4-a8d1-44a0-ae6c-8c429e8a835e\"><\/a><\/h4><p>What do you do if you don't have an image processing background and your boss asks you to produce a music video for a MATLAB-vs-Simulink rap battle? It's probably not a question you've given much thought to. I admit that I hadn't either. But that's the awkward situation I found myself in a little while back.<\/p><p>A MATLAB-vs-Simulink rap battle -- obviously I <b>had<\/b> to be a part of this project! There was just one significant problem: I'm not a video producer. Yes, I have a nice digital camera. I have some very basic editing software. And I certainly have some... \"creative\" ideas. But I don't have any software designed to do serious video editing.<\/p><p>In particular, as soon as the concept was explained to me, I decided that we needed to do a \"green screen\" video:<\/p><div><ul><li>We'd record our rappers busting out their rhymes in front of a green backdrop<\/li><li>I'd make an entertaining background video<\/li><li>Then we'd put the two together, using the background video in place of the green screen:<\/li><\/ul><\/div><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greenscreen.png\" alt=\"\"> <\/p><p>[This background may not seem particularly \"entertaining\", but if I showed the real one we used, lawyers might get involved... And nobody wants that.]<\/p><h4>MATLAB to the Rescue<a name=\"4ce9e85d-4007-4d77-83c0-aea317d34933\"><\/a><\/h4><p>With no purpose-made software to achieve this, and only a few days before it was due (thanks, boss!), I turned to my default standby: MATLAB!<\/p><p>A video is just a sequence of images, so all I needed to do was to write MATLAB code to \"green-screen\" two images, then read in two videos and apply my algorithm frame-by-frame. A quick hunt around in the documentation, and I had the framework code ready to go:<\/p><pre class=\"codeinput\"><span class=\"comment\">% Open the input video files<\/span>\r\nv1 = VideoReader(<span class=\"string\">'background.mp4'<\/span>);\r\nv2 = VideoReader(<span class=\"string\">'foreground.mp4'<\/span>);\r\n<span class=\"comment\">% Determine the number of frames for the final video<\/span>\r\nnFrames = min(v1.NumberOfFrames,v2.NumberOfFrames);\r\n<span class=\"comment\">% Open a video file to write the result to<\/span>\r\nvOut = VideoWriter(<span class=\"string\">'mixedvideo.mp4'<\/span>,<span class=\"string\">'MPEG-4'<\/span>);\r\nvOut.FrameRate = 24;\r\nopen(vOut);\r\n<span class=\"comment\">% Loop over all the frames<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:nFrames\r\n    <span class=\"comment\">% Get the kth frame of both inputs<\/span>\r\n    x = read(v1,k);\r\n    y = read(v2,k);\r\n    <span class=\"comment\">% Mix them together<\/span>\r\n    z = magicHappensHere(x,y);  <span class=\"comment\">% TODO: Write this...<\/span>\r\n    <span class=\"comment\">% Write the result to file<\/span>\r\n    writeVideo(vOut,z)\r\n<span class=\"keyword\">end<\/span>\r\n<span class=\"comment\">% And we're done!<\/span>\r\nclose(vOut);\r\n<\/pre><p>Now all I needed to do was write the magic function that would do the greenscreening...<\/p><h4>Real Life is Complicated<a name=\"35f1533b-a07a-4f05-9015-072538644796\"><\/a><\/h4><p>Actually, I soon discovered that I needed to do a bit more preprocessing first. I had created the background video on my computer, but filmed the foreground on my digital camera. Of course the dimensions didn't match... No worries! That's what Image Processing Toolbox is for. I added a line of code to define the desired dimensions of the output video. Then <tt>imresize<\/tt> took care of the rest:<\/p><pre class=\"codeinput\"><span class=\"comment\">% Loop over all the frames<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:nFrames\r\n    <span class=\"comment\">% Get the kth frame of both inputs<\/span>\r\n    x = imresize(read(v1,k),outDims);\r\n    y = imresize(read(v2,k),outDims);\r\n    <span class=\"comment\">% Mix them together<\/span>\r\n    z = magicHappensHere(x,y);\r\n    <span class=\"comment\">% Write the result to file<\/span>\r\n    writeVideo(vOut,z)\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>Wow, this image processing stuff is <b>easy<\/b>!  What's next, then?<\/p><h4>Developing a Green-Screen Algorithm<a name=\"7734ebe3-6a5e-480f-8876-f68c6d6f4c72\"><\/a><\/h4><p>The first thing to try, of course, when developing an algorithm is to see if one already exists. Sadly, searching the Image Processing Toolbox documentation didn't turn up a <tt>greenscreen<\/tt> function or anything I could see that would do the job. So it looks like I have to do this from scratch.<\/p><p>If I can identify the green pixels, then I can just replace those pixels in the foreground image with the corresponding pixels in the background image. So all of this really just boils down to: how do I find the green pixels? I know that a true-color image in MATLAB is an m-by-n-by-3 array, where each of the three m-by-n \"planes\" represent red, green, and blue intensity, respectively:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/rgb.png\" alt=\"\"> <\/p><p>\"Green\" would, ideally, be [0, 255, 0] (in <tt>uint8<\/tt>). So a simple way to find the green pixels would be to do some kind of thresholding:<\/p><pre class=\"codeinput\">isgreen = (y(:,:,1) &lt; 175) &amp; (y(:,:,2) &gt; 200) &amp; (y(:,:,3) &lt; 175);\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/basicthreshold.png\" alt=\"\"> <\/p><p>This didn't work too badly, but neither was it great. One of the problems is the use of magic threshold values -- three of them, no less! Playing around with three magic parameters to try to find the right combination doesn't sound like fun to me. And besides, using absolute values like this might not be the best way to judge \"greenness\". If you look at the example image above, you can see that there's a fair amount of variation in the green background. (This is mainly because we were using an unevenly lit conference room that happened to have a green wall -- this was clearly a very professional production!)<\/p><p>At this point, I had to ask myself <a href=\"https:\/\/blogs.mathworks.com\/steve\/2010\/12\/17\/what-color-is-green\/\">\"what does it really mean to be 'green'?\"<\/a>. As a human, I can easily distinguish between the various shades of green on the wall, despite the difference in brightness and shading. At first, this lead me into <a href=\"https:\/\/blogs.mathworks.com\/steve\/2010\/12\/23\/two-dimensional-histograms\/\">thinking about color spaces<\/a> and how a color is a point in some 3-D space (RGB, CYM, HSV, etc.). Perhaps I could define \"green\" as \"sufficiently close to [0 1 0] in RGB space\". And so followed a few experiments with rounding, distance calculations, even playing with <tt>rgb2ind<\/tt>. But to no great success:<\/p><pre class=\"codeinput\">yd = round(double(y)\/255);\r\nisgreen = (yd(:,:,1)==0) &amp; (yd(:,:,2)==1) &amp; (yd(:,:,3)==0);\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/8bitrounding.png\" alt=\"\"> <\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/8bitthresholding.png\" alt=\"\"> <\/p><p>As this image shows, the problem with our makeshift green screen was that it wasn't particularly green -- more a yellow-green, with significant variations in brightness.<\/p><p>Just as I was started to get worried that this wasn't going to work at all, I realized that one obvious characteristic of the green pixels is that the green value is significantly <b>greater<\/b> than the blue and red values. This is, of course, embedded in my simple thresholding code above, but I can do away with at least some of the magic constants by looking at <b>relative<\/b> values. So how about calculating a \"greenness\" value from the three color channels:<\/p><pre class=\"codeinput\">yd = double(y)\/255;\r\n<span class=\"comment\">% Greenness = G*(G-R)*(G-B)<\/span>\r\ngreenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greenness.png\" alt=\"\"> <\/p><p>Now <b>that<\/b> looks very promising. If I threshold that, I should get what I'm looking for. But rather than use an absolute threshold value maybe I can use something based on the distribution of the greenness values. In particular, I can exploit one nice feature of the images: they have a lot of green pixels. The histogram of greenness values looks like this:<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/histogram.png\" alt=\"\"> <\/p><p>I tried a simple thresholding based on the average greenness value, ignoring the negative values:<\/p><pre class=\"codeinput\">thresh = 0.3*mean(greenness(greenness&gt;0));\r\nisgreen = greenness &gt; thresh;\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greennessthreshold.png\" alt=\"\"> <\/p><p>Success! OK, so there's one magic constant left in there, but it was pretty quick and easy to tune, and seemed to be quite robust to minor changes. I'm calling that a good result. Now that I know where the green pixels are, I just have to replace them with the corresponding pixels from the background.<\/p><pre class=\"codeinput\"><span class=\"comment\">% Start with the foreground (essentially preallocation)<\/span>\r\nz = y;\r\n<span class=\"comment\">% Find the green pixels<\/span>\r\nyd = double(y)\/255;\r\n<span class=\"comment\">% Greenness = G*(G-R)*(G-B)<\/span>\r\ngreenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\nthresh = 0.3*mean(greenness(greenness&gt;0));\r\nisgreen = greenness &gt; thresh;\r\n<span class=\"comment\">% Blend the images<\/span>\r\n<span class=\"comment\">% Loop over the 3 color planes (RGB)<\/span>\r\n<span class=\"keyword\">for<\/span> j = 1:3\r\n    rgb1 = x(:,:,j);  <span class=\"comment\">% Extract the jth plane of the background<\/span>\r\n    rgb2 = y(:,:,j);  <span class=\"comment\">% Extract the jth plane of the foreground<\/span>\r\n    <span class=\"comment\">% Replace the green pixels of the foreground with the background<\/span>\r\n    rgb2(isgreen) = rgb1(isgreen);\r\n    <span class=\"comment\">% Put the combined image into the output<\/span>\r\n    z(:,:,j) = rgb2;\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><p>I love using logical indexing -- it's probably my favorite construct in the MATLAB language. Here I'm using it to replace the green pixels of the foreground (<tt>rgb2(isgreen)<\/tt>) with the corresponding background pixels (<tt>rgb1(isgreen)<\/tt>).<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/result.png\" alt=\"\"> <\/p><h4>Going Further. Because I Can.<a name=\"397c67b9-3168-4f87-8be8-ba56cfdf18ab\"><\/a><\/h4><p>At this point I had a video that would have sufficed, but there was a noticeable green outline around my rappers. (It's not necessarily apparent in a single still frame.) It was only Sunday morning by this point, which meant that I still had a few hours left to tinker. Could I find a way to \"shave\" a pixel or two from the edge of the black silhouette in the above image of <tt>isgreen<\/tt>? To do that, I'd first need to find that edge. I don't know much about image processing, but I do know that \"edge detection\" is a thing that image people do, so it's time for me to search the doc... Sure enough, there's an <tt>edge<\/tt> function:<\/p><pre class=\"codeinput\">outline = edge(isgreen,<span class=\"string\">'roberts'<\/span>);\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/edge.png\" alt=\"\"> <\/p><p>Easy! Now if I can thicken that edge and combine it with the original <tt>isgreen<\/tt>, I'll have the slightly larger greenscreen mask that I need. Two more Image Processing Toolbox functions help me thicken the edge:<\/p><pre class=\"codeinput\">se = strel(<span class=\"string\">'disk'<\/span>,1);\r\noutline = imdilate(outline,se);\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/thickedge.png\" alt=\"\"> <\/p><p>The edge variable <tt>outline<\/tt> is a logical array, so I can easily combine it with <tt>isgreen<\/tt>, to get my final result:<\/p><pre class=\"codeinput\">isgreen = isgreen | outline;\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/finalthreshold.png\" alt=\"\"> <\/p><p>Put it all together and I have a pretty simple script:<\/p><pre class=\"codeinput\"><span class=\"comment\">% Open the input video files<\/span>\r\nv1 = VideoReader(<span class=\"string\">'background.mp4'<\/span>);\r\nv2 = VideoReader(<span class=\"string\">'foreground.mp4'<\/span>);\r\n<span class=\"comment\">% Determine the number of frames for the final video<\/span>\r\nnFrames = min(v1.NumberOfFrames,v2.NumberOfFrames);\r\n<span class=\"comment\">% Set the output dimensions<\/span>\r\noutDims = [400 640];\r\n<span class=\"comment\">% Open a video file to write the result to<\/span>\r\nvOut = VideoWriter(<span class=\"string\">'mixedvideo.mp4'<\/span>,<span class=\"string\">'MPEG-4'<\/span>);\r\nvOut.FrameRate = 24;\r\nopen(vOut);\r\n<span class=\"comment\">% Loop over all the frames<\/span>\r\n<span class=\"keyword\">for<\/span> k = 1:nFrames\r\n    <span class=\"comment\">% Get the kth frame of both inputs<\/span>\r\n    x = imresize(read(v1,k),outDims);\r\n    y = imresize(read(v2,k),outDims);\r\n    <span class=\"comment\">% Mix them together<\/span>\r\n    z = y;  <span class=\"comment\">% Preallocate space for the result<\/span>\r\n    <span class=\"comment\">% Find the green pixels in the foreground (y)<\/span>\r\n    yd = double(y)\/255;\r\n    <span class=\"comment\">% Greenness = G*(G-R)*(G-B)<\/span>\r\n    greenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\n    <span class=\"comment\">% Threshold the greenness value<\/span>\r\n    thresh = 0.3*mean(greenness(greenness&gt;0));\r\n    isgreen = greenness &gt; thresh;\r\n    <span class=\"comment\">% Thicken the outline to expand the greenscreen mask a little<\/span>\r\n    outline = edge(isgreen,<span class=\"string\">'roberts'<\/span>);\r\n    se = strel(<span class=\"string\">'disk'<\/span>,1);\r\n    outline = imdilate(outline,se);\r\n    isgreen = isgreen | outline;\r\n    <span class=\"comment\">% Blend the images<\/span>\r\n    <span class=\"comment\">% Loop over the 3 color planes (RGB)<\/span>\r\n    <span class=\"keyword\">for<\/span> j = 1:3\r\n        rgb1 = x(:,:,j);  <span class=\"comment\">% Extract the jth plane of the background<\/span>\r\n        rgb2 = y(:,:,j);  <span class=\"comment\">% Extract the jth plane of the foreground<\/span>\r\n        <span class=\"comment\">% Replace the green pixels of the foreground with the background<\/span>\r\n        rgb2(isgreen) = rgb1(isgreen);\r\n        <span class=\"comment\">% Put the combined image into the output<\/span>\r\n        z(:,:,j) = rgb2;\r\n    <span class=\"keyword\">end<\/span>\r\n    <span class=\"comment\">% Write the result to file<\/span>\r\n    writeVideo(vOut,z)\r\n<span class=\"keyword\">end<\/span>\r\n<span class=\"comment\">% And we're done!<\/span>\r\nclose(vOut);\r\n<\/pre><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/final.png\" alt=\"\"> <\/p><h4>The Moral of the Story<a name=\"3b7802d0-7639-46f7-94c1-dcc892387750\"><\/a><\/h4><p>This is not the best greenscreen code ever written (although it <b>is<\/b> the best greenscreen code <b>I've<\/b> ever written). But the real point is that it shouldn't even exist -- I had no prior image processing knowledge, I was working with videos of different dimensions, the greenscreen was unevenly lit (and not even particularly green), and I had a weekend to make it happen. The fact that I was able to do this is one of the main reasons I love MATLAB (and, of course, the associated toolboxes): I was able to spend my time on the heart of <b>my algorithm<\/b>, not all the coding details. Different-sized images? Fixed with one function call (<tt>imresize<\/tt>). Need to find edges of a region? One function call (<tt>edge<\/tt>). Want to thicken that edge? Two function calls (<tt>strel<\/tt> and <tt>imdilate<\/tt>). The final script is less than 50 lines.<\/p><p>My time was spent wrestling with the core of the problem: how I was going to figure out what \"green\" looked like in my images. Because I could spend my time there, I was able to tinker with different ideas. This, to me, is what MATLAB is all about: rapid prototyping -- ideas becoming working code.<\/p><h4>Post Script: Fortune Favors the Brave (and\/or Crazy)<a name=\"10a49a95-463a-45d1-ac32-dd4cad138b26\"><\/a><\/h4><p>Having made the video and feeling pretty good about what I managed to pull together in a few hours of MATLAB, I later realized that I got lucky. The way I calculated greenness (<tt>G*(G-R)*(G-B)<\/tt>) was intended to give negative results for anything that didn't have more green than red or blue. But I simply overlooked the possibility of getting a high positive value if <b>both<\/b> <tt>G-R<\/tt> and <tt>G-B<\/tt> were negative. According to my formula, a light magenta -- e.g. <tt>[1 0.6 1]<\/tt> -- would be very green! Luckily for me, my source videos didn't include anything like that. A more robust approach would have been something like this:<\/p><pre class=\"codeinput\">yd = double(y)\/255;\r\n<span class=\"comment\">% Greenness = G*(G-R)*(G-B)<\/span>\r\ngreenness = yd(:,:,2).*max(yd(:,:,2)-yd(:,:,1),0).*max(yd(:,:,2)-yd(:,:,3),0);\r\n<\/pre><p>Now any pixel with more red or more blue than green will have a greenness of 0. Everything else works the same.<\/p><p>I think there's a nice moral in this postscript, as well. I love using MATLAB to solve a problem -- a real, immediate problem that I need to solve right now. I don't necessarily need the best solution, or the solution to the most general problem. The greenscreen algorithm I hacked together wasn't the best, but it worked for my application. If my image had light magenta in it, I would have discovered this bug and fixed it, but that didn't happen and didn't matter.<\/p><p>If I had approached this assignment like a professional software developer, I'd still be working on it. But as a MATLAB hacker, I was done in a weekend!<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_b2e543e264fb4f7899b3edb10445d7be() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='b2e543e264fb4f7899b3edb10445d7be ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' b2e543e264fb4f7899b3edb10445d7be';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2014 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_b2e543e264fb4f7899b3edb10445d7be()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014a<br><\/p><\/div><!--\r\nb2e543e264fb4f7899b3edb10445d7be ##### SOURCE BEGIN #####\r\n%% It Ain't Easy Seeing Green (Unless You Have MATLAB)\r\n%\r\n% _I'd like to welcome guest blogger and ace MATLAB training content\r\n% developer Matt Tearle for today's post. Thanks, Matt!_\r\n%\r\n%% An Unusual Assignment\r\n% What do you do if you don't have an image processing background and your\r\n% boss asks you to produce a music video for a MATLAB-vs-Simulink rap\r\n% battle? It's probably not a question you've given much thought to. I\r\n% admit that I hadn't either. But that's the awkward situation I found\r\n% myself in a little while back.\r\n%\r\n% A MATLAB-vs-Simulink rap battle REPLACE_WITH_DASH_DASH obviously I *had* to be a part of this\r\n% project! There was just one significant problem: I'm not a video\r\n% producer. Yes, I have a nice digital camera. I have some very basic\r\n% editing software. And I certainly have some... \"creative\" ideas. But I\r\n% don't have any software designed to do serious video editing.\r\n%\r\n% In particular, as soon as the concept was explained to me, I decided that\r\n% we needed to do a \"green screen\" video:\r\n%\r\n% * We'd record our rappers busting out their rhymes in front of a green\r\n% backdrop\r\n% * I'd make an entertaining background video\r\n% * Then we'd put the two together, using the background video in place of\r\n% the green screen:\r\n% \r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greenscreen.png>>\r\n%\r\n% [This background may not seem particularly \"entertaining\", but if I\r\n% showed the real one we used, lawyers might get involved... And nobody\r\n% wants that.]\r\n% \r\n%% MATLAB to the Rescue\r\n% With no purpose-made software to achieve this, and only a few days before\r\n% it was due (thanks, boss!), I turned to my default standby: MATLAB!\r\n%\r\n% A video is just a\r\n% sequence of images, so all I needed to do was to write MATLAB code to\r\n% \"green-screen\" two images, then read in two videos and apply my algorithm\r\n% frame-by-frame. A quick hunt around in the documentation, and I had the\r\n% framework code ready to go:\r\n\r\n%%\r\n\r\n% Open the input video files\r\nv1 = VideoReader('background.mp4');\r\nv2 = VideoReader('foreground.mp4');\r\n% Determine the number of frames for the final video\r\nnFrames = min(v1.NumberOfFrames,v2.NumberOfFrames);\r\n% Open a video file to write the result to\r\nvOut = VideoWriter('mixedvideo.mp4','MPEG-4');\r\nvOut.FrameRate = 24;\r\nopen(vOut);\r\n% Loop over all the frames\r\nfor k = 1:nFrames\r\n    % Get the kth frame of both inputs\r\n    x = read(v1,k);\r\n    y = read(v2,k);\r\n    % Mix them together\r\n    z = magicHappensHere(x,y);  % TODO: Write this...\r\n    % Write the result to file\r\n    writeVideo(vOut,z)\r\nend\r\n% And we're done!\r\nclose(vOut);\r\n\r\n%%\r\n% Now all I needed to do was write the magic function that would do the\r\n% greenscreening...\r\n\r\n%% Real Life is Complicated\r\n% Actually, I soon discovered that I needed to do a bit more preprocessing\r\n% first. I had created the background video on my computer, but filmed the\r\n% foreground on my digital camera. Of course the dimensions didn't match...\r\n% No worries! That's what Image Processing Toolbox is for. I added a line\r\n% of code to define the desired dimensions of the output video. Then\r\n% |imresize| took care of the rest:\r\n\r\n%%\r\n\r\n% Loop over all the frames\r\nfor k = 1:nFrames\r\n    % Get the kth frame of both inputs\r\n    x = imresize(read(v1,k),outDims);\r\n    y = imresize(read(v2,k),outDims);\r\n    % Mix them together\r\n    z = magicHappensHere(x,y);\r\n    % Write the result to file\r\n    writeVideo(vOut,z)\r\nend\r\n\r\n%%\r\n% Wow, this image processing stuff is *easy*!  What's next, then?\r\n\r\n%% Developing a Green-Screen Algorithm\r\n% The first thing to try, of course, when developing an algorithm is to see\r\n% if one already exists. Sadly, searching the Image Processing Toolbox\r\n% documentation didn't turn up a |greenscreen| function or anything I could\r\n% see that would do the job. So it looks like I have to do this from\r\n% scratch.\r\n%\r\n% If I can identify the green pixels, then I can just replace those pixels\r\n% in the foreground image with the corresponding pixels in the background\r\n% image. So all of this really just boils down to: how do I find the green\r\n% pixels? I know that a true-color image in MATLAB is an m-by-n-by-3\r\n% array, where each of the three m-by-n \"planes\" represent red, green, and\r\n% blue intensity, respectively:\r\n%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/rgb.png>>\r\n%\r\n% \"Green\" would, ideally, be [0, 255, 0] (in |uint8|). So a simple way to\r\n% find the green pixels would be to do some kind of thresholding:\r\n% \r\n%%\r\nisgreen = (y(:,:,1) < 175) & (y(:,:,2) > 200) & (y(:,:,3) < 175);\r\n\r\n%%\r\n%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/basicthreshold.png>>\r\n%\r\n% This didn't work too badly, but neither was it great. One of the\r\n% problems is the use of magic threshold values REPLACE_WITH_DASH_DASH three of them, no less!\r\n% Playing around with three magic parameters to try to find the right\r\n% combination doesn't sound like fun to me. And besides, using absolute\r\n% values like this might not be the best way to judge \"greenness\". If you\r\n% look at the example image above, you can see that there's a fair amount\r\n% of variation in the green background. (This is mainly because we were\r\n% using an unevenly lit conference room that happened to have a green wall\r\n% REPLACE_WITH_DASH_DASH this was clearly a very professional production!)\r\n%\r\n% At this point, I had to ask myself\r\n% <https:\/\/blogs.mathworks.com\/steve\/2010\/12\/17\/what-color-is-green\/ \"what\r\n% does it really mean to be 'green'?\">. As a human, I can easily\r\n% distinguish between the various shades of green on the wall, despite the\r\n% difference in brightness and shading. At first, this lead me into\r\n% <https:\/\/blogs.mathworks.com\/steve\/2010\/12\/23\/two-dimensional-histograms\/\r\n% thinking about color spaces> and how a color is a point in some 3-D space\r\n% (RGB, CYM, HSV, etc.). Perhaps I could define \"green\" as \"sufficiently\r\n% close to [0 1 0] in RGB space\". And so followed a few experiments with\r\n% rounding, distance calculations, even playing with |rgb2ind|. But to no\r\n% great success:\r\n\r\nyd = round(double(y)\/255);\r\nisgreen = (yd(:,:,1)==0) & (yd(:,:,2)==1) & (yd(:,:,3)==0);\r\n%%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/8bitrounding.png>>\r\n%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/8bitthresholding.png>>\r\n% \r\n%\r\n% As this image shows, the problem with our makeshift green screen was that\r\n% it wasn't particularly green REPLACE_WITH_DASH_DASH more a yellow-green, with significant\r\n% variations in brightness.\r\n%\r\n% Just as I was started to get worried that this wasn't going to work at\r\n% all, I realized that one obvious characteristic of the green pixels is\r\n% that the green value is significantly *greater* than the blue and red\r\n% values. This is, of course, embedded in my simple thresholding code\r\n% above, but I can do away with at least some of the magic constants by\r\n% looking at *relative* values. So how about calculating a \"greenness\"\r\n% value from the three color channels:\r\n\r\nyd = double(y)\/255;\r\n% Greenness = G*(G-R)*(G-B)\r\ngreenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\n\r\n%%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greenness.png>>\r\n%\r\n% Now *that* looks very promising. If I threshold that, I should get what\r\n% I'm looking for. But rather than use an absolute threshold value maybe I\r\n% can use something based on the distribution of the greenness values. In\r\n% particular, I can exploit one nice feature of the images: they have a lot\r\n% of green pixels. The histogram of greenness values looks like this:\r\n%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/histogram.png>>\r\n%\r\n% I tried a simple thresholding based on the average greenness value,\r\n% ignoring the negative values:\r\n\r\nthresh = 0.3*mean(greenness(greenness>0));\r\nisgreen = greenness > thresh;\r\n%%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/greennessthreshold.png>>\r\n%\r\n% Success! OK, so there's one magic constant left in there, but it was\r\n% pretty quick and easy to tune, and seemed to be quite robust to minor\r\n% changes. I'm calling that a good result. Now that I know where the green\r\n% pixels are, I just have to replace them with the corresponding pixels\r\n% from the background.\r\n%%\r\n\r\n% Start with the foreground (essentially preallocation)\r\nz = y;\r\n% Find the green pixels\r\nyd = double(y)\/255;\r\n% Greenness = G*(G-R)*(G-B)\r\ngreenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\nthresh = 0.3*mean(greenness(greenness>0));\r\nisgreen = greenness > thresh;\r\n% Blend the images\r\n% Loop over the 3 color planes (RGB)\r\nfor j = 1:3\r\n    rgb1 = x(:,:,j);  % Extract the jth plane of the background\r\n    rgb2 = y(:,:,j);  % Extract the jth plane of the foreground\r\n    % Replace the green pixels of the foreground with the background\r\n    rgb2(isgreen) = rgb1(isgreen);\r\n    % Put the combined image into the output\r\n    z(:,:,j) = rgb2;\r\nend\r\n%%\r\n% I love using logical indexing REPLACE_WITH_DASH_DASH it's probably my favorite construct in\r\n% the MATLAB language. Here I'm using it to replace the green pixels of\r\n% the foreground (|rgb2(isgreen)|) with the corresponding background pixels\r\n% (|rgb1(isgreen)|).\r\n%\r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/result.png>>\r\n%\r\n\r\n%% Going Further. Because I Can.\r\n% At this point I had a video that would have sufficed, but there was a\r\n% noticeable green outline around my rappers. (It's not necessarily\r\n% apparent in a single still frame.) It was only Sunday morning by this\r\n% point, which meant that I still had a few hours left to tinker. Could I\r\n% find a way to \"shave\" a pixel or two from the edge of the black\r\n% silhouette in the above image of |isgreen|? To do that, I'd first need\r\n% to find that edge. I don't know much about image processing, but I do\r\n% know that \"edge detection\" is a thing that image people do, so it's time\r\n% for me to search the doc... Sure enough, there's an |edge| function:\r\n\r\n%%\r\noutline = edge(isgreen,'roberts');\r\n%%\r\n% \r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/edge.png>>\r\n% \r\n% Easy! Now if I can thicken that edge and combine it with the original\r\n% |isgreen|, I'll have the slightly larger greenscreen mask that I need.\r\n% Two more Image Processing Toolbox functions help me thicken the edge:\r\n%%\r\nse = strel('disk',1);\r\noutline = imdilate(outline,se);\r\n%%\r\n% \r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/thickedge.png>>\r\n% \r\n% The edge variable |outline| is a logical array, so I can easily combine\r\n% it with |isgreen|, to get my final result:\r\n%%\r\nisgreen = isgreen | outline;\r\n%%\r\n% \r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/finalthreshold.png>>\r\n% \r\n% Put it all together and I have a pretty simple script:\r\n%%\r\n\r\n% Open the input video files\r\nv1 = VideoReader('background.mp4');\r\nv2 = VideoReader('foreground.mp4');\r\n% Determine the number of frames for the final video\r\nnFrames = min(v1.NumberOfFrames,v2.NumberOfFrames);\r\n% Set the output dimensions\r\noutDims = [400 640];\r\n% Open a video file to write the result to\r\nvOut = VideoWriter('mixedvideo.mp4','MPEG-4');\r\nvOut.FrameRate = 24;\r\nopen(vOut);\r\n% Loop over all the frames\r\nfor k = 1:nFrames\r\n    % Get the kth frame of both inputs\r\n    x = imresize(read(v1,k),outDims);\r\n    y = imresize(read(v2,k),outDims);\r\n    % Mix them together\r\n    z = y;  % Preallocate space for the result\r\n    % Find the green pixels in the foreground (y)\r\n    yd = double(y)\/255;\r\n    % Greenness = G*(G-R)*(G-B)\r\n    greenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));\r\n    % Threshold the greenness value\r\n    thresh = 0.3*mean(greenness(greenness>0));\r\n    isgreen = greenness > thresh;\r\n    % Thicken the outline to expand the greenscreen mask a little\r\n    outline = edge(isgreen,'roberts');\r\n    se = strel('disk',1);\r\n    outline = imdilate(outline,se);\r\n    isgreen = isgreen | outline;\r\n    % Blend the images\r\n    % Loop over the 3 color planes (RGB)\r\n    for j = 1:3\r\n        rgb1 = x(:,:,j);  % Extract the jth plane of the background\r\n        rgb2 = y(:,:,j);  % Extract the jth plane of the foreground\r\n        % Replace the green pixels of the foreground with the background\r\n        rgb2(isgreen) = rgb1(isgreen);\r\n        % Put the combined image into the output\r\n        z(:,:,j) = rgb2;\r\n    end\r\n    % Write the result to file\r\n    writeVideo(vOut,z)\r\nend\r\n% And we're done!\r\nclose(vOut);\r\n%%\r\n% \r\n% <<https:\/\/blogs.mathworks.com\/images\/steve\/2014\/final.png>>\r\n% \r\n\r\n%% The Moral of the Story\r\n% This is not the best greenscreen code ever written (although it *is* the\r\n% best greenscreen code *I've* ever written). But the real point is that\r\n% it shouldn't even exist REPLACE_WITH_DASH_DASH I had no prior image processing knowledge, I\r\n% was working with videos of different dimensions, the greenscreen was\r\n% unevenly lit (and not even particularly green), and I had a weekend to\r\n% make it happen. The fact that I was able to do this is one of the main\r\n% reasons I love MATLAB (and, of course, the associated toolboxes): I was\r\n% able to spend my time on the heart of *my algorithm*, not all the coding\r\n% details. Different-sized images? Fixed with one function call\r\n% (|imresize|). Need to find edges of a region? One function call\r\n% (|edge|). Want to thicken that edge? Two function calls (|strel| and\r\n% |imdilate|). The final script is less than 50 lines.\r\n%\r\n% My time was spent wrestling with the core of the problem: how I was going\r\n% to figure out what \"green\" looked like in my images. Because I could\r\n% spend my time there, I was able to tinker with different ideas. This, to\r\n% me, is what MATLAB is all about: rapid prototyping REPLACE_WITH_DASH_DASH ideas becoming\r\n% working code.\r\n\r\n%% Post Script: Fortune Favors the Brave (and\/or Crazy)\r\n% Having made the video and feeling pretty good about what I managed to\r\n% pull together in a few hours of MATLAB, I later realized that I got\r\n% lucky. The way I calculated greenness (|G*(G-R)*(G-B)|) was intended to\r\n% give negative results for anything that didn't have more green than red\r\n% or blue. But I simply overlooked the possibility of getting a high\r\n% positive value if *both* |G-R| and |G-B| were negative. According to my\r\n% formula, a light magenta REPLACE_WITH_DASH_DASH e.g. |[1 0.6 1]| REPLACE_WITH_DASH_DASH would be very green!\r\n% Luckily for me, my source videos didn't include anything like that. A\r\n% more robust approach would have been something like this:\r\n%%\r\nyd = double(y)\/255;\r\n% Greenness = G*(G-R)*(G-B)\r\ngreenness = yd(:,:,2).*max(yd(:,:,2)-yd(:,:,1),0).*max(yd(:,:,2)-yd(:,:,3),0);\r\n%%\r\n% Now any pixel with more red or more blue than green will have a greenness\r\n% of 0. Everything else works the same.\r\n%\r\n% I think there's a nice moral in this postscript, as well. I love using\r\n% MATLAB to solve a problem REPLACE_WITH_DASH_DASH a real, immediate problem that I need to\r\n% solve right now. I don't necessarily need the best solution, or the\r\n% solution to the most general problem. The greenscreen algorithm I hacked\r\n% together wasn't the best, but it worked for my application. If my image\r\n% had light magenta in it, I would have discovered this bug and fixed it,\r\n% but that didn't happen and didn't matter.\r\n%\r\n% If I had approached this assignment like a professional software\r\n% developer, I'd still be working on it. But as a MATLAB hacker, I was\r\n% done in a weekend!\r\n\r\n##### SOURCE END ##### b2e543e264fb4f7899b3edb10445d7be\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/2014\/final.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p><i>I'd like to welcome guest blogger and ace MATLAB training content developer Matt Tearle for today's post. Thanks, Matt!<\/i>... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/steve\/2014\/08\/12\/it-aint-easy-seeing-green-unless-you-have-matlab\/\">read more >><\/a><\/p>","protected":false},"author":42,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[278,402,228,124,156,122,308,380,1089,753,188,106,1087,719,1091],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/1116"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/users\/42"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/comments?post=1116"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/1116\/revisions"}],"predecessor-version":[{"id":1117,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/1116\/revisions\/1117"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/media?parent=1116"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/categories?post=1116"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/tags?post=1116"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}