[ceph.git] / ceph / src / boost / boost / geometry / srs / projections / impl / pj_apply_gridshift.hpp

// Boost.Geometry
// This file is manually converted from PROJ4

// This file was modified by Oracle on 2018, 2019.
// Modifications copyright (c) 2018-2019, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

// This file is converted from PROJ4, http://trac.osgeo.org/proj
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
// PROJ4 is maintained by Frank Warmerdam
// This file was converted to Geometry Library by Adam Wulkiewicz

// Original copyright notice:
// Author:   Frank Warmerdam, warmerdam@pobox.com

// Copyright (c) 2000, Frank Warmerdam

// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.

#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP


#include <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/radian_access.hpp>

#include <boost/geometry/srs/projections/impl/adjlon.hpp>
#include <boost/geometry/srs/projections/impl/function_overloads.hpp>
#include <boost/geometry/srs/projections/impl/pj_gridlist.hpp>

#include <boost/geometry/util/range.hpp>


namespace boost { namespace geometry { namespace projections
{

namespace detail
{

// Originally implemented in nad_intr.c
template <typename CalcT>
inline void nad_intr(CalcT in_lon, CalcT in_lat,
                     CalcT & out_lon, CalcT & out_lat,
                     pj_ctable const& ct)
{
	pj_ctable::lp_t frct;
	pj_ctable::ilp_t indx;
	boost::int32_t in;

	indx.lam = int_floor(in_lon /= ct.del.lam);
	indx.phi = int_floor(in_lat /= ct.del.phi);
	frct.lam = in_lon - indx.lam;
	frct.phi = in_lat - indx.phi;
    // TODO: implement differently
	out_lon = out_lat = HUGE_VAL;
	if (indx.lam < 0) {
		if (indx.lam == -1 && frct.lam > 0.99999999999) {
			++indx.lam;
			frct.lam = 0.;
		} else
			return;
	} else if ((in = indx.lam + 1) >= ct.lim.lam) {
		if (in == ct.lim.lam && frct.lam < 1e-11) {
			--indx.lam;
			frct.lam = 1.;
		} else
			return;
	}
	if (indx.phi < 0) {
		if (indx.phi == -1 && frct.phi > 0.99999999999) {
			++indx.phi;
			frct.phi = 0.;
		} else
			return;
	} else if ((in = indx.phi + 1) >= ct.lim.phi) {
		if (in == ct.lim.phi && frct.phi < 1e-11) {
			--indx.phi;
			frct.phi = 1.;
		} else
			return;
	}
	boost::int32_t index = indx.phi * ct.lim.lam + indx.lam;
	pj_ctable::flp_t const& f00 = ct.cvs[index++];
	pj_ctable::flp_t const& f10 = ct.cvs[index];
	index += ct.lim.lam;
	pj_ctable::flp_t const& f11 = ct.cvs[index--];
	pj_ctable::flp_t const& f01 = ct.cvs[index];
    CalcT m00, m10, m01, m11;
	m11 = m10 = frct.lam;
	m00 = m01 = 1. - frct.lam;
	m11 *= frct.phi;
	m01 *= frct.phi;
	frct.phi = 1. - frct.phi;
	m00 *= frct.phi;
	m10 *= frct.phi;
	out_lon = m00 * f00.lam + m10 * f10.lam +
			  m01 * f01.lam + m11 * f11.lam;
	out_lat = m00 * f00.phi + m10 * f10.phi +
			  m01 * f01.phi + m11 * f11.phi;
}

// Originally implemented in nad_cvt.c
template <bool Inverse, typename CalcT>
inline void nad_cvt(CalcT const& in_lon, CalcT const& in_lat,
                    CalcT & out_lon, CalcT & out_lat,
                    pj_gi const& gi)
{
    static const int max_iterations = 10;
    static const CalcT tol = 1e-12;
    static const CalcT toltol = tol * tol;
    static const CalcT pi = math::pi<CalcT>();

    // horizontal grid expected
    BOOST_GEOMETRY_ASSERT_MSG(gi.format != pj_gi::gtx,
        "Vertical grid cannot be used in horizontal shift.");

    pj_ctable const& ct = gi.ct;

    // TODO: implement differently
    if (in_lon == HUGE_VAL)
    {
        out_lon = HUGE_VAL;
        out_lat = HUGE_VAL;
        return;
    }

    // normalize input to ll origin
    pj_ctable::lp_t tb;
    tb.lam = in_lon - ct.ll.lam;
    tb.phi = in_lat - ct.ll.phi;
    tb.lam = adjlon (tb.lam - pi) + pi;

    pj_ctable::lp_t t;
    nad_intr(tb.lam, tb.phi, t.lam, t.phi, ct);
    if (t.lam == HUGE_VAL)
    {
        out_lon = HUGE_VAL;
        out_lat = HUGE_VAL;
        return;
    }

    if (! Inverse)
    {
        out_lon = in_lon - t.lam;
        out_lat = in_lat - t.phi;
        return;
    }

    t.lam = tb.lam + t.lam;
    t.phi = tb.phi - t.phi;

    int i = max_iterations;
    pj_ctable::lp_t del, dif;
    do
    {
        nad_intr(t.lam, t.phi, del.lam, del.phi, ct);

        // This case used to return failure, but I have
        // changed it to return the first order approximation
        // of the inverse shift.  This avoids cases where the
        // grid shift *into* this grid came from another grid.
        // While we aren't returning optimally correct results
        // I feel a close result in this case is better than
        // no result.  NFW
        // To demonstrate use -112.5839956 49.4914451 against
        // the NTv2 grid shift file from Canada.
        if (del.lam == HUGE_VAL)
        {
            // Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.
            break;
        }

        dif.lam = t.lam - del.lam - tb.lam;
        dif.phi = t.phi + del.phi - tb.phi;
        t.lam -= dif.lam;
        t.phi -= dif.phi;

    }
    while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); // prob. slightly faster than hypot()

