{"id":151,"date":"2007-07-13T10:46:36","date_gmt":"2007-07-13T14:46:36","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/07\/13\/upslope-area-part-6\/"},"modified":"2019-10-23T13:41:15","modified_gmt":"2019-10-23T17:41:15","slug":"upslope-area-part-6","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/07\/13\/upslope-area-part-6\/","title":{"rendered":"Upslope area &#8211; More plateau processing"},"content":{"rendered":"<div xmlns:mwsh=\"https:\/\/www.mathworks.com\/namespace\/mcode\/v1\/syntaxhighlight.dtd\" class=\"content\">\r\n   <introduction>\r\n      <p>In <a href=\"https:\/\/blogs.mathworks.com\/steve\/2007\/05\/31\/upslope-area-part-4\/\">part 4<\/a> of this series I classified DEM plateaus into:\r\n      <\/p>\r\n      <div>\r\n         <ul>\r\n            <li>Max plateaus<\/li>\r\n            <li>Min plateaus<\/li>\r\n            <li>Mid plateaus<\/li>\r\n         <\/ul>\r\n      <\/div>\r\n      <p>In <a href=\"https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/\">part 5<\/a> I outlined the possibility of using <tt>roifill<\/tt> to modify plateau pixels.\r\n      <\/p>\r\n      <p>Today I'll combine these ideas with a little morphological processing to compute pixel flow directions for all DEM pixels,\r\n         except for isolated ones in the center of the min plateaus.\r\n      <\/p>\r\n   <\/introduction>\r\n   <p>First, find all the plateau pixels that have no downhill neighbors. Remember that we can find these using erosion:<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">s = load(<span style=\"color: #A020F0\">'upslope_sample_dem'<\/span>);\r\npond = s.Zm(140:280, 160:230);\r\nplateaus = imerode(pond, ones(3,3)) == pond;\r\n\r\nimshow(plateaus, <span style=\"color: #A020F0\">'initialmag'<\/span>, <span style=\"color: #A020F0\">'fit'<\/span>)\r\ntitle(<span style=\"color: #A020F0\">'Plateau pixels'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/151\/upslope_area_6_01.png\"> <p>Next, find the max plateaus using <tt>imregionalmax<\/tt>:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">max_plateaus = imregionalmax(pond);<\/pre><p>One way to visualize plateau pixels is to overlay them on the original image using another color.  Use <tt>imoverlay<\/tt>, a utility function I wrote last year.  (See my <a href=\"https:\/\/blogs.mathworks.com\/steve\/2006\/03\/28\/image-overlays\/\">March 28, 2006 post<\/a> for details.)\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">vis1 = imoverlay(mat2gray(pond), max_plateaus, [1 0 0]);\r\nimshow(vis1, <span style=\"color: #A020F0\">'initialmag'<\/span>, <span style=\"color: #A020F0\">'fit'<\/span>)\r\ntitle(<span style=\"color: #A020F0\">'Max plateaus in red'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/151\/upslope_area_6_02.png\"> <p>The <tt>bwmorph<\/tt> <i>shrink<\/i> operation reduces connected components to single points (or, for objects with holes, to a ring).\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">max_plateaus_shrunk = bwmorph(max_plateaus, <span style=\"color: #A020F0\">'shrink'<\/span>, Inf);\r\n\r\nvis2 = imoverlay(vis1, max_plateaus_shrunk, [1 1 0]);\r\nimshow(vis2, <span style=\"color: #A020F0\">'initialmag'<\/span>, <span style=\"color: #A020F0\">'fit'<\/span>)\r\ntitle(<span style=\"color: #A020F0\">'Shrunken plateau pixels in yellow'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/151\/upslope_area_6_03.png\"> <p>Bump the shrunken plateau pixels to an elevation slightly above their neighbors.  We'll bump up by a fixed constant here.\r\n       To make a completely general version we'll probably have to calculate the bump based on eps and the maximum plateau height.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">pond2 = pond;\r\npond2(max_plateaus_shrunk) = pond(max_plateaus_shrunk) + 0.1;<\/pre><p>Now let's do the same thing to the min plateaus, except that we'll slightly lower the shrunken plateau pixels.<\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">min_plateaus = imregionalmin(pond);\r\nmin_plateaus_shrunk = bwmorph(min_plateaus, <span style=\"color: #A020F0\">'shrink'<\/span>, Inf);\r\npond2(min_plateaus_shrunk) = pond2(min_plateaus_shrunk) - 0.1;<\/pre><p>As I described in <a href=\"https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/\">part 5<\/a>, I would like to use <tt>roifill<\/tt> here to interpolate the plateau pixels.   However, the logic of <tt>roifill<\/tt> doesn't quite work out right. Specifically, when you pass in a mask to <tt>roifill<\/tt>, it only fills in the <i>interior<\/i> pixels of the mask.  You can see this in the code:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">roifill<\/span> <span style=\"color: #A020F0\">68:74<\/span><\/pre><pre style=\"font-style:oblique\">\r\n68    % Find the perimeter pixels of the specified region; these will\r\n69    % be used to form the boundary conditions for the soap film PDE.\r\n70    perimeter = bwperim(mask);\r\n71    \r\n72    % Find the interior pixels; these are the pixels that will be\r\n73    % replaced in the output image.\r\n74    interior = mask &amp; ~perimeter;\r\n\r\n<\/pre><p>For this application I want to fill all the pixels specified by the mask.  So I'll do what MATLAB users often do: I'll take\r\n      an existing M-file and modify it to suit my needs.  Starting with <tt>roifill<\/tt>, I created a new M-file called <tt>regionfill<\/tt>.  In the new function, <tt>mask<\/tt> specifies the pixels to be filled, and the boundary pixels are computed from the mask by dilation instead of by using <tt>bwperim<\/tt>:\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">dbtype <span style=\"color: #A020F0\">regionfill<\/span> <span style=\"color: #A020F0\">27:29<\/span><\/pre><pre style=\"font-style:oblique\">\r\n27    % Find the N, E, S, and W pixels adjacent to the pixels to be filled. These\r\n28    % will be used to form the boundary conditions for the soap film PDE.\r\n29    perimeter = imdilate(mask, [0 1 0; 1 1 1; 0 1 0]) &amp; ~mask;\r\n\r\n<\/pre><p>Now I can use <tt>regionfill<\/tt> to interpolate the plateau pixels. But not all the plateau pixels. We don't want to fill those that got bumped higher or\r\n      lower.  When those bumped pixels get used by <tt>regionfill<\/tt> as boundary conditions, they'll help eliminate the other plateau pixels.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">mask = plateaus &amp; ~max_plateaus_shrunk &amp; ~min_plateaus_shrunk;\r\n\r\npond2 = regionfill(pond2, mask);<\/pre><p>Using both the original DEM (<tt>pond<\/tt>) and the modified DEM (<tt>pond2<\/tt>), we can now compute pixel flow for every pixel.  We use the modified DEM only for plateau pixels (see the last two paragraphs\r\n      of <a href=\"https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/\">part 5<\/a> for an explanation).\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">[M,N] = size(pond);\r\nr = zeros(M-2, N-2);\r\ns = zeros(M-2, N-2);\r\n<span style=\"color: #0000FF\">for<\/span> i = 1:M-2\r\n    <span style=\"color: #0000FF\">for<\/span> j = 1:N-2\r\n        [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1);\r\n        <span style=\"color: #0000FF\">if<\/span> isnan(r(i,j))\r\n            <span style=\"color: #228B22\">% Plateau pixel; use the modified DEM.<\/span>\r\n            [r(i,j), s(i,j)] = pixelFlow(pond2,i+1,j+1);\r\n        <span style=\"color: #0000FF\">end<\/span>\r\n    <span style=\"color: #0000FF\">end<\/span>\r\n<span style=\"color: #0000FF\">end<\/span><\/pre><p>As the images below show, this procedure has successfully computed pixel flow directions for all pixels except those belonging\r\n      to the shrunken min plateaus.\r\n   <\/p><pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid rgb(200,200,200)\">subplot(1,2,1)\r\nimshow(isnan(r))\r\ntitle(<span style=\"color: #A020F0\">'No pixel flow direction'<\/span>)\r\n\r\nsubplot(1,2,2)\r\nimshow(min_plateaus_shrunk(2:end-1,2:end-1))\r\ntitle(<span style=\"color: #A020F0\">'Min plateaus'<\/span>)<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/steve\/151\/upslope_area_6_04.png\"> <p>Next time I'll start looking at the actual upslope area computation.<\/p><script language=\"JavaScript\">\r\n<!--\r\n\r\n    function grabCode_deac112114944535a29ddee12572c8f6() {\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='deac112114944535a29ddee12572c8f6 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' deac112114944535a29ddee12572c8f6';\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        author = 'Steve Eddins';\r\n        copyright = 'Copyright 2007 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 author and copyright lines at the bottom if specified.\r\n        if ((author.length > 0) || (copyright.length > 0)) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (author.length > 0) {\r\n                d.writeln('% _' + author + '_');\r\n            }\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      \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_deac112114944535a29ddee12572c8f6()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n            the MATLAB code \r\n            <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.4<br><\/p>\r\n<\/div>\r\n<!--\r\ndeac112114944535a29ddee12572c8f6 ##### SOURCE BEGIN #####\r\n%%\r\n% In <https:\/\/blogs.mathworks.com\/steve\/2007\/05\/31\/upslope-area-part-4\/ \r\n% part 4> of this series I classified DEM plateaus into:\r\n%\r\n% * Max plateaus\r\n% * Min plateaus\r\n% * Mid plateaus\r\n%\r\n% In <https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/ \r\n% part 5> I outlined the possibility of using |roifill| to modify plateau\r\n% pixels.\r\n%\r\n% Today I'll combine these ideas with a little morphological processing to\r\n% compute pixel flow directions for all DEM pixels, except for isolated\r\n% ones in the center of the min plateaus.\r\n\r\n%%\r\n% First, find all the plateau pixels that have no downhill neighbors.\r\n% Remember that we can find these using erosion:\r\n\r\ns = load('upslope_sample_dem');\r\npond = s.Zm(140:280, 160:230);\r\nplateaus = imerode(pond, ones(3,3)) == pond;\r\n\r\nimshow(plateaus, 'initialmag', 'fit')\r\ntitle('Plateau pixels')\r\n\r\n%%\r\n% Next, find the max plateaus using |imregionalmax|:\r\n\r\nmax_plateaus = imregionalmax(pond);\r\n\r\n%%\r\n% One way to visualize plateau pixels is to overlay them on the original\r\n% image using another color.  Use |imoverlay|, a utility function I wrote\r\n% last year.  (See my\r\n% <https:\/\/blogs.mathworks.com\/steve\/2006\/03\/28\/image-overlays\/ March 28, \r\n% 2006 post>\r\n% for details.)\r\n\r\nvis1 = imoverlay(mat2gray(pond), max_plateaus, [1 0 0]);\r\nimshow(vis1, 'initialmag', 'fit')\r\ntitle('Max plateaus in red')\r\n\r\n%%\r\n% The |bwmorph| _shrink_ operation reduces connected components to single\r\n% points (or, for objects with holes, to a ring).\r\n\r\nmax_plateaus_shrunk = bwmorph(max_plateaus, 'shrink', Inf);\r\n\r\nvis2 = imoverlay(vis1, max_plateaus_shrunk, [1 1 0]);\r\nimshow(vis2, 'initialmag', 'fit')\r\ntitle('Shrunken plateau pixels in yellow')\r\n\r\n%%\r\n% Bump the shrunken plateau pixels to an elevation slightly above their\r\n% neighbors.  We'll bump up by a fixed constant here.  To make a completely\r\n% general version we'll probably have to calculate the bump based on eps\r\n% and the maximum plateau height.\r\n\r\npond2 = pond;\r\npond2(max_plateaus_shrunk) = pond(max_plateaus_shrunk) + 0.1;\r\n\r\n%%\r\n% Now let's do the same thing to the min plateaus, except that we'll\r\n% slightly lower the shrunken plateau pixels.\r\n\r\nmin_plateaus = imregionalmin(pond);\r\nmin_plateaus_shrunk = bwmorph(min_plateaus, 'shrink', Inf);\r\npond2(min_plateaus_shrunk) = pond2(min_plateaus_shrunk) - 0.1;\r\n\r\n%%\r\n% As I described in \r\n% <https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/ part 5>, \r\n% I would like to use |roifill| here to\r\n% interpolate the plateau pixels.   However, the logic of |roifill| doesn't\r\n% quite work out right. Specifically, when you pass in a mask to |roifill|,\r\n% it only fills in the _interior_ pixels of the mask.  You can see this in\r\n% the code:\r\n\r\ndbtype roifill 68:74\r\n\r\n%%\r\n% For this application I want to fill all the pixels specified\r\n% by the mask.  So I'll do what MATLAB users often do: I'll take an\r\n% existing M-file and modify it to suit my needs.  Starting with |roifill|,\r\n% I created a new M-file called |regionfill|.  In the new function, |mask|\r\n% specifies the pixels to be filled, and the boundary pixels are computed\r\n% from the mask by dilation instead of by using |bwperim|:\r\n\r\ndbtype regionfill 27:29\r\n\r\n%%\r\n% Now I can use |regionfill| to interpolate the plateau pixels. But not all\r\n% the plateau pixels. We don't want to fill those that got bumped higher or\r\n% lower.  When those bumped pixels get used by |regionfill| as boundary\r\n% conditions, they'll help eliminate the other plateau pixels.\r\n\r\nmask = plateaus & ~max_plateaus_shrunk & ~min_plateaus_shrunk;\r\n\r\npond2 = regionfill(pond2, mask);\r\n\r\n%%\r\n% Using both the original DEM (|pond|) and the modified DEM (|pond2|), we\r\n% can now compute pixel flow for every pixel.  We use the modified DEM only\r\n% for plateau pixels (see the last two paragraphs of \r\n% <https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/ part 5>\r\n% for an explanation).\r\n\r\n[M,N] = size(pond);\r\nr = zeros(M-2, N-2);\r\ns = zeros(M-2, N-2);\r\nfor i = 1:M-2\r\n    for j = 1:N-2\r\n        [r(i,j), s(i,j)] = pixelFlow(pond,i+1,j+1);\r\n        if isnan(r(i,j))\r\n            % Plateau pixel; use the modified DEM.\r\n            [r(i,j), s(i,j)] = pixelFlow(pond2,i+1,j+1);\r\n        end\r\n    end\r\nend\r\n\r\n%%\r\n% As the images below show, this procedure has successfully computed pixel\r\n% flow directions for all pixels except those belonging to the shrunken min\r\n% plateaus.\r\n\r\nsubplot(1,2,1)\r\nimshow(isnan(r))\r\ntitle('No pixel flow direction')\r\n\r\nsubplot(1,2,2)\r\nimshow(min_plateaus_shrunk(2:end-1,2:end-1))\r\ntitle('Min plateaus')\r\n\r\n%%\r\n% Next time I'll start looking at the actual upslope area computation.\r\n\r\n##### SOURCE END ##### deac112114944535a29ddee12572c8f6\r\n-->","protected":false},"excerpt":{"rendered":"<p>\r\n   \r\n      In part 4 of this series I classified DEM plateaus into:\r\n      \r\n      \r\n         \r\n        ... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/steve\/2007\/07\/13\/upslope-area-part-6\/\">read more >><\/a><\/p>","protected":false},"author":42,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[9],"tags":[86,124,246,364,366,36,376,362,372,288,374,72,52],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/151"}],"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=151"}],"version-history":[{"count":1,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/151\/revisions"}],"predecessor-version":[{"id":3546,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/posts\/151\/revisions\/3546"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/media?parent=151"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/categories?post=151"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/steve\/wp-json\/wp\/v2\/tags?post=151"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}