{"id":152,"date":"2008-08-25T08:35:15","date_gmt":"2008-08-25T13:35:15","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/2008\/08\/25\/piecewise-linear-interpolation\/"},"modified":"2016-08-01T10:27:53","modified_gmt":"2016-08-01T15:27:53","slug":"piecewise-linear-interpolation","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2008\/08\/25\/piecewise-linear-interpolation\/","title":{"rendered":"Piecewise Linear Interpolation"},"content":{"rendered":"<div class=\"content\">\n<tt>John D'Errico<\/tt> is back today to talk about linear interpolation.<\/p>\n<p>&nbsp;<\/p>\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li><a href=\"#1\">Introduction<\/a><\/li>\n<li><a href=\"#3\">Create Some Data to Interpolate<\/a><\/li>\n<li><a href=\"#5\">histc Solves the Binning Problem<\/a><\/li>\n<li><a href=\"#7\">Binning - A Loop With An Explicit Test<\/a><\/li>\n<li><a href=\"#8\">Binning - A Semi-vectorized Test<\/a><\/li>\n<li><a href=\"#11\">Fully Vectorized Binning<\/a><\/li>\n<li><a href=\"#13\">Interpolation as a Linear Combination<\/a><\/li>\n<li><a href=\"#14\">Do the Interpolation and Plot the Result<\/a><\/li>\n<li><a href=\"#15\">Use interp1 Instead<\/a><\/li>\n<\/ul>\n<\/div>\n<h3>Introduction<a name=\"1\"><\/a><\/h3>\n<p>You saw in my <a href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/06\/11\/interpolation-in-matlab\/\">previous blog<\/a> that high order polynomials can have some problems. Why not go to the opposite extreme? Use a piecewise version of <a href=\"http:\/\/en.wikipedia.org\/wiki\/Linear_interpolation\"><tt>linear interpolation<\/tt><\/a>? I like to call it <a href=\"http:\/\/en.wikipedia.org\/wiki\/Connect_the_dots\"><tt>connect-the-dots<\/tt><\/a>, after the child's game of that name. This is really the simplest interpolation of all.<\/p>\n<p>In MATLAB, given a list of points, sampled from some functional relationship in one dimension, how would we perform piecewise<br \/>\nlinear interpolation? There are really two steps.<\/p>\n<p>For any point u, given a set of (x,y) pairs with a monotonic vector x (by monotonic, I mean that x(k) &lt; x(k+1) ), first find<br \/>\nthe index k, such that<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_eq91533.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Second, perform the linear interpolation to predict the value of y at x=u, between the pair of points (x(k),y(k)) and (x(k+1),y(k+1)).<\/p>\n<p>Each data point in the list of points becomes a point where the slope of the piecewise linear interpolant changes to a new<br \/>\nvalue. However, the function is still continuous across those locations. So one might call these locations \"knots\" because<br \/>\nat those points consecutive polynomial segments are tied together. Knots is a common term applied to splines for these locations;<br \/>\nbreaks is a common alternative name.<\/p>\n<h3>Create Some Data to Interpolate<a name=\"3\"><\/a><\/h3>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">x = linspace(0,2*pi,10)';\r\ny = sin(x);\r\nplot(x,y,<span style=\"color: #a020f0;\">'o'<\/span>)\r\ntitle <span style=\"color: #a020f0;\">'A simple trig function'<\/span>\r\nxlabel <span style=\"color: #a020f0;\">X<\/span>\r\nylabel <span style=\"color: #a020f0;\">Y<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_01.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Suppose I want to interpolate this function at some intermediate point, perhaps u == 1.3?<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">u = 1.3;<\/pre>\n<h3>histc Solves the Binning Problem<a name=\"5\"><\/a><\/h3>\n<p>The first step I describe above is what I call binning. <a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/histc.html\"><tt>histc<\/tt><\/a> solves that problem efficiently.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">[k,k] = histc(u,x);\r\nk<\/pre>\n<pre style=\"font-style: oblique;\">k =\r\n     2\r\n<\/pre>\n<p>As an aside, look at the way I took the output from <tt>histc<\/tt>. Since I only need the second output from <tt>histc<\/tt> but not the first output, rather than clutter up my workspace with a variable that I did not need, the output style <tt>[k,k]<\/tt> returns only the information I need.<\/p>\n<p>Next, it seems instructive to dive a little more deeply into binning, so let me offer a few alternatives to <tt>histc<\/tt>.<\/p>\n<h3>Binning - A Loop With An Explicit Test<a name=\"7\"><\/a><\/h3>\n<p>Just a test inside a loop suffices.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\"><span style=\"color: #0000ff;\">for<\/span> k = 1:(length(x)-1)\r\n  <span style=\"color: #0000ff;\">if<\/span> (x(k) &lt;= u) &amp;&amp; (u &lt; x(k+1))\r\n    <span style=\"color: #0000ff;\">break<\/span>\r\n  <span style=\"color: #0000ff;\">end<\/span>\r\n<span style=\"color: #0000ff;\">end<\/span>\r\nx\r\nk<\/pre>\n<pre style=\"font-style: oblique;\">x =\r\n            0\r\n      0.69813\r\n       1.3963\r\n       2.0944\r\n       2.7925\r\n       3.4907\r\n       4.1888\r\n       4.8869\r\n       5.5851\r\n       6.2832\r\nk =\r\n     2\r\n<\/pre>\n<h3>Binning - A Semi-vectorized Test<a name=\"8\"><\/a><\/h3>\n<p>Do the binning with a single vectorized test. This works reasonably as long as I have only a scalar value for <tt>u<\/tt>, so I call this only a semi-vectorized solution.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">k = sum(u&gt;=x)<\/pre>\n<pre style=\"font-style: oblique;\">k =\r\n     2\r\n<\/pre>\n<p>When I want to place multiple points into their respective bins, this <tt>sum<\/tt> fails.<\/p>\n<p>There are other ways to place points in bins in MATLAB, including a <tt>sort<\/tt>, or with hash tables. You can find several different binning methods in <tt>bindex<\/tt> on the file exchange. It is useful mainly to those with older MATLAB releases, because <tt>histc became<\/tt> available with version 5.3 and later of MATLAB.<\/p>\n<h3>Fully Vectorized Binning<a name=\"11\"><\/a><\/h3>\n<p>Next, I may, and often do, have a list of points to interpolate. This is a common event, where I wish to more finely resample<br \/>\na curve that is sampled only at some short list of points. This time, let me generate 1000 points at which to interpolate<br \/>\nthe sampled function.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">u = linspace(x(1),x(end),1000)';\r\n\r\n[k,k] = histc(u,x);<\/pre>\n<p>This next line handles the very last point. Recall the definition of a bin as<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_eq91533.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Note the strict inequality on the left. So that very last point (at <tt>x(end)<\/tt>) might be technically said to fall into the <tt>n|th bin when I have |n<\/tt> break points.<\/p>\n<p>On the other hand, it is more convenient to put anything that lies exactly at the very last break point into bin <tt>(n-1)<\/tt>.<\/p>\n<p>By the way, any piecewise interpolation should worry about points that fall fully above or below the domain of the data. Will<br \/>\nthe code extrapolate or not? Should you extrapolate?<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">n = length(x);\r\nk(k == n) = n - 1;<\/pre>\n<h3>Interpolation as a Linear Combination<a name=\"13\"><\/a><\/h3>\n<p>The final step is to interpolate between two points. Long ago, I recall from high school what was called a point-slope form<br \/>\nfor a line. If you know a pair of points that a line passes through, as <tt>(x(k),y(k))<\/tt> and <tt>(x(k+1),y(k+1))<\/tt>, then the slope of the line is simple to compute. An equation for our line as a function of the parameter <tt>u<\/tt> is just:<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_eq94249.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>I can also rewrite this in a way that I like, by defining a parameter <tt>t<\/tt> as<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_eq05093.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>Then the interpolant is a simple linear combination of the function values at each end of the interval.<\/p>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_eq98989.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<h3>Do the Interpolation and Plot the Result<a name=\"14\"><\/a><\/h3>\n<p>See how nicely this all worked, in a fully vectorized coding style.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">t = (u - x(k)).\/(x(k+1) - x(k));\r\nyu = (1-t).*y(k) + t.*y(k+1);\r\n\r\nplot(x,y,<span style=\"color: #a020f0;\">'bo'<\/span>,u,yu,<span style=\"color: #a020f0;\">'r-'<\/span>)\r\nxlabel <span style=\"color: #a020f0;\">X<\/span>\r\nylabel <span style=\"color: #a020f0;\">Y<\/span>\r\ntitle <span style=\"color: #a020f0;\">'The connect-the-dots interpolant'<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_02.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<h3>Use interp1 Instead<a name=\"15\"><\/a><\/h3>\n<p>It is always better to use a built-in tool to solve your problem than to do it yourself, so I might just as well have used<br \/>\n<a href=\"https:\/\/www.mathworks.com\/help\/matlab\/ref\/interp1.html\"><tt>interp1<\/tt><\/a> to accomplish this interpolation.<\/p>\n<pre style=\"background: #F9F7F3; padding: 10px; border: 1px solid #c8c8c8;\">yinterp1 = interp1(x,y,u,<span style=\"color: #a020f0;\">'linear'<\/span>);\r\n\r\nplot(x,y,<span style=\"color: #a020f0;\">'bo'<\/span>,u,yinterp1,<span style=\"color: #a020f0;\">'r-'<\/span>)\r\nxlabel <span style=\"color: #a020f0;\">X<\/span>\r\nylabel <span style=\"color: #a020f0;\">Y<\/span>\r\ntitle <span style=\"color: #a020f0;\">'The connect-the-dots interpolant, using interp1'<\/span><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/loren\/152\/interpJD3_03.png\" hspace=\"5\" vspace=\"5\" \/><\/p>\n<p>In my next blog I'll begin talking about piecewise cubic interpolants. Until then please tell me your ideas <a href=\"https:\/\/blogs.mathworks.com\/loren\/?p=152#respond\"><tt>here<\/tt><\/a>. Are there some associated topics that I should cover?<\/p>\n<p><script>\/\/ <![CDATA[\nfunction grabCode_4ee4189285944b34aebd605b231ff953() {\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='4ee4189285944b34aebd605b231ff953 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 4ee4189285944b34aebd605b231ff953';\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 2008 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.6<\/p>\n<\/div>\n<p><!--\n4ee4189285944b34aebd605b231ff953 ##### SOURCE BEGIN #####\n%% Piecewise Linear Interpolation\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadAuthor.do?objectId=801347&objectType=author |John D'Errico|>\n% is back today to talk about linear interpolation.\n\n%% Introduction\n% You saw in my\n% <https:\/\/blogs.mathworks.com\/loren\/2008\/06\/11\/interpolation-in-matlab\/ previous blog>\n% that high order polynomials can have some problems.\n% Why not go to the opposite extreme? Use a piecewise version of\n% <http:\/\/en.wikipedia.org\/wiki\/Linear_interpolation |linear interpolation|>? I like to call it\n% <http:\/\/en.wikipedia.org\/wiki\/Connect_the_dots |connect-the-dots|>,\n% after the child's game of that name. This is really the simplest\n% interpolation of all.\n%\n% In MATLAB, given a list of points, sampled from some functional\n% relationship in one dimension, how would we perform piecewise linear interpolation?\n% There are really two steps.\n\n%%\n% For any point u, given a set of (x,y) pairs\n% with a monotonic vector x (by monotonic, I mean that x(k) < x(k+1) ),\n% first find the index k, such that\n%\n% $$x(k) <= u < x(k+1)$$\n%\n% Second, perform the linear interpolation to predict the value of y\n% at x=u, between the pair of points (x(k),y(k)) and (x(k+1),y(k+1)).\n%\n% Each data point in the list of points becomes a point where the slope\n% of the piecewise linear interpolant changes to a new value. However,\n% the function is still continuous across those locations. So one\n% might call these locations \"knots\" because at those points consecutive polynomial\n% segments are tied together. Knots is a common term applied to splines for\n% these locations; breaks is a common alternative name.\n\n%% Create Some Data to Interpolate\nx = linspace(0,2*pi,10)';\ny = sin(x);\nplot(x,y,'o')\ntitle 'A simple trig function'\nxlabel X\nylabel Y\n\n%%\n% Suppose I want to interpolate this function at some intermediate point,\n% perhaps u == 1.3?\nu = 1.3;\n\n%% histc Solves the Binning Problem\n% The first step I describe above is what I call binning.\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/histc.html |histc|>\n% solves that problem efficiently.\n[k,k] = histc(u,x);\nk\n\n%%\n% As an aside, look at the way I took the output from |histc|. Since I only\n% need the second output from |histc| but not the first output, rather than\n% clutter up my workspace\n% with a variable that I did not need, the output style |[k,k]|\n% returns only the information I need.\n%\n% Next, it seems instructive to dive a little more deeply into binning, so let me\n% offer a few alternatives to |histc|.\n\n%% Binning - A Loop With An Explicit Test\n% Just a test inside a loop suffices.\nfor k = 1:(length(x)-1)\nif (x(k) <= u) && (u < x(k+1)) break end end x k %% Binning - A Semi-vectorized Test % Do the binning with a single vectorized test. This works reasonably as % long as I have only a scalar value for |u|, so I call this only a % semi-vectorized solution. k = sum(u>=x)\n\n%%\n% When I want to place multiple points into their respective bins, this |sum| fails.\n\n%%\n% There are other ways to place points in bins in MATLAB, including\n% a |sort|, or with hash tables. You can find several different binning\n% methods in\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/loadFile.do?objectId=11&objectType=FILE |bindex|>\n% on the file exchange. It is useful mainly to those with older MATLAB\n% releases, because |histc became| available with version\n% 5.3 and later of MATLAB.\n\n%% Fully Vectorized Binning\n% Next, I may, and often do, have a list of points to interpolate. This is a\n% common event, where I wish to more finely resample a curve that is sampled\n% only at some short list of points. This time, let me generate 1000\n% points at which to interpolate the sampled function.\nu = linspace(x(1),x(end),1000)';\n\n[k,k] = histc(u,x);\n%%\n% This next line handles the very last point. Recall the definition of a\n% bin as\n%\n% $$x(k) <= u < x(k+1)$$\n%\n% Note the strict inequality on the left. So that very last point (at |x(end)|)\n% might be technically said to fall into the |n|th bin when I have\n% |n| break points.\n%\n% On the other hand, it is more convenient to put anything\n% that lies exactly at the very last break point into bin |(n-1)|.\n%\n% By the way, any piecewise interpolation\n% should worry about points that fall fully above or below the domain of the data.\n% Will the code extrapolate or not? Should you extrapolate?\nn = length(x);\nk(k == n) = n - 1;\n\n%% Interpolation as a Linear Combination\n% The final step is to interpolate between two points. Long ago, I recall\n% from high school what was called a point-slope form for a line. If you know a pair\n% of points that a line passes through, as |(x(k),y(k))| and |(x(k+1),y(k+1))|,\n% then the slope of the line is simple to compute.\n% An equation for our line as a function of the parameter |u| is just:\n%\n% $$y(u) = y_k + \\frac{y_{k+1} - y_k}{x_{k+1} - x_k}(u - x_k)$$\n%\n% I can also rewrite this in a way that I like, by defining a parameter\n% |t| as\n%\n% $$t_u = \\frac{u - x_k}{x_{k+1} - x_k}$$\n%\n% Then the interpolant is a simple linear combination of the\n% function values at each end of the interval.\n%\n% $$y(u) = (1-t_u)y_k + t_u y_{k+1}$$\n\n%% Do the Interpolation and Plot the Result\n% See how nicely this all worked, in a fully vectorized coding style.\nt = (u - x(k)).\/(x(k+1) - x(k));\nyu = (1-t).*y(k) + t.*y(k+1);\n\nplot(x,y,'bo',u,yu,'r-')\nxlabel X\nylabel Y\ntitle 'The connect-the-dots interpolant'\n\n%% Use interp1 Instead\n% It is always better to use a built-in tool to solve your problem than\n% to do it yourself, so I might just as well have used\n% <https:\/\/www.mathworks.com\/help\/matlab\/ref\/interp1.html |interp1|>\n% to accomplish this interpolation.\nyinterp1 = interp1(x,y,u,'linear');\n\nplot(x,y,'bo',u,yinterp1,'r-')\nxlabel X\nylabel Y\ntitle 'The connect-the-dots interpolant, using interp1'\n\n%%\n% In my next blog\n% I'll begin talking about piecewise cubic interpolants. Until then\n% please tell me your ideas\n% <https:\/\/blogs.mathworks.com\/loren\/?p=152#respond |here|>.\n% Are there some associated topics that I should cover?\n##### SOURCE END ##### 4ee4189285944b34aebd605b231ff953\n--><\/p>\n","protected":false},"excerpt":{"rendered":"<p>\nJohn D'Errico is back today to talk about linear interpolation.<br \/>\n&nbsp;<br \/>\nContents<\/p>\n<p>Introduction<br \/>\nCreate Some Data to Interpolate<br \/>\nhistc Solves the Binning Problem<br \/>\nBinning - A Loop With An Explicit... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/loren\/2008\/08\/25\/piecewise-linear-interpolation\/\">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":[23],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/152"}],"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=152"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/152\/revisions"}],"predecessor-version":[{"id":1869,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/posts\/152\/revisions\/1869"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/media?parent=152"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/categories?post=152"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/loren\/wp-json\/wp\/v2\/tags?post=152"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}