Boost C++ Libraries: Ticket #11930: huiller method is inaccurate

Boost C++ Libraries: Ticket #11930: huiller method is inaccurate https://svn.boost.org/trac10/ticket/11930 <p> The following code illustrates a problem with the "huiller" strategy for computing spherical areas. It computes the area of a sequence of progressively skinny triangles. The relative errors in the computed areas are much larger than they need be. </p> <div class="wiki-code"><div class="code"><pre><span class="cp">#include</span> <span class="cpf"><iostream></span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf"><cmath></span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf"><boost/geometry.hpp></span><span class="cp"></span> <span class="k">namespace</span> <span class="n">bg</span> <span class="o">=</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="p">;</span> <span class="k">typedef</span> <span class="n">bg</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">point</span><span class="o"><</span><span class="kt">double</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">bg</span><span class="o">::</span><span class="n">cs</span><span class="o">::</span><span class="n">spherical_equatorial</span><span class="o"><</span><span class="n">bg</span><span class="o">::</span><span class="n">degree</span><span class="o">></span> <span class="o">></span> <span class="n">pt</span><span class="p">;</span> <span class="kt">int</span> <span class="nf">main</span><span class="p">()</span> <span class="p">{</span> <span class="n">bg</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">polygon</span><span class="o"><</span><span class="n">pt</span><span class="p">,</span> <span class="nb">false</span><span class="p">,</span> <span class="nb">false</span><span class="o">></span> <span class="n">poly</span><span class="p">;</span> <span class="n">bg</span><span class="o">::</span><span class="n">read_wkt</span><span class="p">(</span><span class="s">"POLYGON((0 45,0 -45,0 0))"</span><span class="p">,</span> <span class="n">poly</span><span class="p">);</span> <span class="kt">double</span> <span class="n">c</span> <span class="o">=</span> <span class="mf">0.014458780939650811</span><span class="p">;</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o"><<</span> <span class="s">" area relerr"</span> <span class="o"><<</span> <span class="s">"</span><span class="se">\n</span><span class="s">"</span><span class="p">;</span> <span class="k">for</span> <span class="p">(</span><span class="kt">int</span> <span class="n">i</span> <span class="o">=</span> <span class="mi">3</span><span class="p">;</span> <span class="n">i</span> <span class="o"><=</span> <span class="mi">6</span><span class="p">;</span> <span class="o">++</span><span class="n">i</span><span class="p">)</span> <span class="p">{</span> <span class="kt">double</span> <span class="n">h</span> <span class="o">=</span> <span class="n">std</span><span class="o">::</span><span class="n">pow</span><span class="p">(</span><span class="mf">10.0</span><span class="p">,</span> <span class="o">-</span><span class="n">i</span><span class="p">);</span> <span class="n">poly</span><span class="p">.</span><span class="n">outer</span><span class="p">()[</span><span class="mi">2</span><span class="p">].</span><span class="n">set</span><span class="o"><</span><span class="mi">0</span><span class="o">></span><span class="p">(</span><span class="n">h</span><span class="p">);</span> <span class="kt">double</span> <span class="n">area</span> <span class="o">=</span> <span class="n">bg</span><span class="o">::</span><span class="n">area</span><span class="p">(</span><span class="n">poly</span><span class="p">),</span> <span class="c1">// Relative error in this estimate is 2e-11 for h = 0.001</span> <span class="n">area0</span> <span class="o">=</span> <span class="n">h</span> <span class="o">*</span> <span class="n">c</span><span class="p">;</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o"><<</span> <span class="n">area</span> <span class="o"><<</span> <span class="s">" "</span> <span class="o"><<</span> <span class="p">(</span><span class="n">area</span> <span class="o">-</span> <span class="n">area0</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">area0</span><span class="p">)</span> <span class="o"><<</span> <span class="s">"</span><span class="se">\n</span><span class="s">"</span><span class="p">;</span> <span class="p">}</span> <span class="p">}</span> </pre></div></div><p> Running this code gives </p> <pre class="wiki"> area relerr 1.44588e-05 -2.0871e-06 1.44594e-06 4.55941e-05 1.43961e-07 -0.00433739 0 -1 </pre><p> Discussion: </p> <ul><li>Surely it's a bad idea to embed the algorithm name "huiller" into the API. Computing spherical areas is really straightforward. Users shouldn't need to worry about the underlying method. Maybe "spherical_area" is a better name. </li><li>In any case, you've misspelled the guy's name. It's either "l'Huilier" or "l'Huillier" with two "i"s. </li><li>Furthermore, you're needlessly blackening l'Huilier's name by applying his formula for the area of a spherical triangle in terms of its sides. </li><li>The three-side formula is a terrible choice because if a + b is nearly equal to c (as is the case in my example), the triangle is badly condititioned. </li><li>Much better is to use the formula of the area of a triangle in terms of two sides and the included angle. See the last two formulas in the Wikipedia article on <a class="ext-link" href="https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess"><span class="icon"></span>spherical excess</a>. This nicely reduces to the familar trapezium rule in the plane limit. </li><li>This formula comes with a sign attached, so there's no need to kludge a sign on afterwards as there is with the current implementation. </li><li>This formula together with an ellipsoidal correction is used by <a class="ext-link" href="http://geographiclib.sourceforge.net/html/"><span class="icon"></span>GeographicLib</a> for computing geodesic areas. </li></ul> en-us Boost C++ Libraries /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/11930 Trac 1.4.3 awulkiew Wed, 29 Mar 2017 22:46:32 GMT status, milestone changed; resolution set https://svn.boost.org/trac10/ticket/11930#comment:1 https://svn.boost.org/trac10/ticket/11930#comment:1 <ul> <li><strong>status</strong> <span class="trac-field-old">new</span> → <span class="trac-field-new">closed</span> </li> <li><strong>resolution</strong> → <span class="trac-field-new">fixed</span> </li> <li><strong>milestone</strong> <span class="trac-field-old">To Be Determined</span> → <span class="trac-field-new">Boost 1.64.0</span> </li> </ul> <p> Thanks! </p> <p> Fixed by this PR: <a class="ext-link" href="https://github.com/boostorg/geometry/pull/359"><span class="icon"></span>https://github.com/boostorg/geometry/pull/359</a> </p> <p> The name is changed to <code>area::spherical</code>. </p> <p> The new results are: </p> <pre class="wiki"> area relerr 1.44588e-05 2.10295e-11 1.44588e-06 2.10165e-13 1.44588e-07 2.01378e-15 1.44588e-08 0 </pre> Ticket