Graham scan convex hull algorithm - include all points on boundary I have am implementing the Graham scan algorithm to find the convex hull of a set of (two-dimensional) points. (My implementation is in Haskell in case anyone wants to know.) The problem is that not all of the boundary points are included. In fact, my implementation is highly sensitive to the order of the points in the input (since the input is technically a list, not a set). I am looking for some help from a mathematical perspective in order to fix the problem.
To start, let me describe how I understand the algorithm:

Input: $S \subseteq \mathbb{R}^2$.
Output: The convex hull $H$ of $S$

*

*Define the point $P$ as follows: First find the subset $Y$ of $S$ of all points $(x, y)$ such that $y$ is minimum. Then $P$ is the point in $Y$ such that $x$ is minimum.


*Define $\preceq$ on $S$ as $P_1 \preceq P_2$ iff $\theta_1 \leq \theta_2$ where $(r_i, \theta_i)$ is the polar representation of the vector $\vec{PP_i}$.


*Define the sequence $A = \{P_i\}$ such that $P_i \in S$ and $P_i \preceq P_{i+1}$ for each $i$.


*Define $H_0 = A$. Define the subsequence $H_i=\{Q_i\}$ of $H_{i-1}$ by removing the point $P_j$ from $H_{i-1}$ such that the angle $\angle P_jP_{j+1}P_{j+2}$ makes a right-hand turn and $j$ is minimum.


*Let $H=H_i$ where $H_i$ contains no angle $\angle P_jP_{j+1}P_{j+2}$ which makes a right-hand turn.

Now for a more concrete example to explain the problem:
Let $S = \{(x, y) : x, y \in \mathbb{Z}, 0 \leq x, y \leq 4\}$. (Since my implementation uses a list rather than a set, $S$ is actually a the sequence $S=\{P_i\}$ of these points in lexicographical order.) From step 1, $P = (0, 0)$. My implementation also uses a stable sort, so for $P_i, P_j \in S$, if $\theta_i = \theta_j$ (again using polar coordinates for the points), then $P_i \prec P_j$ if $i < j$. This creates a problem as I consider the points $(0, 1)$, $(0, 2)$, $(0, 3)$, and $(0, 4)$ in step 4 because they appear in the sequence $S$ in that order and they are the vertices for angles with right-handed turns when they are considered. For example, the angle formed by the points $(x, y)$, $(0, 1)$, $(0, 2)$ where $x \not=0$ is a right-handed turn.
How do I modify the ordering of step 2 so that the points on the boundary appear in the correct order to obtain straight angles as I consider them in step 4?
p.s. I have decided to state my question in mathematical terms, but if you think I would get more/better feedback from programmers, I will post it on SO instead and try to frame it in programming terms.
Addendum:
My implementation in Haskell keeps the points
$\{(0, 0), (1, 0), (2, 0), (3, 0), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4), (3, 4), (2, 4), (1, 4)\}$
as it searches for the convex hull. The remaining points to check are $\{(0, 1), (0, 2), (0, 3), (0, 4)\}$, in that order. The angle formed by $(1, 4)$, $(0, 1)$, $(0, 2)$ is a right-hand turn, so$(0, 1)$ is discarded. Similarly, the points $(0, 2)$ and $(0, 3)$ are discarded. Finally, $(0, 4)$ is correctly retained. I need a way to ensure that these points are encountered in the order $(0, 4)$, $(0, 3)$, $(0, 2)$, $(0, 1)$ so that no right-hand turns will be found and all of these boundary points will be retained. This is just an example, so I need to find a general solution to avoid this kind of problem.
 A: If you just want all the degenerate points on the boundary of the convex hull to be included, you can find the convex hull, then test each point individually and see if it lies on an edge of the convex hull, then split edges, and insert points appropriately.
A: The solution that I found after the discussion in the comments to my question is to change the ordering defined in step 2:

Let $(r_i, \theta_i)$ be the polar representation of the vector $\vec{PP_i}$ for each $P_i \in S$. Furthermore, let $\theta_{min}=\min\{\theta_i\}$ and $\theta_{max}=\max\{\theta_i\}$. Define $\prec$ on $S$ as $P_i \prec P_j$ iff $\theta_i < \theta_j$ or ($\theta_i = \theta_j = \theta_{min}$ and $r_i < r_j$) or ($\theta_i = \theta_j = \theta_{max}$ and $r_j < r_i$).

This treats the minimum and maximum values of $\theta_i$ as special cases and (as far as I can tell) sorts the points on the boundary of the convex hull in "counter-clockwise order" so that the rest of the algorithm will now include the boundary points in the output sequence.
