Enhancement: index polygon on the sphere

[kilir]

Hi Ralf, so far, I did not find an easy way to accomplish the following: If I have a polygon on a sphere defined by some vectors, is there some way to tell whether a point is "inside" or "outside"?

I guess "inside" or "outside" is arbitrary, but it should be possible to tell to which area a given vector belongs. Is there a simple way to do that? So far I try to approximate the shape with many small cones, but this is not very elegant nor I think really the solution, although it totally works for what I do and, it's not urgent at all.

Cheers, Rüdiger

[ralfHielscher]

Hi Rüdiger,

so etwas könnte ich tatsächlich auch gut gebrauchen. Ich denke,

https://de.wikipedia.org/wiki/Punkt-in-Polygon-Test_nach_Jordan

ist ein guter Ansatz, dass sich Schnittpunkte auch auf der Sphäre schnell berechnen lassen.

Ralf.

---

Ralf Hielscher Tel: +371-531-38556 Fakultät für Mathematik +371-531-22200 (Sekr.) Technische Universität Chemnitz Fax: +371-531-22109 Reichenhainer Str. 39 E-mail: [email protected] D-09126 Chemnitz http://www.tu-chemnitz.de/~rahi

---

On Wed, Jan 24, 2018 at 12:34 PM, kilir [email protected] wrote:

Hi Ralf, so far, I did not find an easy way to accomplish the following: If I have a polygon on a sphere defined by some vectors, is there some way to tell whether a point is "inside" or "outside"?

I guess "inside" or "outside" is arbitrary, but it should be possible to tell to which area a given vector belongs. Is there a simple way to do that? So far I try to approximate the shape with many small cones, but this not very elegant no I think really the solution, although it totally works for what I do and, it's not urgent at all.

[image: polygon] https://user-images.githubusercontent.com/11697557/35329956-f50ec43c-0101-11e8-9909-463c649d3314.png

Cheers, Rüdiger

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/mtex-toolbox/mtex/issues/299, or mute the thread https://github.com/notifications/unsubscribe-auth/AIOE6QEPKgJ_l2fM_48GdB2EP6TsWJz-ks5tNxU0gaJpZM4RrHfC .

[kilir]

Hi Ralf, I'm not so sure that this would work on the sphere as a great circle would always return to it's origin, and so it would also have also even number of intersections, wouldn't it? Rüdiger

Edit: Apologies for my scepticism, apparently this indeed works somehow. "Locating a point on a spherical surface relative to a spherical polygon of arbitrary shape" Bevis & Chatelain; https://link.springer.com/article/10.1007/BF00894449

[ralfHielscher]

Hi Rüdiger,

what about this: take an arbitrary point which we call "inside" and the connect the point in question with this point by an bigcircle. If this connecting line intersect an odd number of vertices the point is outside otherwise inside.

Ralf.

[kilir]

Hi, Ralf, I guess my problem was that I thought the line should go once around but testing just for the connection between the two points makes a lot sense. I'll try this. Cheers, Rüdiger

[kilir]

