{"id":455,"date":"2012-06-16T21:58:33","date_gmt":"2012-06-16T20:58:33","guid":{"rendered":"http:\/\/blog.blockos.org\/?p=455"},"modified":"2016-05-05T18:19:52","modified_gmt":"2016-05-05T17:19:52","slug":"hole-in-the-sky","status":"publish","type":"post","link":"https:\/\/blog.blockos.org\/?p=455","title":{"rendered":"Hole in the sky."},"content":{"rendered":"<p>I first came across the \u00e0 trous algorithm when I was looking for a good edge avoiding filter for <a href=\"http:\/\/en.wikipedia.org\/wiki\/Screen_Space_Ambient_Occlusion\" target=\"_blank\">SSAO<\/a>. Back then I was trying to code a demo. I had depth of field, some kind of hdr, spherical harmonics ambient lighting and <a href=\"http:\/\/sirkan.iit.bme.hu\/~szirmay\/ambientlight_link.htm\" target=\"_blank\">screen space ambient transfer<\/a>. I tried several techniques like the bilateral filter described in this <a href=\"http:\/\/www.ppsloan.org\/publications\/ProxyPG.pdf\" target=\"_blank\">paper<\/a> (4.2 Bilateral Upsampling). I don&#8217;t remember why but I have issues with this filter so I went on a quest for a better filter. Long story short, thanks to <a href=\"http:\/\/kesen.realtimerendering.com\/\" target=\"_blank\">Ke-Sen Huang&#8217;s page<\/a> I read a paper called Edge-Avoiding A-Trous Wavelet Transform for fast Global Illumination Filtering (<a href=\"http:\/\/www.uni-ulm.de\/fileadmin\/website_uni_ulm\/iui.inst.100\/institut\/Papers\/atrousGIfilter.pdf\" target=\"_blank\">paper<\/a>, <a href=\"http:\/\/www.highperformancegraphics.org\/previous\/www_2010\/media\/RayTracing_I\/HPG2010_RayTracing_I_Dammertz.pdf\" target=\"_blank\">slides<\/a>).<br \/>\nI won&#8217;t make a lecture about wavelet transforms but this is how it roughly works.<br \/>\nWe start with the source image <strong>I<\/strong>. At each iteration <strong>i<\/strong>we will compute 2 images, \\(d_i\\) and \\(c_i\\). The first one will hold the filtered output and the second the details. This can be translated as:<\/p>\n<ul>\n<li>\\(c_0 = I\\)<\/li>\n<li>\\(c_{i+1}(p) = \\sum_{k \\in H} h_i(k) c_i(p+2^ik)\\)<\/li>\n<li>\\(d_i = c_i &#8211; c_{i+1}\\)<\/li>\n<\/ul>\n<p>With \\(h_i\\) being the filter. And \\(H\\) the size of our filter. When the desired iteration \\(n\\) is reached, the set of detail images and the last filtered image form our wavelet transform.<\/p>\n<ul>\n<li>\\(w(I) = \\{d_0, d1, \\cdots, d_{n-1}, c_n\\}\\)<\/li>\n<\/ul>\n<p>As you can see, the space between each sample doubles at each iteration. Hence the name &#8220;\u00e0-trous&#8221; (with holes). Moreover the dimensions of the filtered images do not change. As opposed to Haar wavelet transform where \\(c_{i+1}\\) is half the size of \\(c_i\\).<br \/>\nThe next step is called the synthesis. We simply add the elements of the \\(w(I)\\) to produce the final image.<\/p>\n<ul>\n<li>\\(I&#8217; = c_n + \\sum_{0 \\leq i &lt; n}d_i\\)<\/li>\n<\/ul>\n<p>The wavelet filter can be seen as a pyramidal algorithm. The wavelet transform is the &#8220;pull&#8221; phase. And the synthesis the &#8220;push&#8221; phase (see <a href=\"http:\/\/wwwcg.in.tum.de\/fileadmin\/user_upload\/Lehrstuehle\/Lehrstuhl_XV\/Research\/Publications\/2009\/grapp09.pdf\">this paper<\/a> for an overview of &#8220;push pull&#8221; algorithm). As the scaling function is a \\(B_3\\) spline, the filter \\(h\\) is :<\/p>\n<ul>\n<li>\\(\\left(\\frac{1}{16},\\frac{1}{4},\\frac{3}{8},\\frac{1}{4},\\frac{1}{16}\\right)\\) in \\(\\mathbb{R}\\)<\/li>\n<li>\\(<br \/>\n\\begin{pmatrix}<br \/>\n\\frac{1}{256} &amp; \\frac{1}{64} &amp; \\frac{3}{128} &amp; \\frac{1}{64} &amp; \\frac{1}{256} \\\\<br \/>\n\\frac{1}{64} &amp; \\frac{1}{16} &amp; \\frac{3}{32} &amp; \\frac{1}{16} &amp; \\frac{1}{64} \\\\<br \/>\n\\frac{3}{128} &amp; \\frac{3}{32} &amp; \\frac{9}{64} &amp; \\frac{3}{32} &amp; \\frac{3}{128} \\\\<br \/>\n\\frac{1}{64} &amp; \\frac{1}{16} &amp; \\frac{3}{32} &amp; \\frac{1}{16} &amp; \\frac{1}{64} \\\\<br \/>\n\\frac{1}{256} &amp; \\frac{1}{64} &amp; \\frac{3}{128} &amp; \\frac{1}{64} &amp; \\frac{1}{256}<br \/>\n\\end{pmatrix}<br \/>\n\\) in \\(\\mathbb{R}^2\\)<\/li>\n<\/ul>\n<p>If we apply this to image filtering, we will have:<\/p>\n<ol>\n<li>\\(c_0 = I\\)<\/li>\n<li>\\(\\forall i \\in [0,n[\\) and \\(\\forall p \\in [0,S_x] \\times [0,S_y[\\)\n<ul>\n<li>\\(c_{i+1}(p) = \\sum_{0 \\leq y &lt; 5} \\sum_{0 \\leq x &lt; 5} h_i(x,y) c_i\\left(p+ 2^i(x,y)\\right)\\)<\/li>\n<li>\\(d_i(p) = c_i &#8211; c_{i+1}(p)\\)<\/li>\n<\/ul>\n<\/li>\n<li>\\(I&#8217;=c_n + \\sum_{0 \\leq i &lt; n}d_i\\)<\/li>\n<\/ol>\n<p>Edge awareness is acheived by adding a Gaussian distribution based weighing function of the colorimetric difference between the sample and source pixel (just like any standard bilateral filter).<\/p>\n<ul>\n<li>\\(\\omega_{\\sigma}(u,v) = e^{\\frac{\\left\\|c_i(u) &#8211; c_i(v)\\right\\|^2}{\\sigma}}\\)<\/li>\n<\/ul>\n<p><strong>2<\/strong>will be:<\/p>\n<ul>\n<li>\\(c_{i+1}(p) = \\frac{1}{W} \\sum \\sum \\omega_{\\sigma}\\left(c_i(p),c_i\\left(p+2^i(x,y)\\right)\\right)h_i(x,y) c_i(p+2^i(x,y))\\)<\/li>\n<li>with \\(W = \\sum\\sum \\omega_{\\sigma}\\left(c_i(p),c_i\\left(p+2^i(x,y)\\right)\\right)h_i(x,y)\\)<\/li>\n<\/ul>\n<p>Denoising can be acheived by thresholding the detail images.<\/p>\n<ul>\n<li>\\(d&#8217;_i = max\\left(0, |d_i|-\\tau\\right) \\cdot sign(d_i)\\)<\/li>\n<\/ul>\n<p>More details can be found in &#8220;Edge-Optimized \u00c0-TrousWavelets for Local Contrast<br \/>\nEnhancement with Robust Denoising&#8221; by Johannes Hanika, Holger Dammertz and Hendrik Lensch (<a href=\"http:\/\/www.darktable.org\/wp-content\/uploads\/2011\/11\/hdl11_paper.pdf\">paper<\/a>, <a href=\"http:\/\/www.darktable.org\/wp-content\/uploads\/2011\/11\/hdl11_talk.pdf\">slides<\/a>).<br \/>\nImplementing it in OpenGL\/glsl is pretty straightforward. The detail textures will be stored in a texture array. And as we only need the last filtered image for the synthesis, we will ping-pong between 2 textures.<br \/>\nHere is how the texture array is created:<\/p>\n<pre class=\"brush: cpp; title: ; notranslate\" title=\"\">glBindTexture(GL_TEXTURE_2D_ARRAY, texId);\r\n\tglTexImage3D(GL_TEXTURE_2D_ARRAY, 0, GL_RGB32F, imageWidth, imageHeight, LEVEL_MAX, 0, GL_RGB, GL_UNSIGNED_BYTE, 0);\r\n\tfor(i=0; i&lt;LEVEL_MAX; i++)\r\n        {\r\n\t\tglTexSubImage3D(GL_TEXTURE_2D_ARRAY, 0, 0, 0, i, imageWidth, imageHeight, 1, GL_RGB, GL_UNSIGNED_BYTE, imageData);\r\n\t}\r\n\t\/\/ Set texture wrap and filtering here.\r\nglBindTexture(GL_TEXTURE_2D_ARRAY, 0);\r\n<\/pre>\n<p>You attach a layer just like a standard 3D texture.<\/p>\n<pre class=\"brush: cpp; title: ; notranslate\" title=\"\">glBindFramebuffer(GL_FRAMEBUFFER, framebuffer);\r\nglFramebufferTextureLayer(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, texId, 0, layerIndex);<\/pre>\n<p>You access it in a shader using a <em>sampler2DArray<\/em>. Here is an example based on the synthesis shader:<\/p>\n<pre class=\"brush: cpp; title: ; notranslate\" title=\"\">#extension GL_EXT_texture_array : enable\r\nuniform sampler2D filtered;\r\nuniform sampler2DArray details;\r\n\r\nout vec4 outData;\r\n\r\nvoid main()\r\n{\r\n    int levelMax = int(textureSize(uDetails, 0).z);\r\n    vec4 data = texelFetch(filtered, ivec2(gl_FragCoord.xy), 0);;\r\n\r\n    for(int i=0; i&lt;levelMax; i++)\r\n    {\r\n\t\tdata += texelFetch(details, ivec3(gl_FragCoord.xy, i), 0);\r\n    }\r\n\r\n    outData = data;\r\n}\r\n<\/pre>\n<p>I uploaded the source code of a little demo <a href=\"http:\/\/blockos.org\/releases\/wavelets\/atrous.zip\" target=\"_blank\">here<\/a>. It will filter the infamous mandrill image using a edge preserving \u00e0-trous wavelet transform. I intentionally used &#8220;extreme&#8221; parameters to show the effects of the edge avoiding filter. It comes with a Microsoft Visual Studion 2010 project. Note that you will need <a href=\"http:\/\/www.glfw.org\/\" target=\"_blank\">GLFW<\/a> and <a href=\"http:\/\/glew.sourceforge.net\/\" target=\"_blank\">GLEW<\/a>. You will also have to edit the <strong>ContribDir<\/strong> item in <strong>platform\\windows\\wavelet.vcxproj.props<\/strong>.<br \/>\nOn the left you have the original image. And on the right the filtered image with \\(\\sigma=1.125\\) and \\(\\tau=0.05125\\).<br \/>\n<img decoding=\"async\" loading=\"lazy\" alt=\"\" src=\"http:\/\/blockos.org\/releases\/wavelets\/mandrill.jpg\" title=\"Original mandrill\" class=\"alignnone\" width=\"250\" height=\"240\" \/> <img decoding=\"async\" loading=\"lazy\" alt=\"\" src=\"http:\/\/blockos.org\/releases\/wavelets\/mandrill_filtered.jpg\" title=\"Filtered mandrill\" class=\"alignnone\" width=\"250\" height=\"240\" \/><\/p>\n<h2>Stuffs<\/h2>\n<ul>\n<li><a href=\"http:\/\/blockos.org\/releases\/wavelets\/atrous.zip\" target=\"_blank\">source code<\/a><\/li>\n<\/ul>\n<h2>References<\/h2>\n<ul>\n<li>Gilbert Strang, <a href=\"http:\/\/videolectures.net\/mit18085f07_strang_lec27\/\" target=\"_blank\">Lecture 27: Multiresolution, wavelet transform and scaling function<\/a><\/li>\n<li>Albert Bijaoui, Jean-Luc Starck, Fionn Murtagh, <a href=\"http:\/\/documents.irevues.inist.fr\/handle\/2042\/1863\" target=\"_blank\">Restauration des images multi-\u00e9chelles par l&#8217;algorithme \u00e0 trous<\/a> (french)<\/li>\n<li>Holger Dammertz, Daniel Sewtz, Johannes Hanika, Hendrik P.A. Lensch, <a href=\"http:\/\/www.uni-ulm.de\/fileadmin\/website_uni_ulm\/iui.inst.100\/institut\/Papers\/atrousGIfilter.pdf\" target=\"_blank\">Edge-Avoiding \u00c0-Trous Wavelet Transform for fast Global<br \/>\nIllumination Filtering<\/a>, <a href=\"http:\/\/www.highperformancegraphics.org\/previous\/www_2010\/media\/RayTracing_I\/HPG2010_RayTracing_I_Dammertz.pdf\" target=\"_blank\">(slides)<\/a><\/li>\n<li>Johannes Hanika, Holger Dammertz, Hendrik Lensch, <a href=\"http:\/\/www.darktable.org\/wp-content\/uploads\/2011\/11\/hdl11_paper.pdf\" target=\"_blank\">Edge-Optimized \u00c0-Trous Wavelets for Local Contrast Enhancement with Robust Denoising<\/a>, <a href=\"http:\/\/www.darktable.org\/wp-content\/uploads\/2011\/11\/hdl11_talk.pdf\" target=\"_blank\">(slides)<\/a>\n<li>Peter-Pike Sloan, Naga K. Govindaraju, Derek Nowrouzezahrai, John Snyder,  <a href=\"http:\/\/www.ppsloan.org\/publications\/ProxyPG.pdf\" target=\"_blank\">Image-Based Proxy Accumulation for Real-Time Soft Global Illumination<\/a><\/li>\n<li>Martin Kraus, <a href=\"http:\/\/wwwcg.in.tum.de\/fileadmin\/user_upload\/Lehrstuehle\/Lehrstuhl_XV\/Research\/Publications\/2009\/grapp09.pdf\" target=\"_blank\">The Pull-Push Algorithm Revisited<\/a><\/li>\n<\/ul>\n","protected":false},"excerpt":{"rendered":"<p>I first came across the \u00e0 trous algorithm when I was looking for a good edge avoiding filter for SSAO. Back then I was trying to code a demo. I had depth of field, some kind of hdr, spherical harmonics ambient lighting and screen space ambient transfer. I tried several\u2026 <a class=\"continue-reading-link\" href=\"https:\/\/blog.blockos.org\/?p=455\">Continue reading<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[3,5],"tags":[27,11,29],"_links":{"self":[{"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/posts\/455"}],"collection":[{"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=455"}],"version-history":[{"count":234,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/posts\/455\/revisions"}],"predecessor-version":[{"id":1045,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=\/wp\/v2\/posts\/455\/revisions\/1045"}],"wp:attachment":[{"href":"https:\/\/blog.blockos.org\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=455"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=455"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blog.blockos.org\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=455"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}