    if (i==0)
    {
        // Inverse grid shift iterator failed to converge.
        out_lon = HUGE_VAL;
        out_lat = HUGE_VAL;
        return;
    }

    out_lon = adjlon (t.lam + ct.ll.lam);
    out_lat = t.phi + ct.ll.phi;
}


/************************************************************************/
/*                             find_grid()                              */
/*                                                                      */
/*    Determine which grid is the correct given an input coordinate.    */
/************************************************************************/

// Originally find_ctable()
// here divided into grid_disjoint(), find_grid() and load_grid()

template <typename T>
inline bool grid_disjoint(T const& lam, T const& phi,
                          pj_ctable const& ct)
{
    double epsilon = (fabs(ct.del.phi)+fabs(ct.del.lam))/10000.0;
    return ct.ll.phi - epsilon > phi
        || ct.ll.lam - epsilon > lam
        || (ct.ll.phi + (ct.lim.phi-1) * ct.del.phi + epsilon < phi)
        || (ct.ll.lam + (ct.lim.lam-1) * ct.del.lam + epsilon < lam);
}

template <typename T>
inline pj_gi * find_grid(T const& lam,
                         T const& phi,
                         std::vector<pj_gi>::iterator first,
                         std::vector<pj_gi>::iterator last)
{
    pj_gi * gip = NULL;

    for( ; first != last ; ++first )
    {
        // skip tables that don't match our point at all.
        if (! grid_disjoint(lam, phi, first->ct))
        {
            // skip vertical grids
            if (first->format != pj_gi::gtx)
            {
                gip = boost::addressof(*first);
                break;
            }
        }
    }

    // If we didn't find a child then nothing more to do
    if( gip == NULL )
        return gip;

    // Otherwise use the child, first checking it's children
    pj_gi * child = find_grid(lam, phi, first->children.begin(), first->children.end());
    if (child != NULL)
        gip = child;

    return gip;
}

template <typename T>
inline pj_gi * find_grid(T const& lam,
                         T const& phi,
                         pj_gridinfo & grids,
                         std::vector<std::size_t> const& gridindexes)
{
    pj_gi * gip = NULL;

    // keep trying till we find a table that works
    for (std::size_t i = 0 ; i < gridindexes.size() ; ++i)
    {
        pj_gi & gi = grids[gridindexes[i]];

        // skip tables that don't match our point at all.
        if (! grid_disjoint(lam, phi, gi.ct))
        {
            // skip vertical grids
            if (gi.format != pj_gi::gtx)
            {
                gip = boost::addressof(gi);
                break;
            }
        }
    }

    if (gip == NULL)
        return gip;

    // If we have child nodes, check to see if any of them apply.
    pj_gi * child = find_grid(lam, phi, gip->children.begin(), gip->children.end());
    if (child != NULL)
        gip = child;

    // if we get this far we have found a suitable grid
    return gip;
}


template <typename StreamPolicy>
inline bool load_grid(StreamPolicy const& stream_policy, pj_gi_load & gi)
{
    // load the grid shift info if we don't have it.
    if (gi.ct.cvs.empty())
    {
        typename StreamPolicy::stream_type is;
        stream_policy.open(is, gi.gridname);

        if (! pj_gridinfo_load(is, gi))
        {
            //pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
            return false;
        }
    }

    return true;
}


/************************************************************************/
/*                        pj_apply_gridshift_3()                        */
/*                                                                      */
/*      This is the real workhorse, given a gridlist.                   */
/************************************************************************/

// Generic stream policy and standard grids
template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename Grids>
inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
                                 Range & range,
                                 Grids & grids,
                                 std::vector<std::size_t> const& gridindexes,
                                 grids_tag)
{
    typedef typename boost::range_size<Range>::type size_type;
    
    // If the grids are empty the indexes are as well
    if (gridindexes.empty())
    {
        //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
        //return PJD_ERR_FAILED_TO_LOAD_GRID;
        return false;
    }

    size_type point_count = boost::size(range);

    for (size_type i = 0 ; i < point_count ; ++i)
    {
        typename boost::range_reference<Range>::type
            point = range::at(range, i);

        CalcT in_lon = geometry::get_as_radian<0>(point);
        CalcT in_lat = geometry::get_as_radian<1>(point);
        
        pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes);

        if ( gip != NULL )
        {
            // load the grid shift info if we don't have it.
            if (! gip->ct.cvs.empty() || load_grid(stream_policy, *gip))
            {
                // TODO: use set_invalid_point() or similar mechanism
                CalcT out_lon = HUGE_VAL;
                CalcT out_lat = HUGE_VAL;

                nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
            
                // TODO: check differently
                if ( out_lon != HUGE_VAL )
                {
                    geometry::set_from_radian<0>(point, out_lon);
                    geometry::set_from_radian<1>(point, out_lat);
                }
            }
        }
    }