Hi Ralf, I tried to do this and it seems to work. (It probably could get rid of one of the loops and maybe there's a more reasonable way to find a point 'inside', because the mean(v) will certainly fail in some situations.

function ind = inpolygonSP(v,c,p)
% silly check for a point to be inside a polygon on a sphere
% input:
% v - point to be tested @vector3d
% c - a point to be known inside the poylgon @vector3d
%     for easy shapes mean(v) might do
% p - a list of points @vector3d defining the polygon on the sphere
%
% output:
% ind - 1 if it is considered to be inside
%

% TODO: get rid of the loop


% don't even start
if length(p) < 3
    error('polygon needs at least 3 points')
end    

%close the polygon
if p(1) ~= p(end)
    p(end+1) = p(1);
end

ind = ones(length(v),1);
for i=1:length(v)
    isIntersection = zeros(length(p)-1,1);
    % iterate through polygon
    for j=1:length(p)-1
        % find a potential intersection
        p1 = cross(p(j),p(j+1)); % pole to the plane deined by two vertex vectors
        p2 = cross(c,v(i));  % pole to the plane defined by inside point and testpoint
        pP = cross(p2,p1); % intesection is the pole to plane defined p1 and p2
        
        % some attempt to figure out if the intersection is actually on the segment
        % assuming that the angle between two endpoints of a segment should
        % be larger that the angle between the potential intersection and
        % any end point - needs to be tested for the polygon vertices and for the
        % center and test point pair
        
        if angle(p(j),p(j+1)) > max(angle(pP,p([j j+1]))) && ...
           angle(c,v(i))     > max(angle(pP,[c,v(i)]))
           isIntersection(j)=1; % the interesection is on the edge of the polygon
        else
           isIntersection(j)=0; % no, is is not
        end
        
        % sum of intersections:
        % ifthere are more than zero intersections and the number is odd,
        % point should be outside
        if sum(isIntersection) > 0 && mod(sum(isIntersection),2)==1  
            ind(i) = 0;  % is outside
        end
        
    end
end

Cheers, Rüdiger

[ralfHielscher]

Hi Rüdiger,

it might help that

angle(v1,v2,N)

returns a signed angle relative to the normal N.

Ralf.

[kilir]

After a little bit of figuring out what to do with indenting segments, it seems this is also quite faster:

function isInside = inpolygonSP(v,c,p)
% silly check for a point to be inside a polygon on a sphere
% input:
% v - point to be tested @vector3d
% c - a point to be known inside the poylgon @vector3d
%     for easy shapes mean(v) might do
% p - a list of points @vector3d defining the polygon on the sphere
%
% output:
% isInside - 1 if it is considered to be inside
%

% don't even start
if length(p) < 3
    error('polygon needs at least 3 points')
end

v.antipodal=0;
p = reshape(p,[],1);
v = reshape(v,[],1);

% find a potential intersection
p1 = cross(p(1:end),p([2:end,1])); % pole to the plane deined by two vertex vectors
p2 = cross(c,v);  % pole to the plane defined by inside point and testpoint
pP = cross_outer(p1,p2); % intesections is the poles to plane defined by p1 and p2

% if pP is between p(1) and p(2), the angle around
% p(1) X p(2) should have opposite signes

% only 1 of cp1 or cp2 should be >=pi which means  means they have opposite
% signs

cp1 =  angle(p(1:end),pP,p1) >=pi;
cp2 =  angle(p([2:end, 1]),pP,p1) >=pi;

isOnPolygon = cp1 ~= cp2;

% some special condition to work around "indenting parts of the
% polygon where direction of pP reverses
pP(cp1 == 1)=pP(cp1 == 1).*-1;

% for the testline we need opposite conditions
isOnTestline = angle(c,pP,p2') <= pi & ...
               angle(v',pP,p2') >= pi; 


% intersections should be true for both
% test for inside: no or even number of intersections
isInside =  mod(sum(isOnPolygon & isOnTestline),2) == 0;
end

[kilir]

To select a polygon by hand, there might be an easier way, but that's what I came up with so far:

function pos = clickRegionFromShpere(mtexFig)
% select points on a pole figure
% works only on tight grids or smooth contours
% waits for user click and stops when key s" is pressed

disp('Click on a polefigure to selects points, press "s" key to stop')

set(gcf,'CurrentCharacter',' ' );
disp('Select region, when done, hit "s" ')
pos=[];
datacursormode on
while ~strcmpi(get(gcf,'CurrentCharacter'),'s')
    t=waitforbuttonpress;
    if t==0
        % wait for a click and do something
        [posi] = getDataCursorPos(mtexFig);
        hold on
        plot(posi,'parent',mtexFig.gca)
        hold off
        pos=[pos posi]
    end
    
    if strcmpi(get(gcf,'CurrentCharacter'),'s');
        break
    end
end

end

The only problem so far is that it needs a rather tight grid or a contour plot, but it cannot select points in an otherwise empty plot. Cheers,Rüdiger

[ralfHielscher]

Hi Rüdiger,

what is the general purpose with such region? If it is of bigger interest - it would probably make sense to implement it as a class similar to sphericalRegion. Maybe those could be derived from the same base class and thus share some properties.

Ralf.

[kilir]

Hi Ralf, the purpose was to select an area on a polefigure which cannot be simply described by a cone and and then use this condition to e.g. explore grains with a certain axis within this area or colorcode them. Here's a small example with a cross girdle and one might be wondering if there is any measurable difference between grains in the 'green' or in the ' red' leg. The quickest way to separate the areas of interest in the polefigure was just to click along them.

Cheers, Rüdiger

[kilir]

Hi Ralf, considering the new way to define a spherical region using .byVertices() I thought the above (polygon on a sphere and the check wether somethign is inside or not) would be redundand since the spherical region also has checkInside. So far - apart from some not yet reproducible problems - it seems to work.

Cheers, Rüdiger

[ralfHielscher]

Hi Rüdiger,

indeed when representing spherical polyhedrons by spherical regions. The spherical polyhedron need to be convex. The reason lies in the way those spherical regions are represented.

Ralf.