{"id":140,"date":"2007-06-08T14:09:55","date_gmt":"2007-06-08T18:09:55","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/"},"modified":"2019-10-23T12:57:44","modified_gmt":"2019-10-23T16:57:44","slug":"upslope-area-part-5","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/06\/08\/upslope-area-part-5\/","title":{"rendered":"Upslope area – Plateau processing"},"content":{"rendered":"
\r\n

[NOTE: After I originally wrote this post about handling plateaus, I discovered problems with the method, so I am no longer using it. I revisited the plateau issue in a later post<\/a>. -SE]<\/em>\r\n<\/p>\r\n\r\n

Today I'll continue thinking about plateaus in a DEM and what they mean for calculating pixel flow. Here's a very tiny portion\r\n of my DEM array, just northwest of North Pond<\/a>.\r\n <\/p>

a = [120 118 123; 120 118 125; 119 119 126]<\/pre>
\r\na =\r\n\r\n   120   118   123\r\n   120   118   125\r\n   119   119   126\r\n\r\n<\/pre>
The center pixel doesn't have a downhill neighbor, so I'm calling it a plateau pixel<\/i>.  No matter that it's a very small plateau; we can't compute a pixel flow direction for it, and that will cause problems\r\n      when we get around to computing upslope area.\r\n   <\/p>\r\n   
One technique used in image processing to deal with plateaus is called the lower complete transformation<\/i>.  A lower complete image<\/i> is one where the only pixels having no downhill neighors are the pixels belonging to regional minima.  Roughly speaking,\r\n      the lower complete transformation operates by elevating plateau pixels.  Pixels further from the plateau edge get elevated\r\n      more.  Of course, pixels uphill from the plateau need to be elevated as well so that they remain uphill.  I have read about\r\n      this technique but have never implemented it myself.  (For more information, I recommend P. Soille, Morphological Image Analysis: Principles and Applications<\/i>, 2nd edition, Springer, 2003, pp. 222-226.)\r\n   <\/p>\r\n   
I'm going to use a different technique, one that's based on the algorithm in the Image Processing Toolbox function roifill<\/tt>.  This function fills in pixels in a particular region by smoothly interpolating from the pixels surrounding the region.\r\n       It does so by solving a sparse linear system based on the notion that, in the filled output, each pixel in the filled region\r\n      equals the average of its north, east, south, and west neighbors.  (It's kind of an interesting method.  Maybe that'll be\r\n      a future blog topic.)\r\n   <\/p>\r\n   
Here's the example (slightly modified) from the roifill<\/tt> documentation:\r\n   <\/p>
I = imread('eight.tif'<\/span>);\r\nimshow(I)\r\n\r\nc = [222 272 300 270 221 194 222];\r\nr = [21 21 75 121 121 75 21];\r\nhold on<\/span>\r\nplot(c,r,'r'<\/span>,'LineWidth'<\/span>,4)<\/pre> 
J = roifill(I,c,r);\r\nimshow(J)<\/pre> 
And here's what roifill<\/tt> does to our tiny snippet from the DEM:\r\n   <\/p>
mask = true(3, 3);\r\nb = roifill(a, mask)<\/pre>
\r\nb =\r\n\r\n  120.0000  118.0000  123.0000\r\n  120.0000  120.5000  125.0000\r\n  119.0000  119.0000  126.0000\r\n\r\n<\/pre>
You can see that the center pixel has been replaced by a new value, and that the center pixel now has a downhill neighbor.<\/p>\r\n   
Here's a slightly bigger example, where the two center pixels have no downhill neighbors:<\/p>
a2 = pond(15:17, 19:22)<\/pre>
\r\na2 =\r\n\r\n   116   116   116   120\r\n   117   115   115   117\r\n   119   115   115   116\r\n\r\n<\/pre>
And the output of roifill<\/tt><\/p>
mask2 = true(3, 4);\r\nb2 = roifill(a2, mask2)<\/pre>
\r\nb2 =\r\n\r\n   116   116   116   120\r\n   117   116   116   117\r\n   119   115   115   116\r\n\r\n<\/pre>
Well, the two center pixels now have downhill neighbors, but WAIT!  The (1,2) pixel no longer has a downhill neighbor!  OK,\r\n      time to 'fess up.  I didn't really anticipate that this would be a problem.  I'll have to give this some thought.  My first\r\n      reaction is: It'll be OK because we can use the original, unmodified DEM values to compute the pixel flow for the nonplateau\r\n      pixels.  We only need the modified DEM values to compute pixel flow for the plateau pixels.  But I'll sleep on it.\r\n   <\/p>\r\n   
Assuming I can resolve that problem successfully, we can use roifill<\/tt> to compute pixel flow directions for the midplateau pixels.  It won't help us, though, for the plateau pixels that are part\r\n      of regional maxima or regional minima.  I'll figure out what to do about those next week.\r\n   <\/p>