    return true;
}

// Generic stream policy and shared grids
template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename SharedGrids>
inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
                                 Range & range,
                                 SharedGrids & grids,
                                 std::vector<std::size_t> const& gridindexes,
                                 shared_grids_tag)
{
    typedef typename boost::range_size<Range>::type size_type;
    
    // If the grids are empty the indexes are as well
    if (gridindexes.empty())
    {
        //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
        //return PJD_ERR_FAILED_TO_LOAD_GRID;
        return false;
    }

    size_type point_count = boost::size(range);

    // local storage
    pj_gi_load local_gi;

    for (size_type i = 0 ; i < point_count ; )
    {
        bool load_needed = false;

        CalcT in_lon = 0;
        CalcT in_lat = 0;

        {
            typename SharedGrids::read_locked lck_grids(grids);

            for ( ; i < point_count ; ++i )
            {
                typename boost::range_reference<Range>::type
                    point = range::at(range, i);

                in_lon = geometry::get_as_radian<0>(point);
                in_lat = geometry::get_as_radian<1>(point);

                pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);

                if (gip == NULL)
                {
                    // do nothing
                }
                else if (! gip->ct.cvs.empty())
                {
                    // TODO: use set_invalid_point() or similar mechanism
                    CalcT out_lon = HUGE_VAL;
                    CalcT out_lat = HUGE_VAL;

                    nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);

                    // TODO: check differently
                    if (out_lon != HUGE_VAL)
                    {
                        geometry::set_from_radian<0>(point, out_lon);
                        geometry::set_from_radian<1>(point, out_lat);
                    }
                }
                else
                {
                    // loading is needed
                    local_gi = *gip;
                    load_needed = true;
                    break;
                }
            }
        }

        if (load_needed)
        {
            if (load_grid(stream_policy, local_gi))
            {
                typename SharedGrids::write_locked lck_grids(grids);

                // check again in case other thread already loaded the grid.
                pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);

                if (gip != NULL && gip->ct.cvs.empty())
                {
                    // swap loaded local storage with empty grid
                    local_gi.swap(*gip);
                }
            }
            else
            {
                ++i;
            }
        }
    }

    return true;
}


