diff --git a/model_discretization_api.py b/model_discretization_api.py
new file mode 100644
index 0000000..ed19b89
--- /dev/null
+++ b/model_discretization_api.py
@@ -0,0 +1,209 @@
+'''
+A framework for efficient model image rasterization in astropy. 
+
+This is a GSoC 2015 project and is therefore on a relatively short time frame. 
+We focus first on developing a robust solution to a simple but common use case that 
+is extensible to future implementation of more complex models and solutions. 
+
+Use case: Generate a sky model image of discrete sources, e.g. something like 
+photutils.datasets.make_100gaussians_image(), but faster and convolving a source model
+with an instrument PSF. 
+
+The primary challenge- and motivation- is maximizing speed while maintainging precision. A 
+straightforward approach is to use oversampling within a bounded area around the 
+source containing the majority of its flux (e.g. 5-sigma) for discretization and integration. 
+
+The key steps are:
+1. Model generation
+	- Convolution of a source model with an instrument PSF model
+	- Fitting the convolved model to image data if required
+2. Model discretization/Image rasterization
+	- Integration of the model over pixel areas
+	- Reprojection if required
+
+The pieces necessary are largely in existence. Relevant packages 
+include astropy.modeling, astropy.convolution, astropy.nddata, sherpa, galfit, reproject,
+... others.
+
+Outlined here is a rough draft of the proposed API for community discussion. 
+
+'''
+
+######################
+# ConvolvedModel API #
+######################
+
+from astropy.modeling import ConvolvedModel
+from astropy.modeling.models import Box2D, Lorentz1D, Gaussian2D
+from astropy.convolution import Gaussian2DKernel, Gaussian1DKernel
+
+# Initiate a ConvolvedModel instance, takes source model and convolution kernel
+# like:
+source = Box2D(1, 0, 0, 1, 1)
+psf = Gaussian2DKernel(1)
+convolved_model = ConvolvedModel(source, psf, mode='oversample', factor=10)
+
+# mode will be used both for convolution and later integration
+# if a different mode is desired for convolution it should be passed in the psf
+
+# support of 1D models
+source = Lorentz1D(1, 0, 0, 1, 1)
+psf = Gaussian1DKernel(1)
+convolved_model = ConvolvedModel(source, psf)
+
+# should there be a Convolved1DModel and Convolved2DModel?
+# should this be an extension of the CompoundModel class?
+
+# later we can think about using any model as convolution kernels
+# and even implement the analytical solution for the common Gaussian case
+# or other cases where analytical solutions are known
+source = Gaussian2D(1, 0, 0, 1, 1)
+psf = Gaussian2D(1, 0, 0, 1, 1)
+convolved_model = ConvolvedModel(source, psf)
+
+
+# the convolution method can be set via string, to choose between the astropy.convolution
+# methods
+convolved_model = ConvolvedModel(source, psf, method='fft')
+convolved_model = ConvolvedModel(source, psf, method='standard')
+
+# or by passing a function that should be used. This is mainly nececessary, because
+# the performance of the astropy convolution algorithms is rather weak...
+convolved_model = ConvolvedModel(source, psf, method=scipy.signal.fftconvolve)
+
+# recommendation would be to use scipy.signal.fftconvolve which is very fast
+
+# additional keyword arguments can be passed to the convolution function
+convolved_model = ConvolvedModel(source, psf, method='standard', mode='oversample')
+
+# the convolved_model can be fit 
+
+# Further notes:
+# 1. Implementation of this class should be straight forward and not much effort
+# 2. Testing could be set up using Gaussian models, where the analytical solution is
+# known
+# 3. Later one could think about assigning the corresponding fourier transforms to analytical
+# models and use this for the convolution, which will probably have a better performance.
+# Source model and PSF are evaluated analytically in Fourier space, mutliplied and transformed
+# back to normal space.
+# 4. As convolution is a linear operation the derivative of a convolution behaves like 
+# http://en.wikipedia.org/wiki/Convolution#Differentiation this relation can be used
+# to implement `fit_deriv` methods for convolved models.
+
+
+#########################
+# Model integration API #
+#########################
+
+import numpy as np
+
+# to avoid bias, when the scale of the model is similar to the size of the pixels,
+# models have to be integrated over pixels, for this purpose models could define
+# an integrate function, that takes upper and lower bounds of the integration
+# (see also https://github.com/astropy/astropy/issues/1586)
+model = Gaussian1D(1, 0, 1)
+model.integrate(0, np.inf)
+
+# it should also work with arrays of bounds of course
+x_lo = np.arange(-10, 11) - 0.5
+x_hi = np.arange(-10, 11) + 0.5
+model.integrate((x_lo, x_hi))
+
+# for 2d models it would be
+y_lo = np.arange(-10, 11) - 0.5
+y_hi = np.arange(-10, 11) + 0.5
+model.integrate((x_lo, x_hi, y_lo, y_hi))
+
+# whenever possible (e.g. polynomials) the integrate function would be implemented
+# analytically, for all others a standard numerical integration routine would be used
+# e.g. from scipy.integrate
+# the integration mode can be passed to the integrate method if it is not already an 
+# attribute of model
+model.integrate((x_lo, x_hi, y_lo, y_hi), mode='integrate')
+
+# for convenience one could incorporate the functionality into the standard model API
+# as is done for the convolution kernals 
+model = Gaussian1D(1, 0, 1, mode='center') 
+model = Gaussian1D(1, 0, 1, mode='linear_interp') 
+model = Gaussian1D(1, 0, 1, mode='integrate')
+model = Gaussian1D(1, 0, 1, mode='oversample', factor=10)
+
+# the given mode could also be used when the model is called for evaluation
+x = np.arange(-10, 11)
+y = np.arange(-10, 11)
+model(np.meshgrid(x,y))
+
+
+# models that are spatially confined to a certain region such as Gaussian2D, Disk2D,
+# Delta2D, Box2D, ... should have the possibility to apply a different normalization,
+# where the amplitude parameter corresponds to the integral of the model and not the
+# peak value at maximum (see e.g.: https://gammapy.readthedocs.org/en/latest/api/gammapy.morphology.Shell2D.html#gammapy.morphology.Shell2D)
+# This is very useful for modeling sources, because the amplitude parameter then
+# corresponds to the total flux of the source...
+model = Gaussian2D(1, 0, 0, 1, 1, normalization='peak')
+model = Gaussian2D(1, 0, 0, 1, 1, normalization='integral')
+
+# alternatively one could introduce new NormGaussian2D, NormBox2D, ... models, but that
+# would probably be a lot of code duplication?
+
+# all models that are 'normalisable' should have an `extent` attribute that specifies
+# a rectangular region where the model is nonzero or has 99.99% containment (or a similar
+# criterion). A reasonable value would be chosen as default and would be given as multiples 
+# model's width parameters. Model evaluation can later be limited to the region defined
+# by `extent` for better performance 
+Gaussian2D.extent= dict(x=5 * x_stddev, y=5 * y_stddev)
+Disk2D.extent = dict(x=R_0, y=R_0) 
+
+# models should have an attribute that defines a reasonable sample size for the model,
+# warnings could be raised, when models are undersampled.
+Gaussian2D.sample_size = width / 5.
+
+
+###############################
+# Model Image Computation API #
+###############################
+
+# functionality is added to NDData via NDDataMixin classes
+class SkyImage(NDData, NDSlicingMixin, NDWCSCutoutMixin) : pass
+
+# a sky image object is created as following
+sky_image = SkyImage(np.zeros(501, 10001), wcs=WCS())
+
+# it has a coordinates attribute, that returns a SkyCoord object
+# with every pixel containing the corresponding position
+sky_image.coordinates
+
+# one can make rectangular cutouts of the sky image using
+cutout = sky_image.extract(position, extent)
+cutout = sky_image.extract(postion, extent, copy=True)
+
+# which returns a SkyImage object, but with modified WCS transform 
+# data can be copied optionally
+
+# to render the model on a sky image a new function should be defined
+source_image = render_model_to_sky(model, wcs=WCS())
+
+# which would return again a SkyImage object centered on the nearest pixel of the 
+# source position of size defined by extent. The pixel values are computed according to
+# the model's mode. 
+
+# The source image can then be added to sky_image 
+sky_image.add(source_image)
+
+# later `.add` could use any kind of resampling/reprojection
+# for now we just assume that `sky_image` and `source_image`
+# pixels are aligned (which should be checked...)
+
+# Combining it all, 
+sky_image.add_sources(models)
+
+# models should be an iterable of model instances for each source.
+# then paralleization can speed things up.
+
+sky_image.add_sources(models, distribute=True, N_core='none')
+
+
+
+
+
+