PreviousUpNext

15.4.1030  src/lib/std/2d/geometry2d.pkg

## geometry2d.pkg
#
# Basic geometric types and operations.
#
# The 'X' on the name is for X Windows; this
# file was originally part of x-kit, which is
# still its major user.  However, I now regard
# it as platform-independent code.
#                        -- 2014-06-27 CrT

# Compiled by:
#     src/lib/std/standard.lib


###                      "Much learning does not teach understanding."
###
###                                               -- Heraclitus  (540BC-480BC, On the Universe)


###                      "Logic will get you from A to B. Imagination will take you everywhere."
###
###                                               -- Albert Einstein


                                                                                # Geometry2d            is from   src/lib/std/2d/geometry2d.api
stipulate
    package rc  =  range_check;                                                 # range_check           is from   src/lib/std/2d/range-check.pkg
    package ebf =  eight_byte_float;                                            # eight_byte_float      is from   src/lib/std/eight-byte-float.pkg
    package lms =  list_mergesort;                                              # list_mergesort        is from   src/lib/src/list-mergesort.pkg
    #
    nb =  log::note_on_stderr;                                                  # log                   is from   src/lib/std/src/log.pkg
herein

    package geometry2d {
        #
        stipulate

            fun min (col: Int, row) =   col < row  ??  col  ::  row;
            fun max (col: Int, row) =   col > row  ??  col  ::  row;

        herein

            # Geometric types (from xlib.h)
            #
            Point =   { col:  Int,
                        row:  Int
                      };

            Line =  (Point, Point);

            Size =  { wide:  Int,
                      high:  Int
                    };

            Box  =  { col:  Int,
                      row:  Int,
                      #
                      wide:  Int,
                      high:  Int
                    };

            Arc =   { row:  Int,
                      col:  Int,
                      #
                      wide:  Int,
                      high:  Int,
                      #
                      start_angle:  Float,              # In degrees, with zero angle at 3 o'clock, increasing counterclockwise.  Use positive angles from 0.0 to 360.0.
                      fill_angle:   Float               # Draw a pie-slice of this many degrees starting at start_angle and running counterclockwise from there.
                    };
                                                        # Examples:
                                                        #     Upper-right quadrant ==  { ..., start_angle =>   0.0,  fill_angle =>  90.0 }
                                                        #     Upper-left  quadrant ==  { ..., start_angle =>  90.0,  fill_angle =>  90.0 }
                                                        #     Lower-left  quadrant ==  { ..., start_angle => 180.0,  fill_angle =>  90.0 }
                                                        #     Lower-right quadrant ==  { ..., start_angle => 270.0,  fill_angle =>  90.0 }
                                                        #     Upper half           ==  { ...  start_angle =>   0.0,  fill_angle => 180.0 };
                                                        #     Lower half           ==  { ...  start_angle => 180.0,  fill_angle => 180.0 };
                                                        #     Full disk            ==  { ..., start_angle =>   0.0,  fill_angle => 360.0 }

            Arc64 = { col:     Int,
                      row:     Int,
                      # 
                      wide:    Int,                     # If wide != high the arc drawn is from an ellipse rather than a circle.
                      high:    Int,
                      #
                      angle1:  Int,                     # In degrees * 64, with zero angle at 3 o'clock, increasing counterclockwise.
                      angle2:  Int                      # Arc is drawn from angle1, extending for angle2/64 degrees.
                    };

            # The size and position of a window         # XXX BUGGO FIXME This belongs in xclient, not stdlib.
            # with respect to its parent:
            #
            Window_Site =   { upperleft:        Point,
                              size:             Size,
                              border_thickness: Int
                            };

            package point {

                # Points:
                #
                zero = { col => 0, row => 0 };

                fun col ({ col, ... }: Point) =  col;
                fun row ({ row, ... }: Point) =  row;

                fun add      ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 } ) =  { col=>(col1+col2), row=>(row1+row2) };
                fun subtract ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 } ) =  { col=>(col1-col2), row=>(row1-row2) };

                fun scale ({ col, row }, s ) =  { col=>s*col, row=>s*row };

                fun ne ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 != col2) or  (row1 != row2);
                fun eq ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 == col2) and (row1 == row2);
                fun lt ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 <  col2) and (row1 <  row2);
                fun le ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 <= col2) and (row1 <= row2);
                fun gt ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 >  col2) and (row1 >  row2);
                fun ge ({ col=>col1, row=>row1 }, { col=>col2, row=>row2 }) =  (col1 >= col2) and (row1 >= row2);

                fun add_size ({ col, row }, { wide, high } )
                    =
                    { col =>  col + wide,
                      row =>  row + high
                    };
                #
                fun clip ({ col, row }, { wide, high })
                    =
                    { col =>  if (col <= 0)  0; elif (col < wide)  col; else (wide - 1); fi,
                      row =>  if (row <= 0)  0; elif (row < high)  row; else (high - 1); fi
                    };

                fun in_box ({ col=>px, row=>py }, { col, row, wide, high } )
                    =
                    px >=  col    and
                    py >=  row    and
                    px < col+wide and
                    py < row+high;


                fun compare_xy (p1: Point, p2: Point)                                   # Comparison fn to sort points by X (and by Y when X coords match).
                    =                                                                   # Used in convex_hull (below); generally useful to induce a total ordering on points.
                    if       (p1.col == p2.col)
                        if   (p1.row == p2.row) EQUAL;
                        elif (p1.row >  p2.row) GREATER;
                        else                    LESS;           fi;
                    elif     (p1.col >  p2.col) GREATER;
                    else                        LESS;
                    fi;
                                                                                        # We'll probably want a compare_yx one of these days to sort points along the Y axis.
                fun mean (points: List(Point))
                    =
                    { row =>  int::mean  (map .row points),
                      col =>  int::mean  (map .col points)
                    };
            };

            package size {
                #
                fun add      ({ wide=>w1, high=>h1 }, { wide=>w2, high=>h2 } ) =  { wide=>(w1+w2), high=>(h1+h2) };
                fun subtract ({ wide=>w1, high=>h1 }, { wide=>w2, high=>h2 } ) =  { wide=>(w1-w2), high=>(h1-h2) };
                fun scale    ({ wide,     high     }, s                      ) =  { wide=>wide*s,  high=>high*s  };
                #
                fun eq       ({ wide=>w1, high=>h1 }, { wide=>w2, high=>h2 } ) =  (w1==w2 and h1==h2);
            };

            package box {
                #
                zero = { col  => 0,
                         row  => 0,
                         high => 0,
                         wide => 0
                       };

                fun clone_box_at
                      (
                        { high, wide, ... }:    Box,
                        { row, col }:           Point
                      )
                    =
                    { row, col, high, wide };

                fun ne ({ col=>col1, row=>row1, high=>high1, wide=>wide1 }, { col=>col2, row=>row2, high=>high2, wide=>wide2 }) =  (col1 != col2) or  (row1 != row2) or  (high1 != high2) or  (wide1 != wide2);
                fun eq ({ col=>col1, row=>row1, high=>high1, wide=>wide1 }, { col=>col2, row=>row2, high=>high2, wide=>wide2 }) =  (col1 == col2) and (row1 == row2) and (high1 == high2) and (wide1 == wide2);

                fun make ({ col, row }, { wide, high })
                    =
                    { col, row, wide, high };


                fun upperleft_and_size ({ col, row, wide, high } )
                    =
                    ({ col, row }, { wide, high } );


                fun upperleft ({ col, row, ... }: Box)
                    =
                    { col, row };

                fun lowerright ({ col, row, high, wide }: Box)          # Returns  { col => box.col + box.wide - 1,  row => box.row + box.high - 1 }
                    =
                    { col => col + wide - 1,
                      row => row + high - 1
                    };

                fun lowerright1 r                                       # Returns  { col => box.col + box.wide    ,  row => box.row + box.high     }
                    =
                    point::add_size (upperleft_and_size r);


                fun size ({ wide, high, ... }: Box)
                    =
                    { wide, high };


                fun area ({ wide, high, ... }: Box)
                    =
                    wide * high;

                fun to_points ({ col, row, wide, high }: Box)
                    =
                    [ { row => row       , col => col        },
                      { row => row + high, col => col        },
                      { row => row + high, col => col + wide },
                      { row => row       , col => col + wide }
                    ];

                fun box_corners ({ col, row, wide, high }: Box)
                    =
                    { upper_left  =>  { row => row       , col => col        },
                      lower_left  =>  { row => row + high, col => col        },
                      lower_right =>  { row => row + high, col => col + wide },
                      upper_right =>  { row => row       , col => col + wide }
                    };

                fun clip_point ({ col=>min_col, row=>min_row, wide, high }, { col, row } )
                    =
                    {
                      col => if (col <= min_col)  min_col; elif (col < min_col+wide)  col; else (min_col+wide - 1); fi,
                      row => if (row <= min_row)  min_row; elif (row < min_row+high)  row; else (min_row+high - 1); fi
                    };

                fun  translate ({ col, row, wide, high }, { col=>px, row=>py } ) =  { col=>col+px, row=>row+py, wide, high };
                fun rtranslate ({ col, row, wide, high }, { col=>px, row=>py } ) =  { col=>col-px, row=>row-py, wide, high };

                fun midpoint ({ col, row, wide, high } )
                    =
                    { col => col + (wide / 2),
                      row => row + (high / 2)
                    };

                fun intersect
                        ( { col=>col1, row=>row1, wide=>w1, high=>h1 },
                          { col=>col2, row=>row2, wide=>w2, high=>h2 }
                        )
                    =
                    (   (col1 < (col2+w2)) and (row1 < (row2+h2))
                    and (col2 < (col1+w1)) and (row2 < (row1+h1)));

                fun intersection
                        ( { col=>col1, row=>row1, wide=>w1, high=>h1 },
                          { col=>col2, row=>row2, wide=>w2, high=>h2 } )
                    =
                    {
                        col = max (col1, col2);
                        row = max (row1, row2);

                        cx = min (col1+w1, col2+w2);
                        cy = min (row1+h1, row2+h2);

                        if (col < cx  and  row < cy)
                            #
                            THE { col, row, wide=>(cx-col), high=>(cy-row) };
                        else
                            NULL;
                        fi;
                      };

                fun union (
                      r1 as { col=>col1, row=>row1, wide=>w1, high=>h1 },
                      r2 as { col=>col2, row=>row2, wide=>w2, high=>h2 }
                    )
                    =
                    if   (w1 == 0  or  h1 == 0)    r2;
                    elif (w2 == 0  or  h2 == 0)    r1;
                    else

                        col = min (col1, col2);
                        row = min (row1, row2);

                        cx = max (col1+w1, col2+w2);
                        cy = max (row1+h1, row2+h2);

                        { col, row, wide=>(cx-col), high=>(cy-row) };
                    fi;


                fun point_in_box
                    ( { col, row },
                      { col=> box_col, row=> box_row, wide, high }
                    )
                    =
                    col >= box_col          and                         # The >= < pattern here is intended to ensure that
                    row >= box_row          and                         # if we have a grid of boxes sharing vertices that
                    col <  box_col + wide   and                         # any given point is in exactly one box in the grid.
                    row <  box_row + high;

                fun point_on_box_perimeter
                    ( p as { col, row },
                      b as { col=> box_col, row=> box_row, wide, high }
                    )
                    =
                    point_in_box (p, b)
                    and
                    (
                        col == box_col                  or
                        col == box_col + wide - 1       or
                        row == box_row                  or
                        col == box_row + high - 1
                    );

                fun box_a_in_box_b
                    { a => { col=>col1, row=>row1, wide=>w1, high=>h1 },
                      b => { col=>col2, row=>row2, wide=>w2, high=>h2 }
                    }
                    =
                    col1 >= col2        and
                    row1 >= row2        and
                    col1+w1 <= col2+w2  and
                    row1+h1 <= row2+h2;

                fun make_nested_box
                      ( box as { col, row, wide, high }: Box,
                        by:   Int
                      )
                    =
                    # Create a box nested within given box,
                    # shrunk by given number of pixels:
                    #
                    if   (by   <= 0)   box;
                    elif (high <= 2)   box;
                    elif (wide <= 2)   box;
                    else
                        wide2 = wide / 2;
                        high2 = high / 2;

                        by = if (by > wide2)  wide2;
                             else             by;
                             fi;

                        by = if (by > high2)  high2;
                             else             by;
                             fi;

                        { row => row + by,  high => high - 2*by,
                          col => col + by,  wide => wide - 2*by
                        };
                    fi;



                # The symmetric difference of two sets is essentially
                # a geometric XOR operation;  it contains all elements
                # in either set but not both sets -- in other words,
                # the union minus the intersection:
                #
                #     http://en.wikipedia.org/wiki/Symmetric_difference
                #
                # Here we compute the symmetric difference of two
                # rectangles and return it as a list of rectangles:
                # 
                fun xor (r, r')
                    =
                    difference (r', r, difference (r, r',[]))
                    where
                        fun difference (r as { col=>x, row=>y, wide, high }, r', result_list)
                            =
                            case (intersection (r, r'))
                                #
                                NULL => r ! result_list;
                                #
                                THE { col=>ix, row=>iy, wide=>iwide, high=>ihigh }              # "i-" for "intersection-".
                                    =>
                                    {

                                        icx = ix + iwide;               # Opposite corner of
                                        icy = iy + ihigh;               # intersection box.

                                        # (x,y) is one corner of a rectangle,
                                        # (cx,cy) is the opposite corner.
                                        #       
                                        # Cyclically identify all parts of the rectangle
                                        # which project outside the borders of the above
                                        # intersection rectangle, adding each of them to
                                        # the result list and then shrinking the argument
                                        # rectangle correspondingly:
                                        #
                                        fun pare (x, y, cx, cy, result_list)
                                            =
                                            if   (  x < ix)  pare (ix,  y,  cx,  cy, ({ col=>x,   row=>y,   high=>cy-y,   wide=>ix-x   } ) ! result_list);      # Pare off the part to the left.
                                            elif (  y < iy)  pare ( x, iy,  cx,  cy, ({ col=>x,   row=>y,   high=>iy-y,   wide=>cx-x   } ) ! result_list);      # Pare off the part above. (Assuming y==0 is at top.)
                                            elif (icx < cx)  pare ( x,  y, icx,  cy, ({ col=>icx, row=>y,   high=>cy-y,   wide=>cx-icx } ) ! result_list);      # Pare off the part to the right.
                                            elif (icy < cy)  pare ( x,  y,  cx, icy, ({ col=>x,   row=>icy, high=>cy-icy, wide=>cx-x   } ) ! result_list);      # Pare off the part below.
                                            else
                                                result_list;
                                            fi;

                                        pare (x, y, x+wide, y+high, result_list);
                                    };
                            esac;
                    end;

                                                                                                                                                        # Compute intersection of lists of boxes.
                                                                                                                                                        #
                                                                                                                                                        # This is intended mainly for processing X EXPOSE events
                                                                                                                                                        # which contain lists of boxes.
                                                                                                                                                        #
                                                                                                                                                        # I'm expecting perhaps half a dozen boxes here, so the
                                                                                                                                                        # simple naive O(N**2) algorithm should be sufficient.
                                                                                                                                                        #
                                                                                                                                                        # We make no attempt to (say) merge geometrically adjacent boxes.
                                                                                                                                                        #
                fun intersect_box_with_boxes                                                                                                            # 
                      (                                                                                                                                 # 
                        box:            Box,
                        boxes:          List(Box)
                      ) 
                      :                 List(Box) 
                    =
                    do_boxes (boxes, [])                                                                                                                # 
                    where       
                        fun do_boxes ([],  result: List(Box))
                                =>
                                reverse result;                                                                                                         # Restore original order just in case caller cares.

                            do_boxes (box' ! rest,  result: List(Box))
                                =>
                                case (intersection (box,box'))
                                    #
                                    NULL  => do_boxes (rest,     result);                                                                               # These two boxes do not intersect so they add nothing to our result.
                                    THE i => do_boxes (rest, i ! result);                                                                               # Add the intersection of these two boxes to our result list and continue.
                                esac;
                        end;
                    end;

                fun intersect_boxes_with_boxes
                      (
                        boxes':         List(Box),
                        boxes:          List(Box)
                      ) 
                      :                 List(Box) 
                    =
                    list::cat (map do_box boxes')
                    where
                        fun do_box box
                            =
                            intersect_box_with_boxes (box, boxes);
                    end;


                fun vertical_lineseg_intersects_box (b as { row, col, high, wide }: Box, lower: Point, upper: Point)                                    # 
                    =                                                                                                                                   # Remember that origin is at upper-left.
                    lower.col >= col                    and
                    lower.col <  col + wide;

                fun bisect_box_vertically (b as { row, col, high, wide }: Box, lower: Point, upper: Point)
                    =
                    if (vertical_lineseg_intersects_box (b, lower, upper))
                        #
                        col' = lower.col;
                        #
                        [ { row, col,         high, wide =>        (col' - col) },
                          { row, col => col', high, wide => wide - (col' - col) }
                        ];
                    else
                        [ b ];
                    fi; 

                fun bisect_boxes_vertically (boxes: List(Box), lower: Point, upper: Point)
                    =
                    list::cat (map {. bisect_box_vertically (#b, lower, upper); }  boxes);



                fun horizontal_lineseg_intersects_box (b as { row, col, high, wide }: Box, left: Point, right: Point)                                   # 
                    =
                    left.row  >= row                    and
                    left.row  <  row + high;

                fun bisect_box_horizontally (b as { row, col, high, wide }: Box, left: Point, right: Point)
                    =
                    if (horizontal_lineseg_intersects_box (b, left, right))
                        #
                        row' = left.row;
                        #
                        [ { row,         col, wide, high =>        (row' - row) },
                          { row => row', col, wide, high => high - (row' - row) }
                        ];
                    else
                        [ b ];
                    fi; 

                fun bisect_boxes_horizontally (boxes: List(Box), left: Point, right: Point)
                    =
                    list::cat (map {. bisect_box_horizontally (#b, left, right); }  boxes);


                fun quadsect_box (b as { row, col, high, wide }: Box, p: Point)                                                                         # Divide box b into four nonintersecting subboxes with p as origin of lower-right one, if p is in interior of b.
                    =                                                                                                                                   # 
                    if   (high < 1 or wide < 1)                 [ b ];                                                                                  # Ignore nonsense.
                    elif (not (point_in_box (p,b)))             [ b ];                                                                                  # Nothing to do if p is not even in b.
                    elif (row == p.row and col == p.col)        [ b ];                                                                                  # Nothing to do if p is at origin of b.
                    elif (row == p.row)                                                                                                                 #
                        #
                        [ { row, col,          high, wide =>        (p.col - col) },                                                                    # p is on top row of b (but not at origin) so it splits b into two parts vertically.
                          { row, col => p.col, high, wide => wide - (p.col - col) }
                        ];
                    elif (col == p.col)                                                                                                                 # p is on left column of b (but not at origin) so it splits b into two parts horizontally.
                        #
                        [ { col, row,          wide, high =>        (p.row - row) },
                          { col, row => p.row, wide, high => high - (p.row - row) }
                        ];
                    else                                                                                                                                # Here we have the general case of splitting b into four quadrants with p at origin of lower-right quadrant.
                        [ { col,           row,           wide =>        (p.col - col),  high =>        (p.row - row) },
                          { col => p.col,  row,           wide => wide - (p.col - col),  high =>        (p.row - row) },
                          { col,           row => p.row,  wide =>        (p.col - col),  high => high - (p.row - row) },
                          { col => p.col,  row => p.row,  wide => wide - (p.col - col),  high => high - (p.row - row) }
                        ];
                    fi;

                fun quadsect_boxes (boxes: List(Box), p: Point)                                                                                         # Divide box b into four nonintersecting subboxes with p as origin of lower-right one, if p is in interior of b.
                    =
                    list::cat  (map  {. quadsect_box (#b,p); }  boxes);



                fun subtract_box_b_from_box_a { a: Box, b: Box }                                                                                        # This is intended for light work like EXPOSE event handling with a dozen or
                    =                                                                                                                                   # so rectangles, so we opt for simplicity rather than asymptotic efficiency:
                    {
                        (box_corners b)  -> { upper_left, lower_left, lower_right, upper_right };                                                       # Get corners of 'b'.
                        #
                        a_parts =   {   a_parts = [ a ];                                                                                                # Break 'a' up into sub-boxes such that each sub-box is entirely inside 'b' or entirely outside 'b'.
                                        a_parts = bisect_boxes_horizontally (a_parts, upper_left,  upper_right );
                                        a_parts = bisect_boxes_horizontally (a_parts, lower_left,  lower_right );
                                        a_parts = bisect_boxes_vertically   (a_parts, lower_left,  upper_left  );
                                        a_parts = bisect_boxes_vertically   (a_parts, lower_right, upper_right );
                                        a_parts;
                                    };
                        #
                        list::remove                                                                                                                    # Return all the parts of 'a' which are not inside 'b'.
                            {. box_a_in_box_b { a => #a_part, b }; }
                            a_parts;
                    };

                fun subtract_boxes_b_from_boxes_a { a: List(Box), b => []: List(Box) }
                        =>
                        a;

                    subtract_boxes_b_from_boxes_a { a, b => (b ! rest) }                                                                                # Here we subtract 'b' from all boxes in 'a' and then recursively
                        =>                                                                                                                              # recursively subtract all boxes in 'rest' from the result.
                        subtract_boxes_b_from_boxes_a
                          {
                            a => list::cat   (map   {. subtract_box_b_from_box_a { a => #a, b }; }   a),
                            b => rest
                          };
                 end;
            };

            package line {
                #
                # Find the intersection of two Lines.
                # Return NULL if the lines are parallel.
                #
                fun intersection                                                                                                                        # For background see e.g.  http://en.wikipedia.org/wiki/Line-line_intersection
                    ( (a1 as ({ col=>a1x, row=>a1y } ), a2),
                      (b1 as ({ col=>b1x, row=>b1y } ), b2)
                    )
                    =
                    {
                        my { col=>ax, row=>ay } = point::subtract (a2, a1);
                        my { col=>bx, row=>by } = point::subtract (b2, b1);

                        ax = ebf::from_int ax;          ay = ebf::from_int ay;                                                                          # The following arithmetic can easily exceed 32-bit Int precision so we use Float instead.
                        bx = ebf::from_int bx;          by = ebf::from_int by;

                        a1x = ebf::from_int a1x;        a1y = ebf::from_int a1y;
                        b1x = ebf::from_int b1x;        b1y = ebf::from_int b1y;

                        axby = ax * by;
                        bxay = bx * ay;
                        axbx = ax * bx;
                        ayby = ay * by;

                        if (axby == bxay)
                            #
                            NULL;
                        else 
                            fun solve (p, q)
                                =
                                {
                                   my (p, q)
                                        =
                                        if (q < 0.0)   (-p,-q);
                                        else           ( p, q);
                                        fi;

                                    if (p < 0.0)   -(((-p) + (q / 2.0)) / q);
                                    else            ((( p) + (q / 2.0)) / q);
                                    fi;
                                };

                            col = solve (a1x*bxay - b1x*axby + (b1y - a1y)*axbx, bxay - axby);
                            row = solve (a1y*axby - b1y*bxay + (b1x - a1x)*ayby, axby - bxay);

                            row = ebf::round row;
                            col = ebf::round col;

                            (THE ({ col, row } ));
                        fi;
                    };

                    fun rotate_90_degrees_counterclockwise                                                                                              # For background see   http://en.wikipedia.org/wiki/Rotation_matrix
                        ( a as { col,       row       },                                                                                                # Basically we're multiplying by the usual 2d rotation matrix
                          b as { col=>col', row=>row' }                                                                                                 #       |  cos(a) -sin(a) |
                        )                                                                                                                               #       |  sin(a)  cos(a) |
                        =                                                                                                                               # but for a 90% rotation sin(a)==1 and cos(a)==0
                        {   c = { col => col+(row'-row),                                                                                                # making the computation very simple.
                                  row => row-(col'-col)
                                };
                            (a, c);
                        };
            };


            # Bounding box of a list of points:
            #
            fun bounding_box []
                    =>
                    { col=>0, row=>0, wide=>0, high=>0 };

                bounding_box (({ col, row } ) ! pts)
                    =>
                    bb (col, row, col, row, pts)
                    where
                        fun bb (minx, miny, maxx, maxy, [])
                                => 
                                { col => minx, row => miny, wide => maxx-minx+1, high => maxy-miny+1 };

                            bb (minx, miny, maxx, maxy, ({ col, row } ) ! pts)
                                => 
                               bb (min (minx, col), min (miny, row), max (maxx, col), max (maxy, row), pts);
                        end;
                    end;
            end;


            fun convex_hull (points: List(Point))                                               # http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
                =
                # This is a nice little algorithm by the late                                   # We can check for "turns to the right" by              
                # Alex M Andrew, an interesting fellow working                                  # computing the cross-product of the vectors formed     
                # mostly in early AI, particularly neural nets:                                 # by the last three points, O,A,B. The cross-product    
                #                                                                               # computes the area bounded by the parallelogram        
                #     http://cyberneticzoo.com/tag/dr-alex-andrew/                              # defined by the two vectors:                           
                #                                                                               #                   A ---- B                            
                # The idea: A convex hull is (obviously!) *convex*,                             #                    /  /                               
                # meaning that if we walk counterclockwise around a                             #                 O ----                                
                # convex hull, at each point we will turn left, never                           # The area will be positive for a left turn and         
                # right.                                                                        # negative for a right turn.                            
                #     So if we sort the points in X and then scan                                                                                       
                # linearly through them successively adding each point                          # We could also just compute the slope of each vector   
                # to our candidate convex hull sequence, any time we                            # and verify that the slope always increases, but       
                # find our candidate solution turning to the right,                             # the slope requires a division, which is usually       
                # something is wrong.                                                           # much slower than multiplication on today's hardware.  
                #     We can fix this by, before we add a new point,                            #     (If we multiply out the slope comparison formula  
                # checking to see if doing so will make the convex                              # to eliminate the divisions we arrive back at the      
                # hull candidate sequence turn to the right:  If so,                            # cross-product formula.)                               
                # we drop trailing points off the candidate until the
                # problem goes away.                                                            # This fn is exercised in:
                #     When we're done with our scan, we'll have the                             #     src/lib/x-kit/widget/widget-unit-test.pkg
                # lower half of the convex hull;  repeating the                         
                # the procedure in the opposite direction gives us                              
                # the upper half.  Glue them together and we're done.                           # Note that result from each scan will always contain
                #     We'll never do more work per point than adding                            # the first and last point in the sorted input. Thus
                # it to the candidate solution, computing a cross product,                      # if we naively glue together the two results, the
                # and later removing it again -- all O(1) operations -- so                      # first and last points will be duplicated in the
                # the time to do the scan is O(N).                                              # final result.
                #     The time to sort the points on X is O(N*log(N)),                          #   We fix this simply by dropping the final point
                # so the algorithm as a whole is thus     O(N*log(N)),                          # from each scan result before concatenating them to
                #     log(N) grows so slowly that for all practical                             # produce the final result.
                # purposes O(N*log(N)) == O(N), which is the fastest
                # we can construct a sequence of N points.
                #
                {
                    points =  lms::sort_list_and_drop_duplicates  point::compare_xy  points;    # Do this before next because we might have a lot of duplicate points.
                    #
                    if (list::length points < 3)
                        #
                        points;                                                                 # Convex hull of two points is easy. *grin*
                    else
                        {
                            result1 = build_halfhull ((        points), []);                    # Build left-to-right lower half of convex hull.
                            result2 = build_halfhull ((reverse points), []);                    # Build right-to-left upper half of convex hull.
                            #
                            reverse (tail result2 @ tail result1);                              # 'tail' because otherwise first and last points in input get duplicated.
                        }                                                                       # 'reverse' just to give conventional counter-clockwise result ordering.
                        where                                                                   # Obviously we could fiddle the sort ordering to eliminate the 'reverse'
                                                                                                # if efficiency was critical, but here we stay with simple and intuitive.
                                                                                                # If efficiency was critical we'd probably do it in C with arrays. *grin*

                            nonfix my o ;                                                       # Make 'o' not infix so we can use it as a plain variable.
                            #
                            fun cross_product (o: Point, a: Point, b: Point)                    # 2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product.
                                =                                                               # Returns a positive value if OAB makes a counter-clockwise turn,
                                (a.col - o.col) * (b.row - o.row)                               # negative for clockwise turn, and zero if the points are collinear.
                                -
                                (a.row - o.row) * (b.col - o.col);                              # I'm presuming we're working in pixel coordinates on a screen; if not, these multiplies might overflow.


                            fun drop_kinks (p, r1 as (p1 ! (r2 as (p2 ! rest))))                # If (p ! r1) kinks right, drop points from r1 until it no longer kinks right.
                                    =>                                                          #
                                    if (cross_product (p2, p1, p) <= 0)   drop_kinks (p, r2);   # Drop last result point and loop.
                                    else                                  r1;                   #
                                    fi;                                                         #
                                                                                                #
                                drop_kinks (_, result) =>  result;                              #
                            end;


                            fun build_halfhull ([], result) =>  result;
                                #
                                build_halfhull (input as p ! rest, result)
                                    =>
                                    build_halfhull (rest, p ! (drop_kinks (p, result)));
                            end;
                        end;
                    fi;
                };

            fun point_in_polygon  (_,    []) =>  FALSE;                                         # I don't care much about point-on-boundary and point-on-vertex etc. 
                point_in_polygon  (_,   [_]) =>  FALSE;                                         # If you do, there are lots of wwweb resources to consult like
                point_in_polygon  (_, [_,_]) =>  FALSE;                                         #     http://www.ics.uci.edu/~eppstein/161/960307.html
                #                                                                               # Here I'm satisfied with something simple that handles interior vs exterior points ok.
                point_in_polygon
                    (
                      p:                Point,                                                  # Return TRUE iff this point    
                      polygon:          List(Point)                                             # is inside this polygon.
                    )
                    =>
                    {   polygon = (list::last polygon) ! polygon;                               # Duplicate last point at start so that iterating over all pairs gives us all edges in polygon.
                        #
                        point_in_polygon' (polygon, FALSE);
                    }
                    where
                        fun point_in_polygon' ([],           result) =>  result;                # Shouldn't happen -- we're only called if pointlist is nonempty, in which case we terminate on next line.
                            point_in_polygon' ([ p: Point ], result) =>  result;
                            #
                            point_in_polygon' (p2 ! p1 ! rest, point_is_inside)                 # We're drawing a ray to the right from 'p' and counting crossings.
                                =>
                                if ( (p.row < p1.row) != (p.row < p2.row)                       # If ray crosses this edge...
                                and  (p.col < ( ( p.row - p1.row)
                                              * (p2.col - p1.col)                               # As above, I'm presuming we're working in pixel coordinates on a screen; if not, these multiplies might overflow.
                                              / (p2.row - p1.row)
                                              )
                                              + p1.col
                                   ) )
                                    #
                                     point_in_polygon' (p1 ! rest, not point_is_inside);        # ... then flip our inside/outside flag.
                                else point_in_polygon' (p1 ! rest,     point_is_inside);
                                fi;
                        end;
                    end;                                                                        # This fn is exercised in:
            end;                                                                                #     src/lib/x-kit/widget/widget-unit-test.pkg


            # XXX SUCKO FIXME Remaining stuff belongs in xclient, not stdlib:

            fun site_to_box ({ upperleft => { col, row }, size => { wide, high }, ... }: Window_Site)
                =
                { col, row, wide, high };


            # Validation routines:
            #
            fun valid_point({ col, row } )   =   rc::valid_signed16 col    and   rc::valid_signed16 row;
            fun valid_line ((p1, p2): Line)  =       valid_point    p1     and       valid_point    p2;
            fun valid_size ({ wide, high } ) =   rc::valid16        wide   and   rc::valid16        high;

            fun valid_box ({ col, row, wide, high } )
                =
                rc::valid_signed16 col      and
                rc::valid_signed16 row      and 
                rc::valid16 wide            and
                rc::valid16 high;

            fun valid_arc ({ col, row, wide, high, angle1, angle2 } )
                =
                rc::valid_signed16  col     and
                rc::valid_signed16  row     and 
                rc::valid16  wide           and
                rc::valid16  high           and 
                rc::valid_signed16  angle1  and
                rc::valid_signed16  angle2;

            fun valid_site ({ upperleft, size, border_thickness }: Window_Site)
                = 
                valid_point  upperleft      and
                valid_size   size           and
                rc::valid16  border_thickness;

        end;            # stipulate
    };                  # package geometry2d

end;


Comments and suggestions to: bugs@mythryl.org

PreviousUpNext