/************************************************************************/
/*                        pj_apply_gridshift_2()                        */
/*                                                                      */
/*      This implementation uses the gridlist from a coordinate         */
/*      system definition.  If the gridlist has not yet been            */
/*      populated in the coordinate system definition we set it up      */
/*      now.                                                            */
/************************************************************************/

template <bool Inverse, typename Par, typename Range, typename ProjGrids>
inline bool pj_apply_gridshift_2(Par const& /*defn*/, Range & range, ProjGrids const& proj_grids)
{
    /*if( defn->catalog_name != NULL )
        return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
                                      x, y, z );*/

    /*std::vector<std::size_t> gridindexes;
    pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"),
                              grids.storage_ptr->stream_policy,
                              grids.storage_ptr->grids,
                              gridindexes);*/

    // At this point the grids should be initialized
    if (proj_grids.hindexes.empty())
        return false;

    return pj_apply_gridshift_3
            <
                Inverse, typename Par::type
            >(proj_grids.grids_storage().stream_policy,
              range,
              proj_grids.grids_storage().hgrids,
              proj_grids.hindexes,
              typename ProjGrids::grids_storage_type::grids_type::tag());
}

template <bool Inverse, typename Par, typename Range>
inline bool pj_apply_gridshift_2(Par const& , Range & , srs::detail::empty_projection_grids const& )
{
    return false;
}


} // namespace detail

}}} // namespace boost::geometry::projections

