commit 5bf486ac815ca204dbc314744c00ae462bb2821a
Author: jan <janp@stanford.edu>
Date:   Tue Mar 15 19:12:39 2016 -0700

    Add all files to repository

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..715503a
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,3 @@
+*.pyc
+__pycache__
+*.idea
diff --git a/README b/README
new file mode 100644
index 0000000..ded6941
--- /dev/null
+++ b/README
@@ -0,0 +1,19 @@
+# Masque README
+
+Masque is a Python module for designing lithography masks.
+
+The general idea is to implement something resembling the GDSII file-format, but
+with some vectorized element types (eg. circles, not just polygons), better support for
+E-beam doses, and the ability to output to multiple formats.
+
+## Installation
+
+Requirements:
+* python 3 (written and tested with 3.5)
+* numpy
+* matplotlib (optional, used for visualization functions)
+* python-gdsii (optional, used for gdsii i/o)
+* svgwrite (optional, used for svg output)
+
+
+
diff --git a/masque/__init__.py b/masque/__init__.py
new file mode 100644
index 0000000..25652dc
--- /dev/null
+++ b/masque/__init__.py
@@ -0,0 +1,32 @@
+"""
+ masque 2D CAD library
+
+ masque is an attempt to make a relatively small library for designing lithography
+  masks. The general idea is to implement something resembling the GDSII file-format, but
+  with some vectorized element types (eg. circles, not just polygons), better support for
+  E-beam doses, and the ability to output to multiple formats.
+
+ Pattern is a basic object containing a 2D lithography mask, composed of a list of Shape
+  objects and a list of SubPattern objects.
+
+ SubPattern provides basic support for nesting Pattern objects within each other, by adding
+  offset, rotation, scaling, and other such properties to a Pattern reference.
+
+ Note that the methods for these classes try to avoid copying wherever possible, so unless
+  otherwise noted, assume that arguments are stored by-reference.
+
+
+ Dependencies:
+    - numpy
+    - matplotlib    [Pattern.visualize(...)]
+    - python-gdsii  [masque.file.gdsii]
+    - svgwrite      [masque.file.svgwrite]
+"""
+
+from .error import PatternError
+from .shapes import Shape
+from .subpattern import SubPattern
+from .pattern import Pattern
+
+
+__author__ = 'Jan Petykiewicz'
diff --git a/masque/error.py b/masque/error.py
new file mode 100644
index 0000000..8a67b6e
--- /dev/null
+++ b/masque/error.py
@@ -0,0 +1,9 @@
+class PatternError(Exception):
+    """
+    Simple Exception for Pattern objects and their contents
+    """
+    def __init__(self, value):
+        self.value = value
+
+    def __str__(self):
+        return repr(self.value)
diff --git a/masque/file/__init__.py b/masque/file/__init__.py
new file mode 100644
index 0000000..8de11a7
--- /dev/null
+++ b/masque/file/__init__.py
@@ -0,0 +1,3 @@
+"""
+Functions for reading from and writing to various file formats.
+"""
\ No newline at end of file
diff --git a/masque/file/gdsii.py b/masque/file/gdsii.py
new file mode 100644
index 0000000..3cf1e74
--- /dev/null
+++ b/masque/file/gdsii.py
@@ -0,0 +1,171 @@
+"""
+GDSII file format readers and writers
+"""
+
+import gdsii.library
+import gdsii.structure
+import gdsii.elements
+
+from typing import List, Any, Dict
+import numpy
+
+from .utils import mangle_name, make_dose_table
+from .. import Pattern, SubPattern, PatternError
+from ..shapes import Polygon
+from ..utils import rotation_matrix_2d, get_bit, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+def write_dose2dtype(pattern: Pattern,
+                     filename: str,
+                     meters_per_unit: float):
+    """
+    Write a Pattern to a GDSII file, by first calling .polygonize() on it
+     to change the shapes into polygons, and then writing patterns as GDSII
+     structures, polygons as boundary elements, and subpatterns as structure
+     references (sref).
+
+    Note that this function modifies the Pattern.
+
+    It is often a good idea to run pattern.subpatternize() prior to calling this function,
+     especially if calling .polygonize() will result in very many vertices.
+
+    If you want pattern polygonized with non-default arguments, just call pattern.polygonize()
+     prior to calling this function.
+
+    :param pattern: A Pattern to write to file. Modified by this function.
+    :param filename: Filename to write to.
+    :param meters_per_unit: Written into the GDSII file, meters per length unit.
+    """
+    # Create library
+    lib = gdsii.library.Library(version=600,
+                                name='masque-write_dose2dtype'.encode('ASCII'),
+                                logical_unit=1,
+                                physical_unit=meters_per_unit)
+
+    # Polygonize pattern
+    pattern.polygonize()
+
+    # Get a dict of id(pattern) -> pattern
+    patterns_by_id = {**(pattern.referenced_patterns_by_id()), id(pattern): pattern}
+
+    # Get a table of (id(subpat.pattern), written_dose) for each subpattern
+    sd_table = make_dose_table(pattern)
+
+    # Figure out all the unique doses necessary to write this pattern
+    #  This means going through each row in sd_table and adding the dose values needed to write
+    #  that subpattern at that dose level
+    dose_vals = set()
+    for pat_id, pat_dose in sd_table:
+        pat = patterns_by_id[pat_id]
+        [dose_vals.add(shape.dose * pat_dose) for shape in pat.shapes]
+
+    if len(dose_vals) > 256:
+        raise PatternError('Too many dose values: {}, maximum 256 when using dtypes.'.format(len(dose_vals)))
+
+    dose_vals_list = list(dose_vals)
+
+    # Now create a structure for each row in sd_table (ie, each pattern + dose combination)
+    #  and add in any Boundary and SREF elements
+    for pat_id, pat_dose in sd_table:
+        pat = patterns_by_id[pat_id]
+
+        structure = gdsii.structure.Structure(name=mangle_name(pat, pat_dose).encode('ASCII'))
+        lib.append(structure)
+
+        # Add a Boundary element for each shape
+        for shape in pat.shapes:
+            for polygon in shape.to_polygons():
+                data_type = dose_vals_list.index(polygon.dose * pat_dose)
+                xy_open = numpy.round(polygon.vertices + polygon.offset).astype(int)
+                xy_closed = numpy.vstack((xy_open, xy_open[0, :]))
+                structure.append(gdsii.elements.Boundary(layer=polygon.layer,
+                                                         data_type=data_type,
+                                                         xy=xy_closed))
+        # Add an SREF for each subpattern entry
+        #  strans must be set for angle and mag to take effect
+        for subpat in pat.subpatterns:
+            dose_mult = subpat.dose * pat_dose
+            sref = gdsii.elements.SRef(struct_name=mangle_name(subpat.pattern,  dose_mult).encode('ASCII'),
+                                       xy=numpy.round([subpat.offset]).astype(int))
+            sref.strans = 0
+            sref.angle = subpat.rotation
+            sref.mag = subpat.scale
+            structure.append(sref)
+
+    with open(filename, mode='wb') as stream:
+        lib.save(stream)
+
+    return dose_vals_list
+
+
+def read_dtype2dose(filename: str) -> (List[Pattern], Dict[str, Any]):
+    """
+    Read a gdsii file and translate it into a list of Pattern objects. GDSII structures are
+     translated into Pattern objects; boundaries are translated into polygons, and srefs and arefs
+     are translated into SubPattern objects.
+
+    :param filename: Filename specifying a GDSII file to read from.
+    :return: Tuple: (List of Patterns generated GDSII structures, Dict of GDSII library info)
+    """
+
+    with open(filename, mode='rb') as stream:
+        lib = gdsii.library.Library.load(stream)
+
+    library_info = {'name': lib.name.decode('ASCII'),
+                    'physical_unit': lib.physical_unit,
+                    'logical_unit': lib.logical_unit,
+                    }
+
+    def ref_element_to_subpat(element, offset: vector2) -> SubPattern:
+        # Helper function to create a SubPattern from an SREF or AREF. Sets subpat.pattern to None
+        #  and sets the instance attribute .ref_name to the struct_name.
+        #
+        # BUG: Figure out what "absolute" means in the context of elements and if the current
+        #       behavior is correct
+        # BUG: Need to check STRANS bit 0 to handle x-reflection
+        subpat = SubPattern(pattern=None, offset=offset)
+        subpat.ref_name = element.struct_name
+        if element.strans is not None:
+            if element.mag is not None:
+                subpat.scale = element.mag
+                # Bit 13 means absolute scale
+                if get_bit(element.strans, 13):
+                    subpat.offset *= subpat.scale
+            if element.angle is not None:
+                subpat.rotation = element.angle
+                # Bit 14 means absolute rotation
+                if get_bit(element.strans, 14):
+                    subpat.offset = numpy.dot(rotation_matrix_2d(subpat.rotation), subpat.offset)
+        return subpat
+
+    patterns = []
+    for structure in lib:
+        pat = Pattern(name=structure.name.decode('ASCII'))
+        for element in structure:
+            # Switch based on element type:
+            if isinstance(element, gdsii.elements.Boundary):
+                pat.shapes.append(
+                    Polygon(vertices=element.xy[:-1],
+                            dose=element.data_type,
+                            layer=element.layer))
+
+            elif isinstance(element, gdsii.elements.SRef):
+                pat.subpatterns.append(ref_element_to_subpat(element, element.xy))
+
+            elif isinstance(element, gdsii.elements.ARef):
+                for offset in element.xy:
+                    pat.subpatterns.append(ref_element_to_subpat(element, offset))
+        patterns.append(pat)
+
+    # Create a dict of {pattern.name: pattern, ...}, then fix up all subpattern.pattern entries
+    #  according to the subpattern.ref_name (which is deleted after use).
+    patterns_dict = dict(((p.name, p) for p in patterns))
+    for p in patterns_dict.values():
+        for sp in p.subpatterns:
+            sp.pattern = patterns_dict[sp.ref_name.decode('ASCII')]
+            del sp.ref_name
+
+    return patterns_dict, library_info
diff --git a/masque/file/svg.py b/masque/file/svg.py
new file mode 100644
index 0000000..6a28c25
--- /dev/null
+++ b/masque/file/svg.py
@@ -0,0 +1,139 @@
+"""
+SVG file format readers and writers
+"""
+
+import svgwrite
+import numpy
+
+from .utils import mangle_name
+from .. import Pattern
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+def write(pattern: Pattern,
+          filename: str,
+          custom_attributes: bool=False):
+    """
+    Write a Pattern to an SVG file, by first calling .polygonize() on it
+     to change the shapes into polygons, and then writing patterns as SVG
+     groups (<g>, inside <defs>), polygons as paths (<path>), and subpatterns
+     as <use> elements.
+
+    Note that this function modifies the Pattern.
+
+    If custom_attributes is True, non-standard pattern_layer and pattern_dose attributes
+     are written to the relevant elements.
+
+    It is often a good idea to run pattern.subpatternize() on pattern prior to
+     calling this function, especially if calling .polygonize() will result in very
+     many vertices.
+
+    If you want pattern polygonized with non-default arguments, just call pattern.polygonize()
+     prior to calling this function.
+
+    :param pattern: Pattern to write to file. Modified by this function.
+    :param filename: Filename to write to.
+    :param custom_attributes: Whether to write non-standard pattern_layer and
+            pattern_dose attributes to the SVG elements.
+    """
+
+    # Polygonize pattern
+    pattern.polygonize()
+
+    [bounds_min, bounds_max] = pattern.get_bounds()
+
+    viewbox = numpy.hstack((bounds_min - 1, (bounds_max - bounds_min) + 2))
+    viewbox_string = '{:g} {:g} {:g} {:g}'.format(*viewbox)
+
+    # Create file
+    svg = svgwrite.Drawing(filename, profile='full', viewBox=viewbox_string,
+                           debug=(not custom_attributes))
+
+    # Get a dict of id(pattern) -> pattern
+    patterns_by_id = {**(pattern.referenced_patterns_by_id()), id(pattern): pattern}
+
+    # Now create a group for each row in sd_table (ie, each pattern + dose combination)
+    #  and add in any Boundary and Use elements
+    for pat in patterns_by_id.values():
+        svg_group = svg.g(id=mangle_name(pat), fill='blue', stroke='red')
+
+        for shape in pat.shapes:
+            for polygon in shape.to_polygons():
+                path_spec = poly2path(polygon.vertices + polygon.offset)
+
+                path = svg.path(d=path_spec)
+                if custom_attributes:
+                    path['pattern_layer'] = polygon.layer
+                    path['pattern_dose'] = polygon.dose
+
+                svg_group.add(path)
+
+        for subpat in pat.subpatterns:
+            transform = 'scale({:g}) rotate({:g}) translate({:g},{:g})'.format(
+                subpat.scale, subpat.rotation, subpat.offset[0], subpat.offset[1])
+            use = svg.use(href='#' + mangle_name(subpat.pattern), transform=transform)
+            if custom_attributes:
+                use['pattern_dose'] = subpat.dose
+            svg_group.add(use)
+
+        svg.defs.add(svg_group)
+    svg.add(svg.use(href='#' + mangle_name(pattern)))
+    svg.save()
+
+
+def write_inverted(pattern: Pattern, filename: str):
+    """
+    Write an inverted Pattern to an SVG file, by first calling .polygonize() and
+     .flatten() on it to change the shapes into polygons, then drawing a bounding
+     box and drawing the polygons with reverse vertex order inside it, all within
+     one <path> element.
+
+    Note that this function modifies the Pattern.
+
+    If you want pattern polygonized with non-default arguments, just call pattern.polygonize()
+     prior to calling this function.
+
+    :param pattern: Pattern to write to file. Modified by this function.
+    :param filename: Filename to write to.
+    """
+    # Polygonize and flatten pattern
+    pattern.polygonize().flatten()
+
+    [bounds_min, bounds_max] = pattern.get_bounds()
+
+    viewbox = numpy.hstack((bounds_min - 1, (bounds_max - bounds_min) + 2))
+    viewbox_string = '{:g} {:g} {:g} {:g}'.format(*viewbox)
+
+    # Create file
+    svg = svgwrite.Drawing(filename, profile='full', viewBox=viewbox_string)
+
+    # Draw bounding box
+    slab_edge = [[bounds_min[0] - 1, bounds_max[1] + 1],
+                 [bounds_max[0] + 1, bounds_max[1] + 1],
+                 [bounds_max[0] + 1, bounds_min[1] - 1],
+                 [bounds_min[0] - 1, bounds_min[1] - 1]]
+    path_spec = poly2path(slab_edge)
+
+    # Draw polygons with reversed vertex order
+    for shape in pattern.shapes:
+        for polygon in shape.to_polygons():
+            path_spec += poly2path(polygon.vertices[::-1] + polygon.offset)
+
+    svg.add(svg.path(d=path_spec, fill='blue', stroke='red'))
+    svg.save()
+
+
+def poly2path(vertices: numpy.ndarray) -> str:
+    """
+    Create an SVG path string from an Nx2 list of vertices.
+
+    :param vertices: Nx2 array of vertices.
+    :return: SVG path-string.
+    """
+    commands = 'M{:g},{:g} '.format(vertices[0][0], vertices[0][1])
+    for vertex in vertices[1:]:
+        commands += 'L{:g},{:g}'.format(vertex[0], vertex[1])
+    commands += ' Z   '
+    return commands
diff --git a/masque/file/utils.py b/masque/file/utils.py
new file mode 100644
index 0000000..615d455
--- /dev/null
+++ b/masque/file/utils.py
@@ -0,0 +1,41 @@
+"""
+Helper functions for file reading and writing
+"""
+import re
+from typing import Set, Tuple
+
+from masque.pattern import Pattern
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+def mangle_name(pattern: Pattern, dose_multiplier: float=1.0) -> str:
+    """
+    Create a name using pattern.name, id(pattern), and the dose multiplier.
+
+    :param pattern: Pattern whose name we want to mangle.
+    :param dose_multiplier: Dose multiplier to mangle with.
+    :return: Mangled name.
+    """
+    expression = re.compile('[^A-Za-z0-9_\?\$]')
+    sanitized_name = expression.sub('_', pattern.name)
+    full_name = '{}_{}_{}'.format(sanitized_name, dose_multiplier, id(pattern))
+    return full_name
+
+
+def make_dose_table(pattern: Pattern, dose_multiplier: float=1.0) -> Set[Tuple[int, float]]:
+    """
+    Create a set containing (id(subpat.pattern), written_dose) for each subpattern
+
+    :param pattern: Source Pattern.
+    :param dose_multiplier: Multiplier for all written_dose entries.
+    :return: {(id(subpat.pattern), written_dose), ...}
+    """
+    dose_table = {(id(pattern), dose_multiplier)}
+    for subpat in pattern.subpatterns:
+        subpat_dose_entry = (id(subpat.pattern), subpat.dose * dose_multiplier)
+        if subpat_dose_entry not in dose_table:
+            subpat_dose_table = make_dose_table(subpat.pattern, subpat.dose * dose_multiplier)
+            dose_table = dose_table.union(subpat_dose_table)
+    return dose_table
diff --git a/masque/pattern.py b/masque/pattern.py
new file mode 100644
index 0000000..b187b9e
--- /dev/null
+++ b/masque/pattern.py
@@ -0,0 +1,412 @@
+"""
+ Base object for containing a lithography mask.
+"""
+
+from typing import List, Callable, Tuple, Dict, Union
+import copy
+import itertools
+import pickle
+from collections import defaultdict
+
+import numpy
+# .visualize imports matplotlib and matplotlib.collections
+
+from .subpattern import SubPattern
+from .shapes import Shape, Polygon
+from .utils import rotation_matrix_2d, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+class Pattern(object):
+    """
+    2D layout consisting of some set of shapes and references to other Pattern objects
+     (via SubPattern). Shapes are assumed to inherit from .shapes.Shape or provide equivalent
+     functions.
+
+    :var shapes: List of all shapes in this Pattern. Elements in this list are assumed to inherit
+            from Shape or provide equivalent functions.
+    :var subpatterns: List of all SubPattern objects in this Pattern. Multiple SubPattern objects
+            may reference the same Pattern object.
+    :var name: An identifier for this object. Not necessarily unique.
+    """
+
+    shapes = List[Shape]
+    subpatterns = List[SubPattern]
+    name = str
+
+    def __init__(self,
+                 shapes: List[Shape]=(),
+                 subpatterns: List[SubPattern]=(),
+                 name: str='',
+                 ):
+        """
+        Basic init; arguments get assigned to member variables.
+         Non-list inputs for shapes and subpatterns get converted to lists.
+
+        :param shapes: Initial shapes in the Pattern
+        :param subpatterns: Initial subpatterns in the Pattern
+        :param name: An identifier for the Pattern
+        """
+        if isinstance(shapes, list):
+            self.shapes = shapes
+        else:
+            self.shapes = list(shapes)
+
+        if isinstance(subpatterns, list):
+            self.subpatterns = subpatterns
+        else:
+            self.subpatterns = list(subpatterns)
+
+        self.name = name
+
+    def append(self, other_pattern: 'Pattern') -> 'Pattern':
+        """
+        Appends all shapes and subpatterns from other_pattern to self's shapes and subpatterns.
+
+        :param other_pattern: The Pattern to append
+        :return: self
+        """
+        self.subpatterns += other_pattern.subpatterns
+        self.shapes += other_pattern.shapes
+        return self
+
+    def subset(self,
+               shapes_func: Callable[[Shape], bool]=None,
+               subpatterns_func: Callable[[SubPattern], bool]=None,
+               ) -> 'Pattern':
+        """
+        Returns a Pattern containing only the shapes and subpatterns for which shapes_func or
+         subpatterns_func returns True.
+        Self is _not_ altered, but shapes and subpatterns are _not_ copied.
+
+        :param shapes_func: Given a shape, returns a boolean denoting whether the shape is a member
+                        of the subset
+        :param subpatterns_func: Given a subpattern, returns a boolean denoting if it is a member
+                        of the subset
+        :return: A Pattern containing all the shapes and subpatterns for which the parameter
+                        functions return True
+        """
+        pat = Pattern()
+        if shapes_func is not None:
+            pat.shapes = [s for s in self.shapes if shapes_func(s)]
+        if subpatterns_func is not None:
+            pat.subpatterns = [s for s in self.subpatterns if subpatterns_func(s)]
+        return pat
+
+    def polygonize(self,
+                   poly_num_points: int=None,
+                   poly_max_arclen: float=None
+                   ) -> 'Pattern':
+        """
+        Calls .to_polygons(...) on all the shapes in this Pattern and any referenced patterns,
+         replacing them with the returned polygons.
+        Arguments are passed directly to shape.to_polygons(...).
+
+        :param poly_num_points: Number of points to use for each polygon. Can be overridden by
+                     poly_max_arclen if that results in more points. Optional, defaults to shapes'
+                      internal defaults.
+        :param poly_max_arclen: Maximum arclength which can be approximated by a single line
+                     segment. Optional, defaults to shapes' internal defaults.
+        :return: self
+        """
+        old_shapes = self.shapes
+        self.shapes = list(itertools.chain.from_iterable(
+                        (shape.to_polygons(poly_num_points, poly_max_arclen)
+                         for shape in old_shapes)))
+        for subpat in self.subpatterns:
+            subpat.pattern.polygonize(poly_num_points, poly_max_arclen)
+        return self
+
+    def subpatternize(self,
+                      recursive: bool=True,
+                      norm_value: int=1e6,
+                      exclude_types: Tuple[Shape]=(Polygon,)
+                      ) -> 'Pattern':
+        """
+        Iterates through this Pattern and all referenced Patterns. Within each Pattern, it iterates
+         over all shapes, calling .normalized_form(norm_value) on them to retrieve a scale-,
+         offset-, dose-, and rotation-independent form. Each shape whose normalized form appears
+         more than once is removed and re-added using subpattern objects referencing a newly-created
+         Pattern containing only the normalized form of the shape.
+
+        Note that the default norm_value was chosen to give a reasonable precision when converting
+         to GDSII, which uses integer values for pixel coordinates.
+
+        :param recursive: Whether to call recursively on self's subpatterns. Default True.
+        :param norm_value: Passed to shape.normalized_form(norm_value). Default 1e6 (see function
+                     note about GDSII)
+        :param exclude_types: Shape types passed in this argument are always left untouched, for
+                speed or convenience. Default: (Shapes.Polygon,)
+        :return: self
+        """
+
+        if exclude_types is None:
+            exclude_types = ()
+
+        if recursive:
+            for subpat in self.subpatterns:
+                subpat.pattern.subpatternize(recursive=True,
+                                             norm_value=norm_value,
+                                             exclude_types=exclude_types)
+
+        # Create a dict which uses the label tuple from .normalized_form() as a key, and which
+        #  stores (function_to_create_normalized_shape, [(index_in_shapes, values), ...]), where
+        #  values are the (offset, scale, rotation, dose) values as calculated by .normalized_form()
+        shape_table = defaultdict(lambda: [None, list()])
+        for i, shape in enumerate(self.shapes):
+            if not any((isinstance(shape, t) for t in exclude_types)):
+                label, values, func = shape.normalized_form(norm_value)
+                shape_table[label][0] = func
+                shape_table[label][1].append((i, values))
+
+        # Iterate over the normalized shapes in the table. If any normalized shape occurs more than
+        #  once, create a Pattern holding a normalized shape object, and add self.subpatterns
+        #  entries for each occurrence in self. Also, note down that we should delete the
+        #  self.shapes entries for which we made SubPatterns.
+        shapes_to_remove = []
+        for label in shape_table:
+            if len(shape_table[label][1]) > 1:
+                shape = shape_table[label][0]()
+                pat = Pattern(shapes=[shape])
+
+                for i, values in shape_table[label][1]:
+                    (offset, scale, rotation, dose) = values
+                    subpat = SubPattern(pattern=pat, offset=offset, scale=scale,
+                                        rotation=rotation, dose=dose)
+                    self.subpatterns.append(subpat)
+                    shapes_to_remove.append(i)
+
+        # Remove any shapes for which we have created subpatterns.
+        for i in sorted(shapes_to_remove, reverse=True):
+            del self.shapes[i]
+
+        return self
+
+    def as_polygons(self) -> List[numpy.ndarray]:
+        """
+        Represents the pattern as a list of polygons.
+
+        Deep-copies the pattern, then calls .polygonize() and .flatten() on the copy in order to
+         generate the list of polygons.
+
+        :return: A list of (Ni, 2) numpy.ndarrays specifying vertices of the polygons. Each ndarray
+            is of the form [[x0, y0], [x1, y1],...].
+        """
+        pat = copy.deepcopy(self).polygonize().flatten()
+        return [shape.vertices + shape.offset for shape in pat.shapes]
+
+    def referenced_patterns_by_id(self) -> Dict[int, 'Pattern']:
+        """
+        Create a dictionary of {id(pat): pat} for all Pattern objects referenced by this
+         Pattern (operates recursively on all referenced Patterns as well)
+
+        :return: Dictionary of {id(pat): pat} for all referenced Pattern objects
+        """
+        ids = {}
+        for subpat in self.subpatterns:
+            if id(subpat.pattern) not in ids:
+                ids[id(subpat.pattern)] = subpat.pattern
+                ids.update(subpat.pattern.referenced_patterns_by_id())
+        return ids
+
+    def get_bounds(self) -> Union[numpy.ndarray, None]:
+        """
+        Return a numpy.ndarray containing [[x_min, y_min], [x_max, y_max]], corresponding to the
+         extent of the Pattern's contents in each dimension.
+        Returns None if the Pattern is empty.
+
+        :return: [[x_min, y_min], [x_max, y_max]] or None
+        """
+        entries = self.shapes + self.subpatterns
+        if not entries:
+            return None
+
+        init_bounds = entries[0].get_bounds()
+        min_bounds = init_bounds[0, :]
+        max_bounds = init_bounds[1, :]
+        for entry in entries[1:]:
+            bounds = entry.get_bounds()
+            min_bounds = numpy.minimum(min_bounds, bounds[0, :])
+            max_bounds = numpy.maximum(max_bounds, bounds[1, :])
+        return numpy.vstack((min_bounds, max_bounds))
+
+    def flatten(self) -> 'Pattern':
+        """
+        Removes all subpatterns and adds equivalent shapes.
+
+        :return: self
+        """
+        subpatterns = copy.deepcopy(self.subpatterns)
+        self.subpatterns = []
+        for subpat in subpatterns:
+            subpat.pattern.flatten()
+            self.shapes += subpat.as_pattern().shapes
+        return self
+
+    def translate_elements(self, offset: vector2) -> 'Pattern':
+        """
+        Translates all shapes and subpatterns by the given offset.
+
+        :param offset: Offset to translate by
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.translate(offset)
+        return self
+
+    def scale_elements(self, scale: float) -> 'Pattern':
+        """"
+        Scales all shapes and subpatterns by the given value.
+
+        :param scale: value to scale by
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.scale(scale)
+        return self
+
+    def scale_by(self, c: float) -> 'Pattern':
+        """
+        Scale this Pattern by the given value
+         (all shapes and subpatterns and their offsets are scaled)
+
+        :param c: value to scale by
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.offset *= c
+            entry.scale_by(c)
+        return self
+
+    def rotate_around(self, pivot: vector2, rotation: float) -> 'Pattern':
+        """
+        Rotate the Pattern around the a location.
+
+        :param pivot: Location to rotate around
+        :param rotation: Angle to rotate by (counter-clockwise, radians)
+        :return: self
+        """
+        pivot = numpy.array(pivot)
+        self.translate_elements(-pivot)
+        self.rotate_elements(rotation)
+        self.rotate_element_centers(rotation)
+        self.translate_elements(+pivot)
+        return self
+
+    def rotate_element_centers(self, rotation: float) -> 'Pattern':
+        """
+        Rotate the offsets of all shapes and subpatterns around (0, 0)
+
+        :param rotation: Angle to rotate by (counter-clockwise, radians)
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.offset = numpy.dot(rotation_matrix_2d(rotation), entry.offset)
+        return self
+
+    def rotate_elements(self, rotation: float) -> 'Pattern':
+        """
+        Rotate each shape and subpattern around its center (offset)
+
+        :param rotation: Angle to rotate by (counter-clockwise, radians)
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.rotate(rotation)
+        return self
+
+    def scale_element_doses(self, factor: float) -> 'Pattern':
+        """
+        Multiply all shape and subpattern doses by a factor
+
+        :param factor: Factor to multiply doses by
+        :return: self
+        """
+        for entry in self.shapes + self.subpatterns:
+            entry.dose *= factor
+        return self
+
+    def copy(self) -> 'Pattern':
+        """
+        Return a copy of the Pattern, deep-copying shapes and copying subpattern entries, but not
+         deep-copying any referenced patterns.
+
+        :return: A copy of the current Pattern.
+        """
+        cp = copy.copy(self)
+        cp.shapes = copy.deepcopy(cp.shapes)
+        cp.subpatterns = [copy.copy(subpat) for subpat in cp.subpatterns]
+        return cp
+
+    @staticmethod
+    def load(filename: str) -> 'Pattern':
+        """
+        Load a Pattern from a file
+
+        :param filename: Filename to load from
+        :return: Loaded Pattern
+        """
+        with open(filename, 'rb') as f:
+            tmp_dict = pickle.load(f)
+
+        pattern = Pattern()
+        pattern.__dict__.update(tmp_dict)
+        return pattern
+
+    def save(self, filename: str) -> 'Pattern':
+        """
+        Save the Pattern to a file
+
+        :param filename: Filename to save to
+        :return: self
+        """
+        with open(filename, 'wb') as f:
+            pickle.dump(self.__dict__, f, protocol=2)
+        return self
+
+    def visualize(self,
+                  offset: vector2=(0., 0.),
+                  line_color: str='k',
+                  fill_color: str='none',
+                  overdraw: bool=False):
+        """
+        Draw a picture of the Pattern and wait for the user to inspect it
+
+        Imports matplotlib.
+
+        :param offset: Coordinates to offset by before drawing
+        :param line_color: Outlines are drawn with this color (passed to matplotlib PolyCollection)
+        :param fill_color: Interiors are drawn with this color (passed to matplotlib PolyCollection)
+        :param overdraw: Whether to create a new figure or draw on a pre-existing one
+        """
+        from matplotlib import pyplot
+        import matplotlib.collections
+
+        offset = numpy.array(offset, dtype=float)
+
+        if not overdraw:
+            figure = pyplot.figure()
+            pyplot.axis('equal')
+        else:
+            figure = pyplot.gcf()
+
+        axes = figure.gca()
+
+        polygons = []
+        for shape in self.shapes:
+            polygons += [offset + s.offset + s.vertices for s in shape.to_polygons()]
+
+        mpl_poly_collection = matplotlib.collections.PolyCollection(polygons,
+                                                                    facecolors=fill_color,
+                                                                    edgecolors=line_color)
+        axes.add_collection(mpl_poly_collection)
+        pyplot.axis('equal')
+
+        for subpat in self.subpatterns:
+            subpat.as_pattern().visualize(offset=offset, overdraw=True,
+                                          line_color=line_color, fill_color=fill_color)
+
+        if not overdraw:
+            pyplot.show()
diff --git a/masque/shapes/__init__.py b/masque/shapes/__init__.py
new file mode 100644
index 0000000..cb7e708
--- /dev/null
+++ b/masque/shapes/__init__.py
@@ -0,0 +1,14 @@
+"""
+Shapes for use with the Pattern class, as well as the Shape abstract class from
+ which they are derived.
+"""
+
+from .shape import Shape, normalized_shape_tuple, DEFAULT_POLY_NUM_POINTS
+
+from .polygon import Polygon
+from .circle import Circle
+from .ellipse import Ellipse
+from .arc import Arc
+
+
+
diff --git a/masque/shapes/arc.py b/masque/shapes/arc.py
new file mode 100644
index 0000000..1cba2ef
--- /dev/null
+++ b/masque/shapes/arc.py
@@ -0,0 +1,273 @@
+from typing import List
+import math
+import numpy
+from numpy import pi
+
+from . import Shape, Polygon, normalized_shape_tuple, DEFAULT_POLY_NUM_POINTS
+from .. import PatternError
+from ..utils import is_scalar, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+class Arc(Shape):
+    """
+    An elliptical arc, formed by cutting off an elliptical ring with two rays which exit from its
+     center. It has a position, two radii, a start and stop angle, a rotation, and a width.
+
+    The radii define an ellipse; the ring is formed with radii +/- width/2.
+    The rotation gives the angle from x-axis, counterclockwise, to the first (x) radius.
+    The start and stop angle are measure counterclockwise from the first (x) radius.
+    """
+
+    _radii = None           # type: numpy.ndarray
+    _angles = None          # type: numpy.ndarray
+    _width = 1.0            # type: float
+    _rotation = 0.0         # type: float
+
+    # Defaults for to_polygons
+    poly_num_points = DEFAULT_POLY_NUM_POINTS   # type: int
+    poly_max_arclen = None                      # type: float
+
+    # radius properties
+    @property
+    def radii(self) -> numpy.ndarray:
+        """
+        Return the radii [rx, ry]
+
+        :return: [rx, ry]
+        """
+        return self.radii
+
+    @radii.setter
+    def radii(self, val: vector2):
+        val = numpy.array(val, dtype=float).flatten()
+        if not val.size == 2:
+            raise PatternError('Radii must have length 2')
+        if not val.min() >= 0:
+            raise PatternError('Radii must be non-negative')
+        self.radii = val
+
+    @property
+    def radius_x(self) -> float:
+        return self.radii[0]
+
+    @radius_x.setter
+    def radius_x(self, val: float):
+        if not val >= 0:
+            raise PatternError('Radius must be non-negative')
+        self.radii[0] = val
+
+    @property
+    def radius_y(self) -> float:
+        return self.radii[1]
+
+    @radius_y.setter
+    def radius_y(self, val: float):
+        if not val >= 0:
+            raise PatternError('Radius must be non-negative')
+        self.radii[1] = val
+
+    # arc start/stop angle properties
+    @property
+    def angles(self) -> vector2:
+        """
+        Return the start and stop angles [a_start, a_stop].
+        Angles are measured from x-axis after rotation, and are stored mod 2*pi
+
+        :return: [a_start, a_stop]
+        """
+        return self._angles
+
+    @angles.setter
+    def angles(self, val: vector2):
+        val = numpy.array(val, dtype=float).flatten()
+        if not val.size == 2:
+            raise PatternError('Angles must have length 2')
+        angles = val % (2 * pi)
+        if angles[0] > pi:
+            self.rotation += pi
+            angles -= pi
+        self._angles = angles
+
+    @property
+    def start_angle(self) -> float:
+        return self.angles[0]
+
+    @start_angle.setter
+    def start_angle(self, val: float):
+        self.angles[0] = val % (2 * pi)
+
+    @property
+    def stop_angle(self) -> float:
+        return self.angles[1]
+
+    @stop_angle.setter
+    def stop_angle(self, val: float):
+        self.angles[1] = val % (2 * pi)
+
+    # Rotation property
+    @property
+    def rotation(self) -> float:
+        """
+        Rotation of radius_x from x_axis, counterclockwise, in radians. Stored mod 2*pi
+
+        :return: rotation counterclockwise in radians
+        """
+        return self._rotation
+
+    @rotation.setter
+    def rotation(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Rotation must be a scalar')
+        self._rotation = val % (2 * pi)
+
+    # Width
+    @property
+    def width(self) -> float:
+        """
+        Width of the arc (difference between inner and outer radii)
+
+        :return: width
+        """
+        return self._width
+
+    @width.setter
+    def width(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Width must be a scalar')
+        if not val > 0:
+            raise PatternError('Width must be positive')
+        self._width = val
+
+    def __init__(self,
+                 radii: vector2,
+                 angles: vector2,
+                 rotation: float=0,
+                 poly_num_points: int=DEFAULT_POLY_NUM_POINTS,
+                 poly_max_arclen: float=None,
+                 offset: vector2=(0.0, 0.0),
+                 layer: int=0,
+                 dose: float=1.0):
+        self.offset = offset
+        self.layer = layer
+        self.dose = dose
+        self.radii = radii
+        self.angles = angles
+        self.rotation = rotation
+        self.poly_num_points = poly_num_points
+        self.poly_max_arclen = poly_max_arclen
+
+    def to_polygons(self, poly_num_points: int=None, poly_max_arclen: float=None) -> List[Polygon]:
+        if poly_num_points is None:
+            poly_num_points = self.poly_num_points
+        if poly_max_arclen is None:
+            poly_max_arclen = self.poly_max_arclen
+
+        if (poly_num_points is None) and (poly_max_arclen is None):
+            raise PatternError('Max number of points and arclength left unspecified' +
+                               ' (default was also overridden)')
+
+        rxy = self.radii
+        ang = self.angles
+
+        # Approximate perimeter
+        # Ramanujan, S., "Modular Equations and Approximations to ,"
+        #  Quart. J. Pure. Appl. Math., vol. 45 (1913-1914), pp. 350-372
+        h = ((rxy[1] - rxy[0]) / rxy.sum()) ** 2
+        ellipse_perimeter = pi * rxy.sum() * (1 + 3 * h / (10 + math.sqrt(4 - 3 * h)))
+        perimeter = abs(ang[0] - ang[1]) / (2 * pi) * ellipse_perimeter
+
+        n = []
+        if poly_num_points is not None:
+            n += [poly_num_points]
+        if poly_max_arclen is not None:
+            n += [perimeter / poly_max_arclen]
+        thetas = numpy.linspace(2 * pi, 0, max(n), endpoint=False)
+
+        sin_th, cos_th = (numpy.sin(thetas), numpy.cos(thetas))
+        wh = self.width / 2.0
+
+        xs1 = (rxy[0] + wh) * cos_th - (rxy[1] + wh) * sin_th
+        ys1 = (rxy[0] + wh) * cos_th - (rxy[1] + wh) * sin_th
+        xs2 = (rxy[0] - wh) * cos_th - (rxy[1] - wh) * sin_th
+        ys2 = (rxy[0] - wh) * cos_th - (rxy[1] - wh) * sin_th
+
+        xs = numpy.hstack((xs1, xs2[::-1]))
+        ys = numpy.hstack((ys1, ys2[::-1]))
+        xys = numpy.vstack((xs, ys)).T
+
+        poly = Polygon(xys, dose=self.dose, layer=self.layer, offset=self.offset)
+        poly.rotate(self.rotation)
+        return [poly]
+
+    def get_bounds(self) -> numpy.ndarray:
+        a = self.angles - 0.5 * pi
+
+        mins = []
+        maxs = []
+        for sgn in (+1, -1):
+            wh = sgn * self.width/2
+            rx = self.radius_x + wh
+            ry = self.radius_y + wh
+
+            sin_r = numpy.sin(self.rotation)
+            cos_r = numpy.cos(self.rotation)
+            tan_r = numpy.tan(self.rotation)
+            sin_a = numpy.sin(a)
+            cos_a = numpy.cos(a)
+
+            xpt = numpy.arctan(-ry / rx * tan_r)
+            ypt = numpy.arctan(+ry / rx / tan_r)
+            xnt = numpy.arcsin(numpy.sin(xpt - pi))
+            ynt = numpy.arcsin(numpy.sin(ypt - pi))
+
+            xr = numpy.sqrt((rx * cos_r) ** 2 + (ry * sin_r) ** 2)
+            yr = numpy.sqrt((rx * sin_r) ** 2 + (ry * cos_r) ** 2)
+
+            xn, xp = sorted(rx * cos_r * cos_a - ry * sin_r * sin_a)
+            yn, yp = sorted(rx * sin_r * cos_a - ry * cos_r * sin_a)
+
+            if min(a) < xpt < max(a):
+                xp = xr
+
+            if min(a) < xnt < max(a):
+                xn = -xr
+
+            if min(a) < ypt < max(a):
+                yp = yr
+
+            if min(a) < ynt < max(a):
+                yn = -yr
+
+            mins.append([xn, yn])
+            maxs.append([xp, yp])
+        return numpy.vstack((numpy.min(mins, axis=0) + self.offset,
+                             numpy.max(maxs, axis=0) + self.offset))
+
+    def rotate(self, theta: float) -> 'Arc':
+        self.rotation += theta
+        return self
+
+    def scale_by(self, c: float) -> 'Arc':
+        self.radii *= c
+        self.width *= c
+        return self
+
+    def normalized_form(self, norm_value: float) -> normalized_shape_tuple:
+        if self.radius_x < self.radius_y:
+            radii = self.radii / self.radius_x
+            scale = self.radius_x
+            rotation = self.rotation
+            angles = self.angles
+        else:  # rotate by 90 degrees and swap radii
+            radii = self.radii[::-1] / self.radius_y
+            scale = self.radius_y
+            rotation = self.rotation + pi / 2
+            angles = self.angles - pi / 2
+        return (type(self), radii, angles, self.layer), \
+               (self.offset, scale/norm_value, rotation, self.dose), \
+               lambda: Arc(radii=radii*norm_value, angles=angles, layer=self.layer)
+
+
diff --git a/masque/shapes/circle.py b/masque/shapes/circle.py
new file mode 100644
index 0000000..f1ddb38
--- /dev/null
+++ b/masque/shapes/circle.py
@@ -0,0 +1,95 @@
+from typing import List
+import numpy
+from numpy import pi
+
+from . import Shape, Polygon, normalized_shape_tuple, DEFAULT_POLY_NUM_POINTS
+from .. import PatternError
+from ..utils import is_scalar, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+class Circle(Shape):
+    """
+    A circle, which has a position and radius.
+    """
+
+    _radius = None                              # type: float
+
+    # Defaults for to_polygons
+    poly_num_points = DEFAULT_POLY_NUM_POINTS   # type: int
+    poly_max_arclen = None                      # type: float
+
+    # radius property
+    @property
+    def radius(self) -> float:
+        """
+        Circle's radius (float, >= 0)
+
+        :return: radius
+        """
+        return self._radius
+
+    @radius.setter
+    def radius(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Radius must be a scalar')
+        if not val >= 0:
+            raise PatternError('Radius must be non-negative')
+        self._radius = val
+
+    def __init__(self,
+                 radius: float,
+                 poly_num_points: int=DEFAULT_POLY_NUM_POINTS,
+                 poly_max_arclen: float=None,
+                 offset: vector2=(0.0, 0.0),
+                 layer: int=0,
+                 dose: float=1.0):
+        self.offset = numpy.array(offset, dtype=float)
+        self.layer = layer
+        self.dose = dose
+        self.radius = radius
+        self.poly_num_points = poly_num_points
+        self.poly_max_arclen = poly_max_arclen
+
+    def to_polygons(self, poly_num_points: int=None, poly_max_arclen: float=None) -> List[Polygon]:
+        if poly_num_points is None:
+            poly_num_points = self.poly_num_points
+        if poly_max_arclen is None:
+            poly_max_arclen = self.poly_max_arclen
+
+        if (poly_num_points is None) and (poly_max_arclen is None):
+            raise PatternError('Number of points and arclength left '
+                               'unspecified (default was also overridden)')
+
+        n = []
+        if poly_num_points is not None:
+            n += [poly_num_points]
+        if poly_max_arclen is not None:
+            n += [2 * pi * self.radius / poly_max_arclen]
+        thetas = numpy.linspace(2 * pi, 0, max(n), endpoint=False)
+        xs = numpy.cos(thetas) * self.radius
+        ys = numpy.sin(thetas) * self.radius
+        xys = numpy.vstack((xs, ys)).T
+
+        return [Polygon(xys, offset=self.offset, dose=self.dose, layer=self.layer)]
+
+    def get_bounds(self) -> numpy.ndarray:
+        return numpy.vstack((self.offset - self.radius,
+                             self.offset + self.radius))
+
+    def rotate(self, theta: float) -> 'Circle':
+        return self
+
+    def scale_by(self, c: float) -> 'Circle':
+        self.radius *= c
+        return self
+
+    def normalized_form(self, norm_value) -> normalized_shape_tuple:
+        rotation = 0.0
+        magnitude = self.radius / norm_value
+        return (type(self), self.layer), \
+               (self.offset, magnitude, rotation, self.dose), \
+               lambda: Circle(radius=norm_value, layer=self.layer)
+
diff --git a/masque/shapes/ellipse.py b/masque/shapes/ellipse.py
new file mode 100644
index 0000000..a9d22e7
--- /dev/null
+++ b/masque/shapes/ellipse.py
@@ -0,0 +1,161 @@
+from typing import List
+import math
+import numpy
+from numpy import pi
+
+from . import Shape, Polygon, normalized_shape_tuple, DEFAULT_POLY_NUM_POINTS
+from .. import PatternError
+from ..utils import is_scalar, rotation_matrix_2d, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+class Ellipse(Shape):
+    """
+    An ellipse, which has a position, two radii, and a rotation.
+    The rotation gives the angle from x-axis, counterclockwise, to the first (x) radius.
+    """
+
+    _radii = None       # type: numpy.ndarray
+    _rotation = 0.0     # type: float
+
+    # Defaults for to_polygons
+    poly_num_points = DEFAULT_POLY_NUM_POINTS   # type: int
+    poly_max_arclen = None                      # type: float
+
+    # radius properties
+    @property
+    def radii(self) -> numpy.ndarray:
+        """
+        Return the radii [rx, ry]
+
+        :return: [rx, ry]
+        """
+        return self.radii
+
+    @radii.setter
+    def radii(self, val: vector2):
+        val = numpy.array(val).flatten()
+        if not val.size == 2:
+            raise PatternError('Radii must have length 2')
+        if not val.min() >= 0:
+            raise PatternError('Radii must be non-negative')
+        self.radii = val
+
+    @property
+    def radius_x(self) -> float:
+        return self.radii[0]
+
+    @radius_x.setter
+    def radius_x(self, val: float):
+        if not val >= 0:
+            raise PatternError('Radius must be non-negative')
+        self.radii[0] = val
+
+    @property
+    def radius_y(self) -> float:
+        return self.radii[1]
+
+    @radius_y.setter
+    def radius_y(self, val: float):
+        if not val >= 0:
+            raise PatternError('Radius must be non-negative')
+        self.radii[1] = val
+
+    # Rotation property
+    @property
+    def rotation(self) -> float:
+        """
+        Rotation of rx from the x axis. Uses the interval [0, pi) in radians (counterclockwise
+         is positive)
+
+        :return: counterclockwise rotation in radians
+        """
+        return self._rotation
+
+    @rotation.setter
+    def rotation(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Rotation must be a scalar')
+        self._rotation = val % pi
+
+    def __init__(self,
+                 radii: vector2,
+                 rotation: float=0,
+                 poly_num_points: int=DEFAULT_POLY_NUM_POINTS,
+                 poly_max_arclen: float=None,
+                 offset: vector2=(0.0, 0.0),
+                 layer: int=0,
+                 dose: float=1.0):
+        self.offset = offset
+        self.layer = layer
+        self.dose = dose
+        self.radii = radii
+        self.rotation = rotation
+        self.poly_num_points = poly_num_points
+        self.poly_max_arclen = poly_max_arclen
+
+    def to_polygons(self,
+                    poly_num_points: int=None,
+                    poly_max_arclen: float=None
+                    ) -> List[Polygon]:
+        if poly_num_points is None:
+            poly_num_points = self.poly_num_points
+        if poly_max_arclen is None:
+            poly_max_arclen = self.poly_max_arclen
+
+        if (poly_num_points is None) and (poly_max_arclen is None):
+            raise PatternError('Number of points and arclength left unspecified'
+                               ' (default was also overridden)')
+
+        rxy = self.radii
+
+        # Approximate perimeter
+        # Ramanujan, S., "Modular Equations and Approximations to ,"
+        #  Quart. J. Pure. Appl. Math., vol. 45 (1913-1914), pp. 350-372
+        h = ((rxy[1] - rxy[0]) / rxy.sum()) ** 2
+        perimeter = pi * rxy.sum() * (1 + 3 * h / (10 + math.sqrt(4 - 3 * h)))
+
+        n = []
+        if poly_num_points is not None:
+            n += [poly_num_points]
+        if poly_max_arclen is not None:
+            n += [perimeter / poly_max_arclen]
+        thetas = numpy.linspace(2 * pi, 0, max(n), endpoint=False)
+
+        sin_th, cos_th = (numpy.sin(thetas), numpy.cos(thetas))
+        xs = rxy[0] * cos_th - rxy[1] * sin_th
+        ys = rxy[0] * sin_th - rxy[1] * cos_th
+        xys = numpy.vstack((xs, ys)).T
+
+        poly = Polygon(xys, dose=self.dose, layer=self.layer, offset=self.offset)
+        poly.rotate(self.rotation)
+        return [poly]
+
+    def get_bounds(self) -> numpy.ndarray:
+        rot_radii = numpy.dot(rotation_matrix_2d(self.rotation), self.radii)
+        return numpy.vstack((self.offset - rot_radii[0],
+                             self.offset + rot_radii[1]))
+
+    def rotate(self, theta: float) -> 'Ellipse':
+        self.rotation += theta
+        return self
+
+    def scale_by(self, c: float) -> 'Ellipse':
+        self.radii *= c
+        return self
+
+    def normalized_form(self, norm_value: float) -> normalized_shape_tuple:
+        if self.radius_x < self.radius_y:
+            radii = self.radii / self.radius_x
+            scale = self.radius_x
+            angle = self.rotation
+        else:
+            radii = self.radii[::-1] / self.radius_y
+            scale = self.radius_y
+            angle = (self.rotation + pi / 2) % pi
+        return (type(self), radii, self.layer), \
+               (self.offset, scale/norm_value, angle, self.dose), \
+               lambda: Ellipse(radii=radii*norm_value, layer=self.layer)
+
diff --git a/masque/shapes/polygon.py b/masque/shapes/polygon.py
new file mode 100644
index 0000000..761c1fa
--- /dev/null
+++ b/masque/shapes/polygon.py
@@ -0,0 +1,174 @@
+from typing import List
+import copy
+import numpy
+from numpy import pi
+
+from . import Shape, normalized_shape_tuple
+from .. import PatternError
+from ..utils import is_scalar, rotation_matrix_2d, vector2
+
+__author__ = 'Jan Petykiewicz'
+
+
+class Polygon(Shape):
+    """
+    A polygon, consisting of a bunch of vertices (Nx2 ndarray) along with an offset.
+
+    A normalized_form(...) is available, but can be quite slow with lots of vertices.
+    """
+    _vertices = None        # type: numpy.ndarray
+
+    # vertices property
+    @property
+    def vertices(self) -> numpy.ndarray:
+        """
+        Vertices of the polygon (Nx2 ndarray: [[x0, y0], [x1, y1], ...]
+
+        :return: vertices
+        """
+        return self._vertices
+
+    @vertices.setter
+    def vertices(self, val: numpy.ndarray):
+        val = numpy.array(val, dtype=float)
+        if val.shape[1] != 2:
+            raise PatternError('Vertices must be an Nx2 array')
+        if val.shape[0] < 3:
+            raise PatternError('Must have at least 3 vertices (Nx2, N>3)')
+        self._vertices = val
+
+    # xs property
+    @property
+    def xs(self) -> numpy.ndarray:
+        """
+        All x vertices in a 1D ndarray
+        """
+        return self.vertices[:, 0]
+
+    @xs.setter
+    def xs(self, val: numpy.ndarray):
+        val = numpy.array(val, dtype=float).flatten()
+        if val.size != self.vertices.shape[0]:
+            raise PatternError('Wrong number of vertices')
+        self.vertices[:, 0] = val
+
+    # ys property
+    @property
+    def ys(self) -> numpy.ndarray:
+        """
+        All y vertices in a 1D ndarray
+        """
+        return self.vertices[:, 1]
+
+    @ys.setter
+    def ys(self, val: numpy.ndarray):
+        val = numpy.array(val, dtype=float).flatten()
+        if val.size != self.vertices.shape[0]:
+            raise PatternError('Wrong number of vertices')
+        self.vertices[:, 1] = val
+
+    def __init__(self,
+                 vertices: numpy.ndarray,
+                 offset: vector2=(0.0, 0.0),
+                 layer: int=0,
+                 dose: float=1.0):
+        self.offset = offset
+        self.layer = layer
+        self.dose = dose
+        self.vertices = vertices
+
+    @staticmethod
+    def square(side_length: float,
+               rotation: float=0.0,
+               offset: vector2=(0.0, 0.0),
+               layer: int=0,
+               dose: float=1.0
+               ) -> 'Polygon':
+        """
+        Draw a square given side_length
+
+        :param side_length: Length of one side
+        :param rotation: Rotation counterclockwise, in radians
+        :param offset: Offset, default (0, 0)
+        :param layer: Layer, default 0
+        :param dose: Dose, default 1.0
+        :return: A Polygon object containing the requested square
+        """
+        norm_square = numpy.array([[-1, -1],
+                                   [-1, +1],
+                                   [+1, +1],
+                                   [+1, -1]], dtype=float)
+        vertices = 0.5 * side_length * norm_square
+        poly = Polygon(vertices, offset, layer, dose)
+        poly.rotate(rotation)
+        return poly
+
+    @staticmethod
+    def rectangle(lx: float,
+                  ly: float,
+                  rotation: float=0,
+                  offset: vector2=(0.0, 0.0),
+                  layer: int=0,
+                  dose: float=1.0
+                  ) -> 'Polygon':
+        """
+        Draw a rectangle with side lengths lx and ly
+
+        :param lx: Length along x (before rotation)
+        :param ly: Length along y (before rotation)
+        :param rotation: Rotation counterclockwise, in radians
+        :param offset: Offset, default (0, 0)
+        :param layer: Layer, default 0
+        :param dose: Dose, default 1.0
+        :return: A Polygon object containing the requested rectangle
+        """
+        vertices = 0.5 * numpy.array([[-lx, -ly],
+                                      [-lx, +ly],
+                                      [+lx, +ly],
+                                      [+lx, -ly]], dtype=float)
+        poly = Polygon(vertices, offset, layer, dose)
+        poly.rotate(rotation)
+        return poly
+
+    def to_polygons(self,
+                    _poly_num_points: int=None,
+                    _poly_max_arclen: float=None,
+                    ) -> List['Polygon']:
+        return [copy.deepcopy(self)]
+
+    def get_bounds(self) -> numpy.ndarray:
+        return numpy.vstack((self.offset + numpy.min(self.vertices, axis=0),
+                             self.offset + numpy.max(self.vertices, axis=0)))
+
+    def rotate(self, theta: float) -> 'Polygon':
+        self.vertices = numpy.dot(rotation_matrix_2d(theta), self.vertices.T).T
+        return self
+
+    def scale_by(self, c: float) -> 'Polygon':
+        self.vertices *= c
+        return self
+
+    def normalized_form(self, norm_value: float) -> normalized_shape_tuple:
+        # Note: this function is going to be pretty slow for many-vertexed polygons, relative to
+        #   other shapes
+        offset = self.vertices.mean(axis=0) + self.offset
+        zeroed_vertices = self.vertices - offset
+
+        scale = zeroed_vertices.std()
+        normed_vertices = zeroed_vertices / scale
+
+        _, _, vertex_axis = numpy.linalg.svd(zeroed_vertices)
+        rotation = numpy.arctan2(vertex_axis[0][1], vertex_axis[0][0]) % (2 * pi)
+        rotated_vertices = numpy.vstack([numpy.dot(rotation_matrix_2d(-rotation), v)
+                                        for v in normed_vertices])
+
+        # Reorder the vertices so that the one with lowest x, then y, comes first.
+        x_min = rotated_vertices[:, 0].argmin()
+        if not is_scalar(x_min):
+            y_min = rotated_vertices[x_min, 1].argmin()
+            x_min = x_min[y_min]
+        reordered_vertices = numpy.roll(rotated_vertices, -x_min, axis=0)
+
+        return (type(self), reordered_vertices.data.tobytes(), self.layer), \
+               (offset, scale/norm_value, rotation, self.dose), \
+               lambda: Polygon(reordered_vertices*norm_value, layer=self.layer)
diff --git a/masque/shapes/shape.py b/masque/shapes/shape.py
new file mode 100644
index 0000000..cacac44
--- /dev/null
+++ b/masque/shapes/shape.py
@@ -0,0 +1,182 @@
+from typing import List, Tuple, Callable
+from abc import ABCMeta, abstractmethod
+import numpy
+
+from .. import PatternError
+from ..utils import is_scalar, rotation_matrix_2d, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+# Type definitions
+normalized_shape_tuple = Tuple[Tuple,
+                               Tuple[numpy.ndarray, float, float, float],
+                               Callable[[], 'Shape']]
+
+# ## Module-wide defaults
+# Default number of points per polygon for shapes
+DEFAULT_POLY_NUM_POINTS = 24
+
+
+class Shape(object, metaclass=ABCMeta):
+    """
+    Abstract class specifying functions common to all shapes.
+    """
+
+    # [x_offset, y_offset]
+    _offset = numpy.array([0.0, 0.0])   # type: numpy.ndarray
+
+    # Layer (integer >= 0)
+    _layer = 0                          # type: int
+
+    # Dose
+    _dose = 1.0                         # type: float
+
+    # --- Abstract methods
+    @abstractmethod
+    def to_polygons(self, num_vertices: int, max_arclen: float) -> List['Polygon']:
+        """
+        Returns a list of polygons which approximate the shape.
+
+        :param num_vertices: Number of points to use for each polygon. Can be overridden by
+                     max_arclen if that results in more points. Optional, defaults to shapes'
+                      internal defaults.
+        :param max_arclen: Maximum arclength which can be approximated by a single line
+                     segment. Optional, defaults to shapes' internal defaults.
+        :return: List of polygons equivalent to the shape
+        """
+        pass
+
+    @abstractmethod
+    def get_bounds(self) -> numpy.ndarray:
+        """
+        Returns [[x_min, y_min], [x_max, y_max]] which specify a minimal bounding box for the shape.
+
+        :return: [[x_min, y_min], [x_max, y_max]]
+        """
+        pass
+
+    @abstractmethod
+    def rotate(self, theta: float) -> 'Shape':
+        """
+        Rotate the shape around its center (0, 0), ignoring its offset.
+
+        :param theta: Angle to rotate by (counterclockwise, radians)
+        :return: self
+        """
+        pass
+
+    @abstractmethod
+    def scale_by(self, c: float) -> 'Shape':
+        """
+        Scale the shape's size (eg. radius, for a circle) by a constant factor.
+
+        :param c: Factor to scale by
+        :return: self
+        """
+        pass
+
+    @abstractmethod
+    def normalized_form(self, norm_value: int) -> normalized_shape_tuple:
+        """
+        Writes the shape in a standardized notation, with offset, scale, rotation, and dose
+         information separated out from the remaining values.
+
+        :param norm_value: This value is used to normalize lengths intrinsic to teh shape;
+                eg. for a circle, the returned magnitude value will be (radius / norm_value), and
+                the returned callable will create a Circle(radius=norm_value, ...). This is useful
+                when you find it important for quantities to remain in a certain range, eg. for
+                GDSII where vertex locations are stored as integers.
+        :return: The returned information takes the form of a 3-element tuple,
+                (intrinsic, extrinsic, constructor). These are further broken down as:
+                extrinsic: ([x_offset, y_offset], scale, rotation, dose)
+                intrinsic: A tuple of basic types containing all information about the instance that
+                            is not contained in 'extrinsic'. Usually, intrinsic[0] == type(self).
+                constructor: A callable (no arguments) which returns an instance of type(self) with
+                            internal state equivalent to 'intrinsic'.
+        """
+        pass
+
+    # ---- Non-abstract properties
+    # offset property
+    @property
+    def offset(self) -> numpy.ndarray:
+        """
+        [x, y] offset
+
+        :return: [x_offset, y_offset]
+        """
+        return self._offset
+
+    @offset.setter
+    def offset(self, val: vector2):
+        if not isinstance(val, numpy.ndarray):
+            val = numpy.array(val, dtype=float)
+
+        if val.size != 2:
+            raise PatternError('Offset must be convertible to size-2 ndarray')
+        self._offset = val.flatten()
+
+    # layer property
+    @property
+    def layer(self) -> int:
+        """
+        Layer number (int, >=0)
+
+        :return: Layer
+        """
+        return self._layer
+
+    @layer.setter
+    def layer(self, val: int):
+        if not isinstance(val, int):
+            raise PatternError('Layer must be an integer')
+        if not val >= 0:
+            raise PatternError('Layer must be non-negative')
+        self._layer = val
+
+    # dose property
+    @property
+    def dose(self) -> float:
+        """
+        Dose (float >= 0)
+
+        :return: Dose value
+        """
+        return self._dose
+
+    @dose.setter
+    def dose(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Dose must be a scalar')
+        if not val >= 0:
+            raise PatternError('Dose must be non-negative')
+        self._dose = val
+
+    # ---- Non-abstract methods
+    def translate(self, offset: vector2) -> 'Shape':
+        """
+        Translate the shape by the given offset
+
+        :param offset: [x_offset, y,offset]
+        :return: self
+        """
+        self.offset += offset
+        return self
+
+    def rotate_around(self, pivot: vector2, rotation: float) -> 'Shape':
+        """
+        Rotate the shape around a point.
+
+        :param pivot: Point (x, y) to rotate around
+        :param rotation: Angle to rotate by (counterclockwise, radians)
+        :return: self
+        """
+        pivot = numpy.array(pivot, dtype=float)
+        self.translate(-pivot)
+        self.rotate(rotation)
+        self.offset = numpy.dot(rotation_matrix_2d(rotation), self.offset)
+        self.translate(+pivot)
+        return self
+
diff --git a/masque/shapes/text.py b/masque/shapes/text.py
new file mode 100644
index 0000000..de2eae1
--- /dev/null
+++ b/masque/shapes/text.py
@@ -0,0 +1,57 @@
+
+#
+# class Text(Shape):
+#     _string = ''
+#     _height = 1.0
+#     _rotation = 0.0
+#     font_path = ''
+#
+#     # vertices property
+#     @property
+#     def string(self):
+#         return self._string
+#
+#     @string.setter
+#     def string(self, val):
+#         self._string = val
+#
+#     # Rotation property
+#     @property
+#     def rotation(self):
+#         return self._rotation
+#
+#     @rotation.setter
+#     def rotation(self, val):
+#         if not is_scalar(val):
+#             raise PatternError('Rotation must be a scalar')
+#         self._rotation = val % (2 * pi)
+#
+#     # Height property
+#     @property
+#     def height(self):
+#         return self._height
+#
+#     @height.setter
+#     def height(self, val):
+#         if not is_scalar(val):
+#             raise PatternError('Height must be a scalar')
+#         self._height = val
+#
+#     def __init__(self, text, height, font_path, rotation=0.0, offset=(0.0, 0.0), layer=0, dose=1.0):
+#         self.offset = offset
+#         self.layer = layer
+#         self.dose = dose
+#         self.text = text
+#         self.height = height
+#         self.rotation = rotation
+#         self.font_path = font_path
+#
+#     def to_polygon(self, _poly_num_points=None, _poly_max_arclen=None):
+#
+#         return copy.deepcopy(self)
+#
+#     def rotate(self, theta):
+#         self.rotation += theta
+#
+#     def scale_by(self, c):
+#         self.height *= c
diff --git a/masque/subpattern.py b/masque/subpattern.py
new file mode 100644
index 0000000..8b389b7
--- /dev/null
+++ b/masque/subpattern.py
@@ -0,0 +1,159 @@
+"""
+ SubPattern provides basic support for nesting Pattern objects within each other, by adding
+  offset, rotation, scaling, and other such properties to the reference.
+"""
+
+from typing import Union
+
+import numpy
+from numpy import pi
+
+from .error import PatternError
+from .utils import is_scalar, rotation_matrix_2d, vector2
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+class SubPattern(object):
+    """
+    SubPattern provides basic support for nesting Pattern objects within each other, by adding
+     offset, rotation, scaling, and associated methods.
+    """
+
+    pattern = None          # type: Pattern
+    _offset = (0.0, 0.0)    # type: numpy.ndarray
+    _rotation = 0.0         # type: float
+    _dose = 1.0             # type: float
+    _scale = 1.0            # type: float
+
+    def __init__(self,
+                 pattern: 'Pattern',
+                 offset: vector2=(0.0, 0.0),
+                 rotation: float=0.0,
+                 dose: float=1.0,
+                 scale: float=1.0):
+        self.pattern = pattern
+        self.offset = offset
+        self.rotation = rotation
+        self.dose = dose
+        self.scale = scale
+
+    # offset property
+    @property
+    def offset(self) -> numpy.ndarray:
+        return self._offset
+
+    @offset.setter
+    def offset(self, val: vector2):
+        if not isinstance(val, numpy.ndarray):
+            val = numpy.array(val, dtype=float)
+
+        if val.size != 2:
+            raise PatternError('Offset must be convertible to size-2 ndarray')
+        self._offset = val.flatten()
+
+    # dose property
+    @property
+    def dose(self) -> float:
+        return self._dose
+
+    @dose.setter
+    def dose(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Dose must be a scalar')
+        if not val >= 0:
+            raise PatternError('Dose must be non-negative')
+        self._dose = val
+
+    # scale property
+    @property
+    def scale(self) -> float:
+        return self._scale
+
+    @scale.setter
+    def scale(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Scale must be a scalar')
+        if not val > 0:
+            raise PatternError('Scale must be positive')
+        self._scale = val
+
+    # Rotation property [ccw]
+    @property
+    def rotation(self) -> float:
+        return self._rotation
+
+    @rotation.setter
+    def rotation(self, val: float):
+        if not is_scalar(val):
+            raise PatternError('Rotation must be a scalar')
+        self._rotation = val % (2 * pi)
+
+    def as_pattern(self) -> 'Pattern':
+        """
+        Returns a copy of self.pattern which has been scaled, rotated, etc. according to this
+         SubPattern's properties.
+        :return: Copy of self.pattern that has been altered to reflect the SubPattern's properties.
+        """
+        pattern = self.pattern.copy()
+        pattern.scale_by(self.scale)
+        pattern.rotate_around((0.0, 0.0), self.rotation)
+        pattern.translate_elements(self.offset)
+        pattern.scale_element_doses(self.dose)
+        return pattern
+
+    def translate(self, offset: vector2) -> 'SubPattern':
+        """
+        Translate by the given offset
+
+        :param offset: Translate by this offset
+        :return: self
+        """
+        self.offset += offset
+        return self
+
+    def rotate_around(self, pivot: vector2, rotation: float) -> 'SubPattern':
+        """
+        Rotate around a point
+
+        :param pivot: Point to rotate around
+        :param rotation: Angle to rotate by (counterclockwise, radians)
+        :return: self
+        """
+        pivot = numpy.array(pivot, dtype=float)
+        self.translate(-pivot)
+        self.offset = numpy.dot(rotation_matrix_2d(rotation), self.offset)
+        self.rotate(rotation)
+        self.translate(+pivot)
+        return self
+
+    def rotate(self, rotation: float) -> 'SubPattern':
+        """
+        Rotate around (0, 0)
+
+        :param rotation: Angle to rotate by (counterclockwise, radians)
+        :return: self
+        """
+        self.rotation += rotation
+        return self
+
+    def get_bounds(self) -> numpy.ndarray or None:
+        """
+        Return a numpy.ndarray containing [[x_min, y_min], [x_max, y_max]], corresponding to the
+         extent of the SubPattern in each dimension.
+        Returns None if the contained Pattern is empty.
+
+        :return: [[x_min, y_min], [x_max, y_max]] or None
+        """
+        return self.as_pattern().get_bounds()
+
+    def scale_by(self, c: float) -> 'SubPattern':
+        """
+        Scale the subpattern by a factor
+
+        :param c: scaling factor
+        """
+        self.scale *= c
+        return self
+
diff --git a/masque/utils.py b/masque/utils.py
new file mode 100644
index 0000000..59e8169
--- /dev/null
+++ b/masque/utils.py
@@ -0,0 +1,41 @@
+"""
+Various helper functions
+"""
+
+from typing import Any, Union, Tuple
+
+import numpy
+
+# Type definitions
+vector2 = Union[numpy.ndarray, Tuple[float, float]]
+
+
+def is_scalar(var: Any) -> bool:
+    """
+    Alias for 'not hasattr(var, "__len__")'
+
+    :param var: Checks if var has a length.
+    """
+    return not hasattr(var, "__len__")
+
+
+def get_bit(bit_string: Any, bit_id: int) -> bool:
+    """
+    Returns true iff bit number 'bit_id' from the right of 'bitstring' is 1
+
+    :param bit_string: st
+    :param bit_id:
+    :return: value of the requested bit (bool)
+    """
+    return bit_string & (1 << bit_id) != 0
+
+
+def rotation_matrix_2d(theta: float) -> numpy.ndarray:
+    """
+    2D rotation matrix for rotating counterclockwise around the origin.
+
+    :param theta: Angle to rotate, in radians
+    :return: rotation matrix
+    """
+    return numpy.array([[numpy.cos(theta), -numpy.sin(theta)],
+                        [numpy.sin(theta), +numpy.cos(theta)]])
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..6189400
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+
+from distutils.core import setup
+
+setup(name='masque',
+      version='0.1',
+      description='Lithography mask library',
+      author='Jan Petykiewicz',
+      author_email='anewusername@gmail.com',
+      url='https://mpxd.net/gogs/jan/masque',
+      packages=['masque'],
+     )
+