{"id":193,"date":"2009-07-22T13:48:24","date_gmt":"2009-07-22T13:48:24","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/22\/computational-geometry-in-matlab-r2009a-part-ii\/"},"modified":"2021-12-09T21:16:25","modified_gmt":"2021-12-10T02:16:25","slug":"computational-geometry-in-matlab-r2009a-part-ii","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/22\/computational-geometry-in-matlab-r2009a-part-ii\/","title":{"rendered":"Computational Geometry in MATLAB R2009a, Part II"},"content":{"rendered":"<div class=\"content\">\n<p>I am pleased to welcome back Damian Sheehy for the sequel to his earlier Computational Geometry <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p192\">post<\/a>.<\/p>\n<p>&nbsp;<\/p>\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li><a href=\"#1\">What's your problem! ?<\/a><\/li>\n<li><a href=\"#2\">griddata Produces \"Funny\" Results<\/a><\/li>\n<li><a href=\"#9\">Scattered Data Interpolation in the Real-World<\/a><\/li>\n<li><a href=\"#20\">What's the Bottom Line?<\/a><\/li>\n<li><a href=\"#24\">Feedback<\/a><\/li>\n<\/ul>\n<\/div>\n<h3>What's your problem! ?<a name=\"1\"><\/a><\/h3>\n<p>Whenever I have a question about MATLAB that I can't figure out from the help or the doc, I only need to stick my head out<br \/>\nof my office and ask my neighbor. Sometimes I wonder how I would get answers to my most challenging how-do-I questions if<br \/>\nI were an isolated MATLAB user tucked away in a lab. I guess the chances are I would get them from the same source, but the<br \/>\nmedium would be via the <a>MATLAB newsgroup<\/a> (comp.soft-sys.matlab).<\/p>\n<p>When I was new to MATLAB I would browse the newsgroup to learn. I found that learning from other people's questions was even<br \/>\nmore powerful than learning from their mistakes. As my knowledge of MATLAB has improved, I browse to get an insight into the<br \/>\nusecases and usability problems in my subject area and try to help out if I can. There aren't many posts on computational<br \/>\ngeometry, but I have clicked on quite a few <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/griddata.html\"><tt>griddata<\/tt><\/a> subject lines and I was able to learn that for some users, <tt>griddata<\/tt> didn't quite \"click\". Now this wasn't because they were slow; it was because <tt>griddata<\/tt> was slow or deficient in some respect.<\/p>\n<p>Despite its shortcomings, <tt>griddata<\/tt> is fairly popular, but it does trip up users from time to time. A new user can make a quick pass through the help example<br \/>\nand they are up and running. They use it to interpolate their own data, get good results, and they are hooked. In fact, it<br \/>\nmay work so well in black-box mode that users don't concern themselves about the internals; that is, until it starts returning<br \/>\n<a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/nan.html\"><tt>NaNs<\/tt><\/a>, or the interpolated surface looks bad close to the boundary, or maybe it looks bad overall. Then there's the user who fully<br \/>\nunderstands what's going on inside the black box, but lacks the functionality to do the job efficiently. This user shares<br \/>\nmy frustrations; the black box is too confined and needs to be reusable and flexible so it can be applied to general interpolation<br \/>\nscenarios.<\/p>\n<h3>griddata Produces \"Funny\" Results<a name=\"2\"><\/a><\/h3>\n<p><tt>griddata<\/tt> is useful, but it's not magic. You cannot throw in any old data and expect to get back a nice looking smooth surface like<br \/>\n<a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/peaks.html\"><tt>peaks<\/tt><\/a>. Like everything else in MATLAB it doesn't work like that. Sometimes our lack of understanding of what goes on in the inside<br \/>\ndistorts our perception of reality and our hopes to get lucky become artificially elevated. It's tempting to go down the path<br \/>\nof diagnosing unexpected results from <tt>griddata<\/tt>, but I will have to leave that battle for another day. Besides, if you are reading this blog you are probably a savvy MATLAB<br \/>\nuser who is already aware of the usual pitfalls. But, I would like to provide one tip; if 2D scattered data interpolation<br \/>\nproduces unexpected results, create a Delaunay triangulation of the <tt>(x, y)<\/tt> locations and plot the triangulation. Ideally we would like to see triangles that are close to equilateral. If you see sliver<br \/>\nshaped triangles then you cannot expect to get good results in that neighborhood. This example shows what I mean. I will use<br \/>\ntwo new features from R2009a; <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html\"><tt>DelaunayTri<\/tt><\/a> and <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html\"><tt>TriScatteredInterp<\/tt><\/a>, but you can use the functions <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/delaunay.html\"><tt>delaunay<\/tt><\/a> and <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/griddata.html\"><tt>griddata<\/tt><\/a> if you are running an earlier version.<\/p>\n<p>I will create a distribution of points (red stars) that would produce concave boundaries if we were to join them like a dot-to-dot<br \/>\npuzzle. I will then create a Delaunay triangulation from these points and plot it (blue triangles).<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">t= 0.4*pi:0.02:0.6*pi;\r\nx1 = cos(t)';\r\ny1 = sin(t)'-1.02;\r\nx2 = x1;\r\ny2 = y1*(-1);\r\nx3 = linspace(-0.3,0.3,16)';\r\ny3 = zeros(16,1);\r\nx = [x1;x2;x3];\r\ny = [y1;y2;y3];\r\ndt = DelaunayTri(x,y);\r\nplot(x,y, <span style=\"color: #a020f0;\">'*r'<\/span>)\r\naxis <span style=\"color: #a020f0;\">equal<\/span>\r\nhold <span style=\"color: #a020f0;\">on<\/span>\r\ntriplot(dt)\r\nplot(x1,y1,<span style=\"color: #a020f0;\">'-r'<\/span>)\r\nplot(x2,y2,<span style=\"color: #a020f0;\">'-r'<\/span>)\r\nhold <span style=\"color: #a020f0;\">off<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/193\/CgBlogPart2_01.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>The triangles within the red boundaries are relatively well shaped; they are constructed from points that are in close proximity.<br \/>\nWe would expect the interpolation to work well in this region. Outside the red boundary, the triangles are sliver-like and<br \/>\nconnect points that are remote from each other. We do not have sufficient sampling in here to accurately capture the surface;<br \/>\nwe would expect the results to be poor in these regions.<\/p>\n<p>Let's see how the interpolation would look if we lifted the points to the surface of a paraboloid.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">z = x.^2 + y.^2;\r\nF = TriScatteredInterp(x,y,z);\r\n[xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688);\r\nzi = F(xi,yi);\r\nmesh(xi,yi,zi)\r\nxlabel(<span style=\"color: #a020f0;\">'Interpolated surface'<\/span>, <span style=\"color: #a020f0;\">'fontweight'<\/span>,<span style=\"color: #a020f0;\">'b'<\/span>);<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/193\/CgBlogPart2_02.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Now let's see what the actual surface looks like in this domain.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">zi = xi.^2 + yi.^2;\r\nmesh(xi,yi,zi)\r\nxlabel(<span style=\"color: #a020f0;\">'Actual surface'<\/span>, <span style=\"color: #a020f0;\">'fontweight'<\/span>,<span style=\"color: #a020f0;\">'b'<\/span>);<\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/193\/CgBlogPart2_03.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>OK, you get the idea. You can see the parts of the interpolated surface that correspond to the sliver triangles. In 3D, visual<br \/>\ninspection of the triangulation gets a bit trickier, but looking at the point distribution can often help illustrate potential<br \/>\nproblems.<\/p>\n<h3>Scattered Data Interpolation in the Real-World<a name=\"9\"><\/a><\/h3>\n<p>Have you ever taken a vacation overseas and experienced the per-transaction exchange-rate penalty for credit card purchases?<br \/>\nOh yeah, you're on your way out the door of the gift shop with your bag of goodies, you realize you forgot someone and . .<br \/>\n. . catchinggg catchinggg; you're hit with the exchange rate penalty again. But hey, it's a vacation, forget it. Besides,<br \/>\nI'm not Mr Savvy Traveler. Now when you're back at work in front of MATLAB, per-transaction penalties have a completely different<br \/>\nmeaning. We have to be savvy to avoid them, but in reality we would simply prefer not to pay them.<\/p>\n<p>The scattered data interpolation tools in MATLAB are very useful for analyzing scientific data. When we collect sample-data<br \/>\nfrom measurements, we often leverage the same equipment to gather various metrics of interest and sample over a period of<br \/>\ntime. Collecting weather data is a good example. We typically collect temperature, rainfall, and snowfall etc, at fixed locations<br \/>\nacross the country and we do this continuously. When we import this data into MATLAB we get to do the interesting part \u2013 analyze<br \/>\nit. In this scenario we could use interpolation to compute weather data at any (longitude, latitude) location in the country.<br \/>\nAssuming we have access to weather archives we could also make historical comparisons. For example, we could compute the total<br \/>\nrainfall at (longitude, latitude) for years 2008, 2007, 2006 etc, and likewise we could do the same for total snowfall. If<br \/>\nwe were to use <tt>griddata<\/tt> to perform the interpolation, then catchinggg catchinggg, we would pay the unnecessary overhead of computing the underlying<br \/>\nDelaunay triangulation each time a query is evaluated.<\/p>\n<p>For example, suppose we have the following annual rainfall data for 2006-2008:<\/p>\n<pre>     (xlon, ylat, annualRainfall06)\r\n     (xlon, ylat, annualRainfall07)\r\n     (xlon, ylat, annualRainfall08)<\/pre>\n<p><tt>xlon<\/tt>, <tt>ylat<\/tt>, and <tt>annualRainfall0*<\/tt> are vectors that represent the sample locations and the corresponding annual rainfall at these locations. We can use <tt>griddata<\/tt> to query by interpolation the annual rainfall <tt>zq<\/tt> at location <tt>(xq, yq)<\/tt>.<\/p>\n<pre>     zq06 = griddata(xlon, ylat, annualRainfall06, xq, yq)\r\n     zq07 = griddata(xlon, ylat, annualRainfall07, xq, yq)\r\n     zq08 = griddata(xlon, ylat, annualRainfall08, xq, yq)<\/pre>\n<p>Note that a Delaunay triangulation of the points <tt>(xlon, ylat)<\/tt> is computed each time the interpolation is evaluated at <tt>(xq, yq)<\/tt>. The same triangulation could have been reused for the 07 and 08 queries, since we only changed the values at each location.<br \/>\nIn this little illustrative example the penalty may not be significant because the number of samples is likely to be in the<br \/>\nhundreds or less. However, when we interpolate large data sets this represents a serious transaction penalty. <tt>TriScatteredInterp<\/tt>, the new scattered data interpolation tool released in R2009a, allows us to avoid this overhead. We can use it to construct<br \/>\nthe underlying triangulation once and make repeated low-cost queries. In addition, for fixed sample locations the sample values<br \/>\nand the interpolation method can be changed independently without re-triangulation. Repeating the previous illustration using<br \/>\n<tt>TriScatteredInterp<\/tt> we would get the following:<\/p>\n<pre>     F = TriScatteredInterp(xlon, ylat, annualRainfall06)\r\n     zq06 = F(xq, yq)\r\n     F.V = annualRainfall07;\r\n     zq07 = F(xq, yq)\r\n     F.V = annualRainfall08;\r\n     zq08 = F(xq, yq)<\/pre>\n<p>We can also change the interpolation method on-the-fly; this switches to the new Natural Neighbor interpolation method.<\/p>\n<pre>     F.Method = 'natural';<\/pre>\n<p>Again, no transaction penalty for switching the method. To evaluate using natural neighbor interpolation we do the following:<\/p>\n<pre>     zq08 = F(xq, yq)<\/pre>\n<p>As you can see, the interpolant works like a function handle making the calling syntax more intuitive.<\/p>\n<h3>What's the Bottom Line?<a name=\"20\"><\/a><\/h3>\n<p>The new scattered data interpolation function, <tt>TriScatteredInterp<\/tt>, allows you to interpolate large data sets very efficiently. Because the triangulation is cached within the black-box (the<br \/>\ninterpolant) you do not need to pay the overhead of recomputing the triangulation each time you evaluate, change the interpolation<br \/>\nmethod, or change the values at the sample locations.<\/p>\n<p>What's more, as in the case of <tt>DelaunayTri<\/tt>, this interpolation tool uses the <a href=\"http:\/\/www.cgal.org\">CGAL<\/a> Delaunay triangulations which are robust against numerical problems. The triangulation algorithms are also faster and more<br \/>\nmemory efficient. So if you use MATLAB to interpolate large scattered data sets, this is great news for you.<\/p>\n<p>The other bonus feature is the availability of the Natural Neighbor interpolation method \u2013 not be confused with Nearest Neighbor<br \/>\ninterpolation. Natural neighbor interpolation is a triangulation-based method that has an area-of-influence weighting associated<br \/>\nwith each sample point. It is widely used in the area of geostatistics because it is C1 continuous (except at the sample locations)<br \/>\nand it performs well in both clustered and sparse data locations.<\/p>\n<p>For more information on this feature and the new computational geometry features shipped in R2009a check out the overview <tt>video<\/tt>, and Delaunay triangulation demo, or follow these links to the documentation<\/p>\n<div>\n<ul>\n<li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/trirep.html\"><tt>TriRep<\/tt><\/a>,<\/li>\n<li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html\"><tt>DelaunayTri<\/tt><\/a>,<\/li>\n<li><a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html\"><tt>TriScatteredInterp<\/tt><\/a>.<\/li>\n<\/ul>\n<\/div>\n<h3>Feedback<a name=\"24\"><\/a><\/h3>\n<p>If you use MATLAB to perform scattered data interpolation, let me know what you think <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=193#respond\">here<\/a>. It's the feedback from users like you that prompts us to provide enhancements that are relevant to you.<\/p>\n<p><script>\/\/ <![CDATA[\nfunction grabCode_190f374fa98c4dc39fdcf0ac326de59d() {\n        \/\/ Remember the title so we can use it in the new page\n        title = document.title;\n\n        \/\/ Break up these strings so that their presence\n        \/\/ in the Javascript doesn't mess up the search for\n        \/\/ the MATLAB code.\n        t1='190f374fa98c4dc39fdcf0ac326de59d ' + '##### ' + 'SOURCE BEGIN' + ' #####';\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 190f374fa98c4dc39fdcf0ac326de59d';\n    \n        b=document.getElementsByTagName('body')[0];\n        i1=b.innerHTML.indexOf(t1)+t1.length;\n        i2=b.innerHTML.indexOf(t2);\n \n        code_string = b.innerHTML.substring(i1, i2);\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\n\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \n        \/\/ in the XML parser.\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\n        \/\/ doesn't go ahead and substitute the less-than character. \n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\n\n        author = 'Loren Shure';\n        copyright = 'Copyright 2009 The MathWorks, Inc.';\n\n        w = window.open();\n        d = w.document;\n        d.write('\n\n\n\n\n\n<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\n\n\n\n\n\n\\n');\n      \n      d.title = title + ' (MATLAB code)';\n      d.close();\n      }\n\/\/ ]]><\/script><\/p>\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\"><a><span style=\"font-size: x-small; font-style: italic;\">Get<br \/>\nthe MATLAB code<br \/>\n<noscript>(requires JavaScript)<\/noscript><\/span><\/a><\/p>\n<p>Published with MATLAB\u00ae 7.8<\/p>\n<\/div>\n<p><!--\n190f374fa98c4dc39fdcf0ac326de59d ##### SOURCE BEGIN #####\n%% Computational Geometry in MATLAB R2009a, Part II\n% I am pleased to welcome back Damian Sheehy for the sequel to his\n% earlier Computational Geometry <https:\/\/blogs.mathworks.com\/loren\/?p192 post>.\n\n%% What's your problem! ?\n% Whenever I have a question about MATLAB that I can't figure out from the\n% help or the doc, I only need to stick my head out of my office and ask my\n% neighbor. Sometimes I wonder how I would get answers to my most challenging\n% how-do-I questions if I were an isolated MATLAB user tucked away in a lab.\n% I guess the chances are I would get them from the same source, but the\n% medium would be via the\n% <http:\/\/ MATLAB newsgroup> (comp.soft-sys.matlab).\n%\n% When I was new to MATLAB I would browse the newsgroup to learn. I found that\n% learning from other people's questions was even more powerful than learning\n% from their mistakes. As my knowledge of MATLAB has improved, I browse to get an\n% insight into the usecases and usability problems in my subject area and try\n% to help out if I can. There aren't many posts on computational geometry,\n% but I have clicked on quite a few\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/griddata.html |griddata|>\n% subject lines and I was able to\n% learn that for some users, |griddata| didn't quite \"click\". Now this wasn't\n% because they were slow; it was because |griddata| was slow or deficient in\n% some respect.\n%\n% Despite its shortcomings, |griddata| is fairly popular, but it does trip up\n% users from time to time. A new user can make a quick pass through the help\n% example and they are up and running. They use it to interpolate their own\n% data, get good results, and they are hooked.  In fact, it may work so well\n% in black-box mode that users don't concern themselves about the internals;\n% that is, until it starts returning\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/nan.html |NaNs|>,\n% or the interpolated surface looks\n% bad close to the boundary, or maybe it looks bad overall. Then there's the\n% user who fully understands what's going on inside the black box, but lacks\n% the functionality to do the job efficiently. This user shares my frustrations;\n% the black box is too confined and needs to be reusable and flexible so it\n% can be applied to general interpolation scenarios.\n\n%% griddata Produces \"Funny\" Results\n% |griddata| is useful, but it's not magic. You cannot throw in any old data\n% and expect to get back a nice looking smooth surface like\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/peaks.html |peaks|>.\n% Like everything else in MATLAB it doesn't work like that.  Sometimes our\n% lack of understanding of what goes on in the inside distorts our perception\n% of reality and our hopes to get lucky become artificially elevated. It's\n% tempting to go down the path of diagnosing unexpected results from |griddata|,\n% but I will have to leave that battle for another day. Besides, if you are\n% reading this blog you are probably a savvy MATLAB user who is already aware\n% of the usual pitfalls. But, I would like to provide one tip; if 2D scattered\n% data interpolation produces unexpected results, create a Delaunay\n% triangulation of the |(x, y)| locations and plot the triangulation. Ideally\n% we would like to see triangles that are close to equilateral. If you see\n% sliver shaped triangles then you cannot expect to get good results in that\n% neighborhood. This example shows what I mean. I will use two new features\n% from R2009a;\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html |DelaunayTri|>\n% and\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html |TriScatteredInterp|>,\n% but you can use the functions\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/delaunay.html |delaunay|>\n% and\n% <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/ref\/griddata.html |griddata|>\n% if you are running an earlier version.\n%%\n% I will create a distribution of\n% points (red stars) that would produce concave boundaries if we were to join them like\n% a dot-to-dot puzzle. I will then create a Delaunay triangulation from\n% these points and plot it (blue triangles).\nt= 0.4*pi:0.02:0.6*pi;\nx1 = cos(t)';\ny1 = sin(t)'-1.02;\nx2 = x1;\ny2 = y1*(-1);\nx3 = linspace(-0.3,0.3,16)';\ny3 = zeros(16,1);\nx = [x1;x2;x3];\ny = [y1;y2;y3];\ndt = DelaunayTri(x,y);\nplot(x,y, '*r')\naxis equal\nhold on\ntriplot(dt)\nplot(x1,y1,'-r')\nplot(x2,y2,'-r')\nhold off\n%%\n% The triangles within the red boundaries are relatively well shaped; they\n% are constructed from points that are in close proximity. We would expect\n% the interpolation to work well in this region. Outside the red boundary,\n% the triangles are sliver-like and connect points that are remote from each\n% other. We do not have sufficient sampling in here to accurately capture\n% the surface; we would expect the results to be poor in these regions.\n%%\n% Let's see how the interpolation would look if we lifted the points to the\n% surface of a paraboloid.\n%%\nz = x.^2 + y.^2;\nF = TriScatteredInterp(x,y,z);\n[xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688);\nzi = F(xi,yi);\nmesh(xi,yi,zi)\nxlabel('Interpolated surface', 'fontweight','b');\n%%\n% Now let's see what the actual surface looks like in this domain.\nzi = xi.^2 + yi.^2;\nmesh(xi,yi,zi)\nxlabel('Actual surface', 'fontweight','b');\n%%\n% OK, you get the idea. You can see the parts of the interpolated surface\n% that correspond to the sliver triangles. In 3D, visual inspection of the\n% triangulation gets a bit trickier, but looking at the point distribution\n% can often help illustrate potential problems.\n\n%% Scattered Data Interpolation in the Real-World\n% Have you ever taken a vacation overseas and experienced the per-transaction\n% exchange-rate penalty for credit card purchases? Oh yeah, you're on your\n% way out the door of the gift shop with your bag of goodies, you realize\n% you forgot someone and . . . . catchinggg catchinggg; you're hit with the\n% exchange rate penalty again. But hey, it's a vacation, forget it. Besides,\n% I'm not Mr Savvy Traveler.\n% Now when you're back at work in front of MATLAB, per-transaction penalties\n% have a completely different meaning. We have to be savvy to avoid them,\n% but in reality we would simply prefer not to pay them.\n%%\n% The scattered data interpolation tools in MATLAB are very useful for\n% analyzing scientific data. When we collect sample-data from measurements,\n% we often leverage the same equipment to gather various metrics of interest\n% and sample over a period of time. Collecting weather data is a good example.\n% We typically collect temperature, rainfall, and snowfall etc, at fixed\n% locations across the country and we do this continuously. When we import\n% this data into MATLAB we get to do the interesting part \u00e2\u20ac\u201c analyze it.\n% In this scenario we could use interpolation to compute weather data at any\n% (longitude, latitude) location in the country. Assuming we have access to\n% weather archives we could also make historical comparisons. For example,\n% we could compute the total rainfall at (longitude, latitude)\n% for years 2008, 2007, 2006 etc, and likewise we could do the same for total\n% snowfall. If we were to use |griddata| to perform the interpolation, then\n% catchinggg catchinggg, we would pay the unnecessary overhead of computing\n% the underlying Delaunay triangulation each time a query is evaluated.\n%%\n% For example, suppose we have the following annual rainfall data for\n% 2006-2008:\n%%\n%       (xlon, ylat, annualRainfall06)\n%       (xlon, ylat, annualRainfall07)\n%       (xlon, ylat, annualRainfall08)\n%%\n% |xlon|, |ylat|, and |annualRainfall0*| are vectors that represent the sample\n% locations and the corresponding annual rainfall at these locations.\n% We can use |griddata| to query by interpolation the annual rainfall |zq| at\n% location |(xq, yq)|.\n%%\n%\n%       zq06 = griddata(xlon, ylat, annualRainfall06, xq, yq)\n%       zq07 = griddata(xlon, ylat, annualRainfall07, xq, yq)\n%       zq08 = griddata(xlon, ylat, annualRainfall08, xq, yq)\n%\n\n%%\n%\n% Note that a Delaunay triangulation of the points |(xlon, ylat)| is computed\n% each time the interpolation is evaluated at |(xq, yq)|. The same triangulation\n% could have been reused for the 07 and 08 queries, since we only changed\n% the values at each location. In this little illustrative example the penalty\n% may not be significant because the number of samples is likely to be in\n% the hundreds or less. However, when we interpolate large data sets this\n% represents a serious transaction penalty. |TriScatteredInterp|, the new\n% scattered data interpolation tool released in R2009a, allows us to avoid\n% this overhead. We can use it to construct the underlying triangulation once\n% and make repeated low-cost queries. In addition, for fixed sample locations\n% the sample values and the interpolation method can be changed independently\n% without re-triangulation. Repeating the previous illustration using\n% |TriScatteredInterp| we would get the following:\n%%\n%\n%       F = TriScatteredInterp(xlon, ylat, annualRainfall06)\n%       zq06 = F(xq, yq)\n%       F.V = annualRainfall07;\n%       zq07 = F(xq, yq)\n%       F.V = annualRainfall08;\n%       zq08 = F(xq, yq)\n%%\n% We can also change the interpolation method on-the-fly; this switches to\n% the new Natural Neighbor interpolation method.\n%\n%       F.Method = 'natural';\n%\n% Again, no transaction penalty for switching the method. To evaluate using\n% natural neighbor interpolation we do the following:\n%%\n%\n%       zq08 = F(xq, yq)\n%\n%%\n% As you can see, the interpolant works like a function handle making the\n% calling syntax more intuitive.\n%\n%% What's the Bottom Line?\n% The new scattered data interpolation function, |TriScatteredInterp|, allows\n% you to interpolate large data sets very efficiently. Because the\n% triangulation is cached within the black-box (the interpolant) you do not\n% need to pay the overhead of recomputing the triangulation each time you\n% evaluate, change the interpolation method, or change the values at the\n% sample locations.\n%%\n% What's more, as in the case of |DelaunayTri|,\n% this interpolation tool uses the\n% <http:\/\/www.cgal.org CGAL> Delaunay triangulations which are robust\n% against numerical problems. The triangulation algorithms are also faster\n% and more memory efficient. So if you use MATLAB to interpolate large\n% scattered data sets, this is great news for you.\n%%\n% The other bonus feature is the availability\n% of the Natural Neighbor interpolation method \u00e2\u20ac\u201c not be confused with Nearest\n% Neighbor interpolation. Natural neighbor interpolation is a\n% triangulation-based method that has an area-of-influence weighting\n% associated with each\n% sample point. It is widely used in the area of geostatistics because it\n% is C1 continuous (except at the sample locations) and it performs well in\n% both clustered and sparse data locations.\n%%\n% For more information on this\n% feature and the new computational geometry features shipped in R2009a\n% check out the overview\n% <https:\/\/www.mathworks.com\/products.htmldemos\/shipping\/matlab\/New-MATLAB-Mathematics-Features-in-R2009a.html |video|>, and Delaunay triangulation\n% <https:\/\/www.mathworks.com\/products.htmldemos\/shipping\/matlab\/demoDelaunayTri.html demo>,\n% or follow these links to the documentation\n%\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/trirep.html |TriRep|>,\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/delaunaytri.html |DelaunayTri|>,\n% * <https:\/\/www.mathworks.com\/help\/releases\/R2009a\/techdoc\/index.html?\/access\/helpdesk\/help\/releases\/R2009a\/techdoc\/ref\/triscatteredinterp.html |TriScatteredInterp|>.\n%% Feedback\n% If you use MATLAB to perform scattered data interpolation, let me know\n% what you think <https:\/\/blogs.mathworks.com\/loren\/?p193#respond here>.\n% It's the feedback from users like you that prompts us to\n% provide enhancements that are relevant to you.\n\n##### SOURCE END ##### 190f374fa98c4dc39fdcf0ac326de59d\n--><\/p>\n","protected":false},"excerpt":{"rendered":"<p>\nI am pleased to welcome back Damian Sheehy for the sequel to his earlier Computational Geometry post.<br \/>\n&nbsp;<br \/>\nContents<\/p>\n<p>What's your problem! ?<br \/>\ngriddata Produces \"Funny\" Results<br \/>\nScattered Data... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2009\/07\/22\/computational-geometry-in-matlab-r2009a-part-ii\/\">read more >><\/a><\/p>\n","protected":false},"author":39,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[32,6,11],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/193"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/users\/39"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/comments?post=193"}],"version-history":[{"count":4,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/193\/revisions"}],"predecessor-version":[{"id":4966,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/193\/revisions\/4966"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=193"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=193"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=193"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}