@C_R  Rouben forgot to include:

addcoords(z_cylindrical,[z, r,theta],[r*cos(theta),r*sin(theta),z]);

For a unique solution, you will need to define a well-order (see wiki). 

The attachment contains an image hint. (For the moment mapleprimes refuses to inline it).

pic1.zip

@acer  For testing (without modifying the polygon) it is indeed more convenient to use clicks. In this case, it is enough to replace 

Explore(T(x,y), x=0..80, y=0..100);

with

Explore(T(x,y), x=0..80, y=0..100, markers=[[x,y]]);

@Carl Love I always enjoy your optimized but too cryptic (for my taste) procedures!

You cannot omit the early return when P is in L, otherwise the result will be wrong in most cases.

@acer  Nice. Unfortunately the blue point cannot be dragged over the boundary, so, the undefined value cannot be tested.

@mmcdara Here is a test using Explore. For older versions of Maple, use the previous Wind.

restart;
Wind:=proc(L::listlist, P::list)
  local p,u,k, x0:=P[1],y0:=P[2], t;
  if member(P,L) then return undefined fi;
  u:=seq(arctan(p[2]-y0,p[1]-x0),p=L); if L[1]<>L[-1] then u:=u,u[1] fi;
  add(piecewise(abs((t:=u[k+1]-u[k]))=Pi,undefined, t<-Pi,+1, t>Pi,-1, 0), k=1..nops([u])-1)
end:
K:=[[10,10], [10,90], [70,90], [70,50], [30,50], [30,30], [60,30], [60,80], [20,80], [20,70], [40,70], [40,40], [20,40], [20,60], [50,60], [50,10], [10,10]]:
T:=proc(x,y)  local w:=Wind(K, [x,y]);
  uses plots;
  display(plottools:-polygon(K),color=yellow,linestyle=solid,thickness=2, 
    pointplot(K, color=red, symbolsize=16, symbol=solidcircle),
    pointplot([x,y], color=blue, symbolsize=16, symbol=solidbox),
    view=[0..80, 0..100], title=typeset([x,y],'wind'=w));
end:
Explore(T(x,y), x=0..80, y=0..100);

@Rouben Rostamian  

alpha_max / Pi  is one of the real roots of the polynomial

9*t^6 - 54*t^5 + 48*t^4 + 168*t^3 - 336*t^2 + 192*t - 32

@nm A much larger maxslope:=100000:  is probably acceptable in most non-pathological cases.

@mmcdara The free and generating families are useful notions. Usually it is simpler to check these two, instead of a single condition. E.g., if a family F is free in a vector space E, then F is a basis for the vector subspace generated by F in E.

@Christopher2222  The strange fact is that evalb and select do not confirm numboccur.
numboccur(L, 1.) = nops(select(`=`, L, 1.));  # 2 = 3

@Carl Love Yes, of course.

Q:=eval([x,y], solve([x+2*y-1, 2*(x-6) - (y+3)])): # projection onto directrix
F:=2*[6,-3]-Q:
P:=(x-F[1])^2 + (y-F[2])^2 = (x+2*y-1)^2 / 5:
with(plots): display(implicitplot([P, x+2*y-1], x=-10..10,y=-10..10), pointplot([[6,-3], F], color=[red,blue]));

@JAMET You have already an answer. What did you not understand?

@mmcdara The unicity (uniqueness) of the linear combination makes the family free. So, the definitions are the same.