#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
Commit	Line	Data
92f5a8d4 TL	1	// Boost.Geometry
	2	// This file is manually converted from PROJ4
	3
	4	// This file was modified by Oracle on 2018, 2019.
	5	// Modifications copyright (c) 2018-2019, Oracle and/or its affiliates.
	6	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	7
	8	// Use, modification and distribution is subject to the Boost Software License,
	9	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	10	// http://www.boost.org/LICENSE_1_0.txt)
	11
	12	// This file is converted from PROJ4, http://trac.osgeo.org/proj
	13	// PROJ4 is originally written by Gerald Evenden (then of the USGS)
	14	// PROJ4 is maintained by Frank Warmerdam
	15	// This file was converted to Geometry Library by Adam Wulkiewicz
	16
	17	// Original copyright notice:
	18	// Author: Frank Warmerdam, warmerdam@pobox.com
	19
	20	// Copyright (c) 2000, Frank Warmerdam
	21
	22	// Permission is hereby granted, free of charge, to any person obtaining a
	23	// copy of this software and associated documentation files (the "Software"),
	24	// to deal in the Software without restriction, including without limitation
	25	// the rights to use, copy, modify, merge, publish, distribute, sublicense,
	26	// and/or sell copies of the Software, and to permit persons to whom the
	27	// Software is furnished to do so, subject to the following conditions:
	28
	29	// The above copyright notice and this permission notice shall be included
	30	// in all copies or substantial portions of the Software.
	31
	32	// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
	33	// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	34	// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	35	// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	36	// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
	37	// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
	38	// DEALINGS IN THE SOFTWARE.
	39
	40	#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
	41	#define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
	42
	43
	44	#include <boost/geometry/core/assert.hpp>
	45	#include <boost/geometry/core/radian_access.hpp>
	46
	47	#include <boost/geometry/srs/projections/impl/adjlon.hpp>
	48	#include <boost/geometry/srs/projections/impl/function_overloads.hpp>
	49	#include <boost/geometry/srs/projections/impl/pj_gridlist.hpp>
	50
	51	#include <boost/geometry/util/range.hpp>
	52
	53
	54	namespace boost { namespace geometry { namespace projections
	55	{
	56
	57	namespace detail
	58	{
	59
	60	// Originally implemented in nad_intr.c
	61	template <typename CalcT>
	62	inline void nad_intr(CalcT in_lon, CalcT in_lat,
	63	CalcT & out_lon, CalcT & out_lat,
	64	pj_ctable const& ct)
65	{
66	pj_ctable::lp_t frct;
67	pj_ctable::ilp_t indx;
68	boost::int32_t in;
69
70	indx.lam = int_floor(in_lon /= ct.del.lam);
71	indx.phi = int_floor(in_lat /= ct.del.phi);
72	frct.lam = in_lon - indx.lam;
73	frct.phi = in_lat - indx.phi;
74	// TODO: implement differently
75	out_lon = out_lat = HUGE_VAL;
76	if (indx.lam < 0) {
77	if (indx.lam == -1 && frct.lam > 0.99999999999) {
78	++indx.lam;
79	frct.lam = 0.;
80	} else
81	return;
82	} else if ((in = indx.lam + 1) >= ct.lim.lam) {
83	if (in == ct.lim.lam && frct.lam < 1e-11) {
84	--indx.lam;
85	frct.lam = 1.;
86	} else
87	return;
88	}
89	if (indx.phi < 0) {
90	if (indx.phi == -1 && frct.phi > 0.99999999999) {
91	++indx.phi;
92	frct.phi = 0.;
93	} else
94	return;
95	} else if ((in = indx.phi + 1) >= ct.lim.phi) {
96	if (in == ct.lim.phi && frct.phi < 1e-11) {
97	--indx.phi;
98	frct.phi = 1.;
99	} else
100	return;
101	}
102	boost::int32_t index = indx.phi * ct.lim.lam + indx.lam;
103	pj_ctable::flp_t const& f00 = ct.cvs[index++];
104	pj_ctable::flp_t const& f10 = ct.cvs[index];
105	index += ct.lim.lam;
106	pj_ctable::flp_t const& f11 = ct.cvs[index--];
107	pj_ctable::flp_t const& f01 = ct.cvs[index];
108	CalcT m00, m10, m01, m11;
109	m11 = m10 = frct.lam;
110	m00 = m01 = 1. - frct.lam;
111	m11 *= frct.phi;
112	m01 *= frct.phi;
113	frct.phi = 1. - frct.phi;
114	m00 *= frct.phi;
115	m10 *= frct.phi;
116	out_lon = m00 * f00.lam + m10 * f10.lam +
117	m01 * f01.lam + m11 * f11.lam;
118	out_lat = m00 * f00.phi + m10 * f10.phi +
119	m01 * f01.phi + m11 * f11.phi;
120	}
121
122	// Originally implemented in nad_cvt.c
123	template <bool Inverse, typename CalcT>
124	inline void nad_cvt(CalcT const& in_lon, CalcT const& in_lat,
125	CalcT & out_lon, CalcT & out_lat,
126	pj_gi const& gi)
127	{
128	static const int max_iterations = 10;
129	static const CalcT tol = 1e-12;
130	static const CalcT toltol = tol * tol;
131	static const CalcT pi = math::pi<CalcT>();
132
133	// horizontal grid expected
134	BOOST_GEOMETRY_ASSERT_MSG(gi.format != pj_gi::gtx,
135	"Vertical grid cannot be used in horizontal shift.");
136
137	pj_ctable const& ct = gi.ct;
138
139	// TODO: implement differently
140	if (in_lon == HUGE_VAL)
141	{
142	out_lon = HUGE_VAL;
143	out_lat = HUGE_VAL;
144	return;
145	}
146
147	// normalize input to ll origin
148	pj_ctable::lp_t tb;
149	tb.lam = in_lon - ct.ll.lam;
150	tb.phi = in_lat - ct.ll.phi;
151	tb.lam = adjlon (tb.lam - pi) + pi;
152
153	pj_ctable::lp_t t;
154	nad_intr(tb.lam, tb.phi, t.lam, t.phi, ct);
155	if (t.lam == HUGE_VAL)
156	{
157	out_lon = HUGE_VAL;
158	out_lat = HUGE_VAL;
159	return;
160	}
161
162	if (! Inverse)
163	{
164	out_lon = in_lon - t.lam;
165	out_lat = in_lat - t.phi;
166	return;
167	}
168
169	t.lam = tb.lam + t.lam;
170	t.phi = tb.phi - t.phi;
171
172	int i = max_iterations;
173	pj_ctable::lp_t del, dif;
174	do
175	{
176	nad_intr(t.lam, t.phi, del.lam, del.phi, ct);
177
178	// This case used to return failure, but I have
179	// changed it to return the first order approximation
180	// of the inverse shift. This avoids cases where the
181	// grid shift into this grid came from another grid.
182	// While we aren't returning optimally correct results
183	// I feel a close result in this case is better than
184	// no result. NFW
185	// To demonstrate use -112.5839956 49.4914451 against
186	// the NTv2 grid shift file from Canada.
187	if (del.lam == HUGE_VAL)
188	{
189	// Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.
190	break;
191	}
192
193	dif.lam = t.lam - del.lam - tb.lam;
194	dif.phi = t.phi + del.phi - tb.phi;
195	t.lam -= dif.lam;
196	t.phi -= dif.phi;
197
198	}
199	while (--i && (dif.lamdif.lam + dif.phidif.phi > toltol)); // prob. slightly faster than hypot()
200
201	if (i==0)
202	{
203	// Inverse grid shift iterator failed to converge.
204	out_lon = HUGE_VAL;
205	out_lat = HUGE_VAL;
206	return;
207	}
208
209	out_lon = adjlon (t.lam + ct.ll.lam);
210	out_lat = t.phi + ct.ll.phi;
211	}
212
213
214	/************************************************************************/
215	/* find_grid() */
216	/* */
217	/* Determine which grid is the correct given an input coordinate. */
218	/************************************************************************/
219
220	// Originally find_ctable()
221	// here divided into grid_disjoint(), find_grid() and load_grid()
222
223	template <typename T>
224	inline bool grid_disjoint(T const& lam, T const& phi,
225	pj_ctable const& ct)
226	{
227	double epsilon = (fabs(ct.del.phi)+fabs(ct.del.lam))/10000.0;
228	return ct.ll.phi - epsilon > phi
229	\|\| ct.ll.lam - epsilon > lam
230	\|\| (ct.ll.phi + (ct.lim.phi-1) * ct.del.phi + epsilon < phi)
231	\|\| (ct.ll.lam + (ct.lim.lam-1) * ct.del.lam + epsilon < lam);
232	}
233
234	template <typename T>
235	inline pj_gi * find_grid(T const& lam,
236	T const& phi,
237	std::vector<pj_gi>::iterator first,
238	std::vector<pj_gi>::iterator last)
239	{
240	pj_gi * gip = NULL;
241
242	for( ; first != last ; ++first )
243	{
244	// skip tables that don't match our point at all.
245	if (! grid_disjoint(lam, phi, first->ct))
246	{
247	// skip vertical grids
248	if (first->format != pj_gi::gtx)
249	{
250	gip = boost::addressof(*first);
251	break;
252	}
253	}
254	}
255
256	// If we didn't find a child then nothing more to do
257	if( gip == NULL )
258	return gip;
259
260	// Otherwise use the child, first checking it's children
261	pj_gi * child = find_grid(lam, phi, first->children.begin(), first->children.end());
262	if (child != NULL)
263	gip = child;
264
265	return gip;
266	}
267
268	template <typename T>
269	inline pj_gi * find_grid(T const& lam,
270	T const& phi,
271	pj_gridinfo & grids,
272	std::vector<std::size_t> const& gridindexes)
273	{
274	pj_gi * gip = NULL;
275
276	// keep trying till we find a table that works
277	for (std::size_t i = 0 ; i < gridindexes.size() ; ++i)
278	{
279	pj_gi & gi = grids[gridindexes[i]];
280
281	// skip tables that don't match our point at all.
282	if (! grid_disjoint(lam, phi, gi.ct))
283	{
284	// skip vertical grids
285	if (gi.format != pj_gi::gtx)
286	{
287	gip = boost::addressof(gi);
288	break;
289	}
290	}
291	}
292
293	if (gip == NULL)
294	return gip;
295
296	// If we have child nodes, check to see if any of them apply.
297	pj_gi * child = find_grid(lam, phi, gip->children.begin(), gip->children.end());
298	if (child != NULL)
299	gip = child;
300
301	// if we get this far we have found a suitable grid
302	return gip;
303	}
304
305
306	template <typename StreamPolicy>
307	inline bool load_grid(StreamPolicy const& stream_policy, pj_gi_load & gi)
308	{
309	// load the grid shift info if we don't have it.
310	if (gi.ct.cvs.empty())
311	{
312	typename StreamPolicy::stream_type is;
313	stream_policy.open(is, gi.gridname);
314
315	if (! pj_gridinfo_load(is, gi))
316	{
317	//pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
318	return false;
319	}
320	}
321
322	return true;
323	}
324
325
326	/************************************************************************/
327	/* pj_apply_gridshift_3() */
328	/* */
329	/* This is the real workhorse, given a gridlist. */
330	/************************************************************************/
331
332	// Generic stream policy and standard grids
333	template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename Grids>
334	inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
335	Range & range,
336	Grids & grids,
337	std::vector<std::size_t> const& gridindexes,
338	grids_tag)
339	{
340	typedef typename boost::range_size<Range>::type size_type;
341
342	// If the grids are empty the indexes are as well
343	if (gridindexes.empty())
344	{
345	//pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
346	//return PJD_ERR_FAILED_TO_LOAD_GRID;
347	return false;
348	}
349
350	size_type point_count = boost::size(range);
351
352	for (size_type i = 0 ; i < point_count ; ++i)
353	{
354	typename boost::range_reference<Range>::type
355	point = range::at(range, i);
356
357	CalcT in_lon = geometry::get_as_radian<0>(point);
358	CalcT in_lat = geometry::get_as_radian<1>(point);
359
360	pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes);
361
362	if ( gip != NULL )
363	{
364	// load the grid shift info if we don't have it.
365	if (! gip->ct.cvs.empty() \|\| load_grid(stream_policy, *gip))
366	{
367	// TODO: use set_invalid_point() or similar mechanism
368	CalcT out_lon = HUGE_VAL;
369	CalcT out_lat = HUGE_VAL;
370
371	nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
372
373	// TODO: check differently
374	if ( out_lon != HUGE_VAL )
375	{
376	geometry::set_from_radian<0>(point, out_lon);
377	geometry::set_from_radian<1>(point, out_lat);
378	}
379	}
380	}
381	}
382
383	return true;
384	}
385
386	// Generic stream policy and shared grids
387	template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename SharedGrids>
388	inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
389	Range & range,
390	SharedGrids & grids,
391	std::vector<std::size_t> const& gridindexes,
392	shared_grids_tag)
393	{
394	typedef typename boost::range_size<Range>::type size_type;
395
396	// If the grids are empty the indexes are as well
397	if (gridindexes.empty())
398	{
399	//pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
400	//return PJD_ERR_FAILED_TO_LOAD_GRID;
401	return false;
402	}
403
404	size_type point_count = boost::size(range);
405
406	// local storage
407	pj_gi_load local_gi;
408
409	for (size_type i = 0 ; i < point_count ; )
410	{
411	bool load_needed = false;
412
413	CalcT in_lon = 0;
414	CalcT in_lat = 0;
415
416	{
417	typename SharedGrids::read_locked lck_grids(grids);
418
419	for ( ; i < point_count ; ++i )
420	{
421	typename boost::range_reference<Range>::type
422	point = range::at(range, i);
423
424	in_lon = geometry::get_as_radian<0>(point);
425	in_lat = geometry::get_as_radian<1>(point);
426
427	pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
428
429	if (gip == NULL)
430	{
431	// do nothing
432	}
433	else if (! gip->ct.cvs.empty())
434	{
435	// TODO: use set_invalid_point() or similar mechanism
436	CalcT out_lon = HUGE_VAL;
437	CalcT out_lat = HUGE_VAL;
438
439	nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
440
441	// TODO: check differently
442	if (out_lon != HUGE_VAL)
443	{
444	geometry::set_from_radian<0>(point, out_lon);
445	geometry::set_from_radian<1>(point, out_lat);
446	}
447	}
448	else
449	{
450	// loading is needed
451	local_gi = *gip;
452	load_needed = true;
453	break;
454	}
455	}
456	}
457
458	if (load_needed)
459	{
460	if (load_grid(stream_policy, local_gi))
461	{
462	typename SharedGrids::write_locked lck_grids(grids);
463
464	// check again in case other thread already loaded the grid.
465	pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
466
467	if (gip != NULL && gip->ct.cvs.empty())
468	{
469	// swap loaded local storage with empty grid
470	local_gi.swap(*gip);
471	}
472	}
473	else
474	{
475	++i;
476	}
477	}
478	}
479
480	return true;
481	}
482
483
484	/************************************************************************/
485	/* pj_apply_gridshift_2() */
486	/* */
487	/* This implementation uses the gridlist from a coordinate */
488	/* system definition. If the gridlist has not yet been */
489	/* populated in the coordinate system definition we set it up */
490	/* now. */
491	/************************************************************************/
492
493	template <bool Inverse, typename Par, typename Range, typename ProjGrids>
494	inline bool pj_apply_gridshift_2(Par const& /defn/, Range & range, ProjGrids const& proj_grids)
495	{
496	/*if( defn->catalog_name != NULL )
497	return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
498	x, y, z );*/
499
500	/*std::vector<std::size_t> gridindexes;
501	pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"),
502	grids.storage_ptr->stream_policy,
503	grids.storage_ptr->grids,
504	gridindexes);*/
505
506	// At this point the grids should be initialized
507	if (proj_grids.hindexes.empty())
508	return false;
509
510	return pj_apply_gridshift_3
511	<
512	Inverse, typename Par::type
513	>(proj_grids.grids_storage().stream_policy,
514	range,
515	proj_grids.grids_storage().hgrids,
516	proj_grids.hindexes,
517	typename ProjGrids::grids_storage_type::grids_type::tag());
518	}
519
520	template <bool Inverse, typename Par, typename Range>
521	inline bool pj_apply_gridshift_2(Par const& , Range & , srs::detail::empty_projection_grids const& )
522	{
523	return false;
524	}
525
526
527	} // namespace detail
528
529	}}} // namespace boost::geometry::projections
530
531	#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP