diff --git a/.gitignore b/.gitignore
index 3ff90a6..deba437 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,8 @@ PKG-INFO
 requires.txt
 SOURCES.txt
 top_level.txt
+build/
+dist/
 
 # Python bytecode
 *.pyc
diff --git a/sphviewer/Camera.py b/sphviewer/Camera.py
index 688a16b..5299cd4 100644
--- a/sphviewer/Camera.py
+++ b/sphviewer/Camera.py
@@ -27,11 +27,11 @@
 class Camera(object):
     def __init__(self, x=None, y=None, z=None, r=None,
                  t=None, p=None, zoom=None, roll=None,
-                 xsize=None, ysize=None, extent=None):
+                 xsize=None, ysize=None, extent=None, projection=None):
         self._name = 'CAMERA'
         self.__params = {'x': x, 'y': y, 'z': z, 'r': r,
                          't': t, 'p': p, 'zoom': zoom, 'roll': roll,
-                         'xsize': xsize, 'ysize': ysize, 'extent': extent}
+                         'xsize': xsize, 'ysize': ysize, 'extent': extent, 'projection': projection}
 
     def get_params(self):
         return self.__params
@@ -120,4 +120,4 @@ def set_autocamera(self, Particles, mode='minmax'):
 
         self.__params = {'x': xmean, 'y': ymean, 'z': zmean, 'r': r,
                          't': 0, 'p': 0, 'zoom': 1, 'roll': 0,
-                         'xsize': 500, 'ysize': 500, 'extent': None}
+                         'xsize': 500, 'ysize': 500, 'extent': None, 'projection': None}
diff --git a/sphviewer/Render.py b/sphviewer/Render.py
index e84a477..f200147 100644
--- a/sphviewer/Render.py
+++ b/sphviewer/Render.py
@@ -42,7 +42,7 @@ def import_code(filename):
 class Render(object):
     def __init__(self, Scene):
         """
-        Produces a render of the Scene. It uses a kernel interpolation 
+        Produces a render of the Scene. It uses a kernel interpolation
         method. Its setting and getting methods are:
 
         getting methods:
@@ -61,7 +61,7 @@ def __init__(self, Scene):
 
         Other methods are:
 
-        - save 
+        - save
         - histogram
 
         """
@@ -97,6 +97,12 @@ def __init__(self, Scene):
     def __make_render(self, x, y, t, kview, xsize, ysize):
         from .extensions import render
 
+        if(self.Scene.Camera.get_params()['r'] == 'infinity'):
+            projection = 0
+        elif(self.Scene.Camera.get_params()['projection'] == 'fisheye'):
+            projection = 2
+        else:
+            projection = 1
         image = render.render(
             self.Scene._x,
             self.Scene._y,
@@ -104,12 +110,14 @@ def __make_render(self, x, y, t, kview, xsize, ysize):
             self.Scene._m,
             np.int32(xsize),
             np.int32(ysize),
+            np.int32(projection),
+            np.float32(self.Scene.get_extent()),
         )
         return np.reshape(image, [ysize, xsize])
 
     def get_image(self):
         """
-        - get_image(): Returns the image, with dimension [xsize,ysize], 
+        - get_image(): Returns the image, with dimension [xsize,ysize],
         where xsize and ysize are the number of pixels of the image defined by the Camera.
         """
         extent = self.Scene.get_extent()
@@ -155,8 +163,8 @@ def get_logscale(self):
 
     def histogram(self, axis=None, **kargs):
         """
-        - histogram(axis=None, **kargs): It computes and shows the histogram of the image. This is 
-        usefull for choosing a proper scale to the output, or for clipping some values. If 
+        - histogram(axis=None, **kargs): It computes and shows the histogram of the image. This is
+        usefull for choosing a proper scale to the output, or for clipping some values. If
         axis is None, it selects the current axis to plot the histogram.
 
         Keyword arguments:
@@ -273,33 +281,33 @@ def histogram(self, axis=None, **kargs):
         :class:`~matplotlib.patches.Patch` instances returned by *hist*:
 
         agg_filter: unknown
-        alpha: float or None         
-        animated: [True | False]         
-        antialiased or aa: [True | False]  or None for default         
-        axes: an :class:`~matplotlib.axes.Axes` instance         
-        clip_box: a :class:`matplotlib.transforms.Bbox` instance         
-        clip_on: [True | False]         
-        clip_path: [ (:class:`~matplotlib.path.Path`,         :class:`~matplotlib.transforms.Transform`) |         :class:`~matplotlib.patches.Patch` | None ]         
+        alpha: float or None
+        animated: [True | False]
+        antialiased or aa: [True | False]  or None for default
+        axes: an :class:`~matplotlib.axes.Axes` instance
+        clip_box: a :class:`matplotlib.transforms.Bbox` instance
+        clip_on: [True | False]
+        clip_path: [ (:class:`~matplotlib.path.Path`,         :class:`~matplotlib.transforms.Transform`) |         :class:`~matplotlib.patches.Patch` | None ]
         color: matplotlib color spec
-        contains: a callable function         
-        edgecolor or ec: mpl color spec, or None for default, or 'none' for no color         
-        facecolor or fc: mpl color spec, or None for default, or 'none' for no color         
-        figure: a :class:`matplotlib.figure.Figure` instance         
-        fill: [True | False]         
-        gid: an id string         
-        hatch: [ '/' | '\\' | '|' | '-' | '+' | 'x' | 'o' | 'O' | '.' | '*' ]         
-        label: any string         
-        linestyle or ls: ['solid' | 'dashed' | 'dashdot' | 'dotted']         
-        linewidth or lw: float or None for default         
-        lod: [True | False]         
+        contains: a callable function
+        edgecolor or ec: mpl color spec, or None for default, or 'none' for no color
+        facecolor or fc: mpl color spec, or None for default, or 'none' for no color
+        figure: a :class:`matplotlib.figure.Figure` instance
+        fill: [True | False]
+        gid: an id string
+        hatch: [ '/' | '\\' | '|' | '-' | '+' | 'x' | 'o' | 'O' | '.' | '*' ]
+        label: any string
+        linestyle or ls: ['solid' | 'dashed' | 'dashdot' | 'dotted']
+        linewidth or lw: float or None for default
+        lod: [True | False]
         path_effects: unknown
-        picker: [None|float|boolean|callable]         
-        rasterized: [True | False | None]         
+        picker: [None|float|boolean|callable]
+        rasterized: [True | False | None]
         snap: unknown
-        transform: :class:`~matplotlib.transforms.Transform` instance         
-        url: a url string         
-        visible: [True | False]         
-        zorder: any number 
+        transform: :class:`~matplotlib.transforms.Transform` instance
+        url: a url string
+        visible: [True | False]
+        zorder: any number
         """
         if axis == None:
             axis = plt.gca()
@@ -307,9 +315,9 @@ def histogram(self, axis=None, **kargs):
 
     def save(self, outputfile, **kargs):
         """
-        - Save the image in some common image formats. It uses the pyplot.save 
-        method.  
-        outputfile is a string containing a path to a filename, 
+        - Save the image in some common image formats. It uses the pyplot.save
+        method.
+        outputfile is a string containing a path to a filename,
         of a Python file-like object. If *format* is *None* and
         *fname* is a string, the output format is deduced from
         the extension of the filename.
diff --git a/sphviewer/Scene.py b/sphviewer/Scene.py
index d454000..6680d3c 100644
--- a/sphviewer/Scene.py
+++ b/sphviewer/Scene.py
@@ -44,14 +44,14 @@ def rotate(angle, axis, pos):
 class Scene(object):
     def __init__(self, Particles, Camera=None):
         """
-        Scene class takes a sphviewer.Particles class and computes the 
-        coordinates of the particles as seen from a Camera. It is to say, 
-        for a given particle, whose coordinates are x,y and z, Scene 
-        computes a new set of coordinates x' and y', which are the aparent 
-        coordinates of the particles as seen from a camera. 
-        By default, Camera=None means that Scene have to use an autocamera. It 
+        Scene class takes a sphviewer.Particles class and computes the
+        coordinates of the particles as seen from a Camera. It is to say,
+        for a given particle, whose coordinates are x,y and z, Scene
+        computes a new set of coordinates x' and y', which are the aparent
+        coordinates of the particles as seen from a camera.
+        By default, Camera=None means that Scene have to use an autocamera. It
         will try to define a Camera whose parameters have some sence according
-        to your data set. In case you want to prevent Scene from using an autocamera, 
+        to your data set. In case you want to prevent Scene from using an autocamera,
         you should provide a valid Camera Class.
         As Particles class, Scene has its own getting and setting methods.
 
@@ -66,7 +66,7 @@ def __init__(self, Particles, Camera=None):
         - get_scene
         - get_extent
 
-        other methods are: 
+        other methods are:
 
         - plot
 
@@ -116,9 +116,9 @@ def __init__(self, Particles, Camera=None):
 
     def set_autocamera(self, mode='density'):
         """
-        - set_autocamera(mode='density'): By default, Scene defines its 
-        own Camera. However, there is no a general way for doing so. Scene 
-        uses a density criterion for getting the point of view. If this is 
+        - set_autocamera(mode='density'): By default, Scene defines its
+        own Camera. However, there is no a general way for doing so. Scene
+        uses a density criterion for getting the point of view. If this is
         not a good option for your problem, you can choose among:
         |'minmax'|'density'|'median'|'mean'|. If None of the previous methods
         work well, you may define the camera params by yourself.
@@ -130,27 +130,27 @@ def set_autocamera(self, mode='density'):
 
     def get_scene(self):
         """
-        - get_scene(): It return the x and y position, the smoothing length 
-        of the particles and the index of the particles that are active in 
-        the scene. In principle this is an internal function and you don't 
-        need this data. 
+        - get_scene(): It return the x and y position, the smoothing length
+        of the particles and the index of the particles that are active in
+        the scene. In principle this is an internal function and you don't
+        need this data.
         """
         return self._x, self._y, self._hsml, self._m, self._kview
 
     def get_extent(self):
         """
-        - get_extent(): It returns the extent array needed for converting 
-        the image coordinates (given in pixels units) into physical coordinates. 
-        It is an array like the following one: [xmin,xmax,ymin,ymax]; it is to say, 
-        an array that contains the extreme values of the scene. 
+        - get_extent(): It returns the extent array needed for converting
+        the image coordinates (given in pixels units) into physical coordinates.
+        It is an array like the following one: [xmin,xmax,ymin,ymax]; it is to say,
+        an array that contains the extreme values of the scene.
         """
         return self.__extent
 
     def update_camera(self, **kargs):
         """
-        - update_camera(**kwarg): By using this method you can define all 
-        the new paramenters of the camera. Read the available **kwarg in 
-        the sphviewer.Camera documentation. 
+        - update_camera(**kwarg): By using this method you can define all
+        the new paramenters of the camera. Read the available **kwarg in
+        the sphviewer.Camera documentation.
         """
         self.Camera.set_params(**kargs)
         self._x, self._y, self._hsml, self._kview = self.__compute_scene()
@@ -189,9 +189,46 @@ def __compute_scene(self):
             else:
                 self.__extent = np.float32(self._camera_params['extent'])
         else:
-            projection = 1
             rcam = np.float32(self._camera_params['r'])
             self.__extent = np.array([0, 0, 0, 0], dtype=np.float32)
+            if(isinstance(self._camera_params['projection'],str)):
+                if('fisheye' in self._camera_params['projection']):
+                    projection = 2
+                    try:
+                        shift = float(self._camera_params['projection'][7:])
+                        theta_rad=np.pi*np.array(theta)/180.
+                        phi_rad=np.pi*np.array(phi)/180.
+                        roll_rad=np.pi*np.array(roll)/180.
+                        kaxis = np.array([np.cos(phi_rad)*np.cos(roll_rad),
+                                          np.cos(theta_rad)*np.sin(roll_rad)+np.sin(theta_rad)*np.sin(phi_rad)*np.cos(roll_rad),
+                                          np.sin(theta_rad)*np.sin(roll_rad)-np.cos(theta_rad)*np.sin(phi_rad)*np.cos(roll_rad)]).T
+                        Kchange = np.cross(kaxis,np.eye(3))
+                        Rchange = np.eye(3)-np.sin(shift*np.pi/180)*Kchange+(1-np.cos(shift*np.pi/180))*np.dot(Kchange,Kchange)
+                        zaxis = np.array([np.sin(phi_rad),-np.sin(theta_rad)*np.cos(phi_rad),np.cos(theta_rad)*np.cos(phi_rad)]).T
+                        jaxis = np.matmul(Rchange,zaxis)
+                        xcopy,ycopy,zcopy = np.array([xcam,ycam,zcam])-rcam*(zaxis-jaxis)
+                        iaxis = np.cross(jaxis,kaxis)
+                        tcopy_rad = np.arctan2(-jaxis[1],jaxis[2]); tcopy = tcopy_rad/np.pi*180
+                        cossintcopy=np.array([np.cos(tcopy_rad),np.sin(tcopy_rad)])
+                        rollcopy_rad = np.arctan2(-iaxis[0],kaxis[0]); rollcopy = rollcopy_rad/np.pi*180
+                        cossinrollcopy=np.array([np.cos(rollcopy_rad),np.sin(rollcopy_rad)])
+                        pcopy_rad=np.arctan2(jaxis[0],np.nanmean([([kaxis[0],-iaxis[0]]/cossinrollcopy)[np.argmax(np.abs(cossintcopy))],
+                                                                  ([jaxis[2],-jaxis[1]]/cossintcopy)[np.argmax(np.abs(cossintcopy))]]))
+                        pcopy=pcopy_rad/np.pi*180
+                        xcam,ycam,zcam,theta,phi,roll = xcopy,ycopy,zcopy,tcopy,pcopy,rollcopy
+                    except:
+                        None
+                else:
+                    raise ValueError("'"+self._camera_params['projection']+"' isn't a valid projection type")
+            else:
+                projection = 1
+
+        xx,yy,hh,kk = scene.scene(pos[:,0],pos[:,1],
+                                  pos[:,2],hsml,
+                                  xcam, ycam, zcam, rcam, theta, phi, roll,
+                                  zoom, self.__extent, xsize, ysize, projection)
+
+        return xx,yy,hh,kk
 
         xx, yy, hh, kk = scene.scene(pos[:, 0], pos[:, 1],
                                      pos[:, 2], hsml,
@@ -298,7 +335,7 @@ def __compute_scene_old(self):
 
     def plot(self, axis=None, **kargs):
         """
-        - plot(axis=None, **kwarg): Finally, sphviewer.Scene class has its own plotting method. 
+        - plot(axis=None, **kwarg): Finally, sphviewer.Scene class has its own plotting method.
         It shows the scene as seen by the camera. It is to say, it plots the particles according
         to their aparent coordinates; axis makes a reference to an existing axis. In case axis is None,
         the plot is made on the current axis.
@@ -306,46 +343,46 @@ def plot(self, axis=None, **kargs):
         The kwargs are :class:`~matplotlib.lines.Line2D` properties:
 
         agg_filter: unknown
-        alpha: float (0.0 transparent through 1.0 opaque)         
-        animated: [True | False]         
-        antialiased or aa: [True | False]         
-        axes: an :class:`~matplotlib.axes.Axes` instance         
-        clip_box: a :class:`matplotlib.transforms.Bbox` instance         
-        clip_on: [True | False]         
-        clip_path: [ (:class:`~matplotlib.path.Path`,         :class:`~matplotlib.transforms.Transform`) |         :class:`~matplotlib.patches.Patch` | None ]         
-        color or c: any matplotlib color         
-        contains: a callable function         
-        dash_capstyle: ['butt' | 'round' | 'projecting']         
-        dash_joinstyle: ['miter' | 'round' | 'bevel']         
-        dashes: sequence of on/off ink in points         
-        data: 2D array (rows are x, y) or two 1D arrays         
-        drawstyle: [ 'default' | 'steps' | 'steps-pre' | 'steps-mid' | 'steps-post' ]         
-        figure: a :class:`matplotlib.figure.Figure` instance         
-        fillstyle: ['full' | 'left' | 'right' | 'bottom' | 'top']         
-        gid: an id string         
-        label: any string         
-        linestyle or ls: [ ``'-'`` | ``'--'`` | ``'-.'`` | ``':'`` | ``'None'`` | ``' '`` | ``''`` ]         and any drawstyle in combination with a linestyle, e.g. ``'steps--'``.         
-        linewidth or lw: float value in points         
-        lod: [True | False]         
+        alpha: float (0.0 transparent through 1.0 opaque)
+        animated: [True | False]
+        antialiased or aa: [True | False]
+        axes: an :class:`~matplotlib.axes.Axes` instance
+        clip_box: a :class:`matplotlib.transforms.Bbox` instance
+        clip_on: [True | False]
+        clip_path: [ (:class:`~matplotlib.path.Path`,         :class:`~matplotlib.transforms.Transform`) |         :class:`~matplotlib.patches.Patch` | None ]
+        color or c: any matplotlib color
+        contains: a callable function
+        dash_capstyle: ['butt' | 'round' | 'projecting']
+        dash_joinstyle: ['miter' | 'round' | 'bevel']
+        dashes: sequence of on/off ink in points
+        data: 2D array (rows are x, y) or two 1D arrays
+        drawstyle: [ 'default' | 'steps' | 'steps-pre' | 'steps-mid' | 'steps-post' ]
+        figure: a :class:`matplotlib.figure.Figure` instance
+        fillstyle: ['full' | 'left' | 'right' | 'bottom' | 'top']
+        gid: an id string
+        label: any string
+        linestyle or ls: [ ``'-'`` | ``'--'`` | ``'-.'`` | ``':'`` | ``'None'`` | ``' '`` | ``''`` ]         and any drawstyle in combination with a linestyle, e.g. ``'steps--'``.
+        linewidth or lw: float value in points
+        lod: [True | False]
         marker: [ ``7`` | ``4`` | ``5`` | ``6`` | ``'o'`` | ``'D'`` | ``'h'`` | ``'H'`` | ``'_'`` | ``''`` | ``'None'`` | ``' '`` | ``None`` | ``'8'`` | ``'p'`` | ``','`` | ``'+'`` | ``'.'`` | ``'s'`` | ``'*'`` | ``'d'`` | ``3`` | ``0`` | ``1`` | ``2`` | ``'1'`` | ``'3'`` | ``'4'`` | ``'2'`` | ``'v'`` | ``'<'`` | ``'>'`` | ``'^'`` | ``'|'`` | ``'x'`` | ``'$...$'`` | *tuple* | *Nx2 array* ]
-        markeredgecolor or mec: any matplotlib color         
-        markeredgewidth or mew: float value in points         
-        markerfacecolor or mfc: any matplotlib color         
-        markerfacecoloralt or mfcalt: any matplotlib color         
-        markersize or ms: float         
+        markeredgecolor or mec: any matplotlib color
+        markeredgewidth or mew: float value in points
+        markerfacecolor or mfc: any matplotlib color
+        markerfacecoloralt or mfcalt: any matplotlib color
+        markersize or ms: float
         markevery: None | integer | (startind, stride)
-        picker: float distance in points or callable pick function         ``fn(artist, event)``         
-        pickradius: float distance in points         
-        rasterized: [True | False | None]         
+        picker: float distance in points or callable pick function         ``fn(artist, event)``
+        pickradius: float distance in points
+        rasterized: [True | False | None]
         snap: unknown
-        solid_capstyle: ['butt' | 'round' |  'projecting']         
-        solid_joinstyle: ['miter' | 'round' | 'bevel']         
-        transform: a :class:`matplotlib.transforms.Transform` instance         
-        url: a url string         
-        visible: [True | False]         
-        xdata: 1D array         
-        ydata: 1D array         
-        zorder: any number         
+        solid_capstyle: ['butt' | 'round' |  'projecting']
+        solid_joinstyle: ['miter' | 'round' | 'bevel']
+        transform: a :class:`matplotlib.transforms.Transform` instance
+        url: a url string
+        visible: [True | False]
+        xdata: 1D array
+        ydata: 1D array
+        zorder: any number
 
         kwargs *scalex* and *scaley*, if defined, are passed on to
         :meth:`~matplotlib.axes.Axes.autoscale_view` to determine
diff --git a/sphviewer/extensions/rendermodule.c b/sphviewer/extensions/rendermodule.c
index ff8c47f..093126c 100644
--- a/sphviewer/extensions/rendermodule.c
+++ b/sphviewer/extensions/rendermodule.c
@@ -1,6 +1,6 @@
 /*////////////////////////////////////////////////////////////////////
     This file is part of Py-SPHViewer
-    
+
     <Py-SPHVIewer is a framework for rendering particles in Python
     using the SPH interpolation scheme.>
     Copyright (C) <2013>  <Alejandro Benitez Llambay>
@@ -33,7 +33,7 @@
 
 float *get_double_array(PyArrayObject *array_obj, int n){
   /* This function returns the data stored in a double PyArrayObject*/
-  double *local_array = (double *)array_obj->data;  
+  double *local_array = (double *)array_obj->data;
   float *output = (float *)malloc( n * sizeof(float) );
 
 #pragma omp parallel for firstprivate(n)
@@ -49,35 +49,45 @@ float cubic_kernel(float r, float h){
   //2D Dome-shapeed quadratic Kernel (1-R^2) (Hicks and Liebrock 2000).
   float func;
   float sigma;
-  
+
   sigma = 15.0/(8.0*3.141592);
-  if(r/h <= 1.0)  
+  if(r/h <= 1.0)
     func = 4.0/3.0*h*pow(sqrt(1.0-(r/h)*(r/h)),3);
-  if(r/h > 1.0)  
+  if(r/h > 1.0)
     func = 0;
   return sigma*func/(h*h*h);
 }
 
+float sinc(float x){
+  float value;
+
+  if(x==0)
+    value = 1.0;
+  else
+    value = sin(x)/x;
+  return value;
+}
+
+
+void c_render(float *x, float *y, float *t, float *mass,
+	      int xsize, int ysize, int n, int projection, float *extent, float *image){
 
-void c_render(float *x, float *y, float *t, float *mass, 
-	      int xsize, int ysize, int n, float *image){ 
-  
   // C function calculating the image of the particles convolved with our kernel
   int size_lim;
- 
+
   if(xsize >= ysize){
     size_lim = xsize;
   }
   else{
     size_lim = ysize;
   }
-  
+
 #pragma omp parallel
   {
     float *local_image;
     int i,j,k,l;
     int xx, yy, tt;
-    float mm;
+    float tt_f, mm;
     int r, nth, ppt, thread_id;
 
   nth = omp_get_num_threads();                // Get number of threads
@@ -85,8 +95,8 @@ void c_render(float *x, float *y, float *t, float *mass,
   ppt = n/nth;                                //  Number or Particles per thread (ppt)
   r = n-ppt*nth;                              // Remainder
 
-  
-  local_image = (float *)malloc( xsize * ysize * sizeof(float) ); 
+
+  local_image = (float *)malloc( xsize * ysize * sizeof(float) );
 
   // Let's initialize the image to zero
 
@@ -95,60 +105,88 @@ void c_render(float *x, float *y, float *t, float *mass,
       local_image[k*xsize+j] = 0.0;
     }
   }
-  
-
-  // Let's compute the local image 
-  //  for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){
-  for(l=0;l<ppt;l++){
-    i = thread_id+nth*l;
-    xx = (int) x[i];
-    yy = (int) y[i];
-    tt = (int) t[i];
-    mm = mass[i];
-
-    if(tt <= 1) {
-      local_image[yy*xsize+xx] += mm;      
-      continue;
-    }
-    
-    if(tt > size_lim) tt = size_lim;
-
-    
-    // Let's compute the convolution with the Kernel
-    for(j=-tt; j<tt+1; j++){
-      for(k=-tt; k<tt+1; k++){
-	if( ( (xx+j) >= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){
-	  local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt); 
-	}
+  // Dome projection/Fish-eye camera : Azimuth Equidistant projection distortions
+  if(projection==2){
+    float xxrad, yyrad, ttrad, rho, sincrho, cosrho;
+    if((r-thread_id) > 0) ppt+=1; //to include the remainder particle
+    // Let's compute the local image
+    //  for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){
+    for(l=0;l<ppt;l++){
+      i = thread_id+nth*l;
+      xx = (int) x[i];
+      yy = (int) y[i];
+      xxrad = (((float)xx/(float)xsize)*(extent[1]-extent[0])+extent[0])*3.141592;
+      yyrad = (((float)yy/(float)ysize)*(extent[3]-extent[2])+extent[2])*3.141592;
+      rho = sqrt(xxrad*xxrad+yyrad*yyrad);
+      sincrho = sinc(rho);
+      cosrho = cos(rho);
+      tt_f = t[i];
+      ttrad = ((tt_f/(float)xsize)*(extent[1]-extent[0]))*3.141592;
+      tt = (int)ceil(tt_f/sincrho);
+      mm = mass[i];
+
+      if(tt <= 1) {
+        if( (rho <= 3.141592) && ( xx >= 0) && ( xx < xsize) && ( yy >= 0) && ( yy < ysize)){
+          local_image[yy*xsize+xx] += mm;
+        }
+      }
+      else{
+        if(tt > size_lim) tt = size_lim;
+
+        int jxx,kyy;
+        float jxxrad, kyyrad, rhopix, cosd;
+        // Let's compute the convolution with the Kernel
+        for(j=-tt; j<tt+1; j++){
+          for(k=-tt; k<tt+1; k++){
+            jxx = j+xx;
+            kyy = k+yy;
+            jxxrad = ((((float)jxx)/(float)xsize)*(extent[1]-extent[0])+extent[0])*3.141592;
+            kyyrad = ((((float)kyy)/(float)ysize)*(extent[3]-extent[2])+extent[2])*3.141592;
+            rhopix = sqrt(jxxrad*jxxrad+kyyrad*kyyrad);
+            if( (rhopix <= 3.141592) && (jxx>=0) && (jxx<xsize) && (kyy>=0) && (kyy<ysize) ){
+              cosd = (xxrad*jxxrad+yyrad*kyyrad)*sincrho*sinc(rhopix)+cosrho*cos(rhopix);
+              if(cosd>=1) cosd = 1; else if(cosd<=-1) cosd = -1;
+              local_image[kyy*xsize+jxx] += mm*cubic_kernel(acos(cosd), ttrad); //I should normalize here by the surface area
+            }
+          }
+        }
       }
     }
-    
   }
-  
-  // Let's compute the image for the remainder particles...
-  if((r-thread_id) > 0){
-    i  = nth*ppt+thread_id;
-    xx = (int) x[i];
-    yy = (int) y[i];
-    tt = (int) t[i];
-    mm = mass[i];
-    
-    if(tt <= 1){
-      local_image[yy*xsize+xx] += mm;
-    }
-  
-    if(tt > size_lim) tt = size_lim;
-    
-    for(j=-tt; j<tt+1; j++){
-      for(k=-tt; k<tt+1; k++){
-	if( ( (xx+j) >= 0) && ( (xx+j) < xsize) && ( (yy+k) >=0) && ( (yy+k) < ysize)){
-	  local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt);
-	}
+  //
+  else{
+    if((r-thread_id) > 0) ppt+=1; //to include the remainder particle
+    // Let's compute the local image
+    //  for(i=(thread_id*ppt); i<(thread_id+1)*ppt; i++){
+    for(l=0;l<ppt;l++){
+      i = thread_id+nth*l;
+      xx = (int) x[i];
+      yy = (int) y[i];
+      tt_f = t[i];
+      tt = (int)ceil(tt_f);
+      mm = mass[i];
+
+      if(tt <= 1) {
+        if( ( xx >= 0) && ( xx < xsize) && ( yy >= 0) && ( yy < ysize)){
+          local_image[yy*xsize+xx] += mm;
+        }
+      }
+      else{
+        if(tt > size_lim) tt = size_lim;
+
+        // Let's compute the convolution with the Kernel
+        for(j=-tt; j<tt+1; j++){
+          for(k=-tt; k<tt+1; k++){
+  	        if( ( (xx+j) >= 0) && ( (xx+j) < xsize) && ( (yy+k) >= 0) && ( (yy+k) < ysize)){
+              local_image[(yy+k)*xsize+(xx+j)] += mm*cubic_kernel(sqrt((float)j*(float)j+(float)k*(float)k), tt_f);
+            }
+          }
+        }
       }
     }
   }
   // Let's merge the local images
-  
+
 
   #pragma omp critical
   {
@@ -166,24 +204,25 @@ void c_render(float *x, float *y, float *t, float *mass,
 
 void test_C(){
   // This function if for testing purposes only. It writes a file called image_test.bin
-  int *x, *y, *t;
-  float *mass;
-  int xsize, ysize, n;
-  float *image;
+  int *x, *y;
+  float *t, *mass;
+  int xsize, ysize, n, projection;
+  float *extent, *image;
   int i;
 
   xsize = 1000;
   ysize = 1000;
   n = 10000;
 
-  x = (int *)malloc( n * sizeof(int) ); 
-  y = (int *)malloc( n * sizeof(int) ); 
-  t = (int *)malloc( n * sizeof(int) ); 
-  mass = (float *)malloc( n * sizeof(float) ); 
-  image = (float *)malloc( xsize * ysize * sizeof(float) ); 
+  x = (int *)malloc( n * sizeof(int) );
+  y = (int *)malloc( n * sizeof(int) );
+  t = (float *)malloc( n * sizeof(int) );
+  mass = (float *)malloc( n * sizeof(float) );
+  image = (float *)malloc( 4 * sizeof(float) );
+  image = (float *)malloc( xsize * ysize * sizeof(float) );
 
   srand( time(NULL) );
-  
+
   for(i=0;i<n;i++){
     x[i] = rand() % xsize;
     y[i] = rand() % ysize;
@@ -191,12 +230,12 @@ void test_C(){
     mass[i] = rand() % 499;
   }
 
-  c_render(x,y,t,mass,xsize,ysize,n,image);
+  c_render(x,y,t,mass,xsize,ysize,n,projection,extent,image);
 
   FILE *output;
 
   output = fopen("image_test.bin","wb");
-  fwrite(image, sizeof(float), xsize*ysize, output);   
+  fwrite(image, sizeof(float), xsize*ysize, output);
   fclose(output);
 }
 
@@ -204,23 +243,23 @@ void test_C(){
 
 static PyObject *rendermodule(PyObject *self, PyObject *args){
   PyArrayObject *x_obj, *y_obj, *t_obj;
-  PyArrayObject *m_obj;
+  PyArrayObject *m_obj, *extent_obj;
   float *x, *y, *t;
   float *mass;
   int xsize, ysize;
-  int n;
-  float *image;
+  int n, projection;
+  float *extent, *image;
   int DOUBLE = 0;
 
-  if(!PyArg_ParseTuple(args, "OOOOii",&x_obj, &y_obj, &t_obj, &m_obj, &xsize, &ysize))
+  if(!PyArg_ParseTuple(args, "OOOOiiiO",&x_obj, &y_obj, &t_obj, &m_obj, &xsize, &ysize, &projection, &extent_obj))
     return NULL;
-    
+
   // Let's check the size of the 1-dimensions arrays.
   n = (int) m_obj->dimensions[0];
 
   // Let's point the C arrays to the numpy arrays
   x    = (float *)x_obj->data;
-  y    = (float *)y_obj->data; 
+  y    = (float *)y_obj->data;
   t    = (float *)t_obj->data; /* These are always floats, as they come from Scene */
 
 
@@ -237,6 +276,8 @@ static PyObject *rendermodule(PyObject *self, PyObject *args){
     return NULL;
   }
 
+  extent = (float *)extent_obj->data;
+
   image = (float *)malloc( xsize * ysize * sizeof(float) );
 
   int i, j;
@@ -247,17 +288,17 @@ static PyObject *rendermodule(PyObject *self, PyObject *args){
   }
 
   // Here we do the work
-  c_render(x,y,t,mass,xsize,ysize,n,image);
+  c_render(x,y,t,mass,xsize,ysize,n,projection,extent,image);
 
   if(DOUBLE) free(mass);
-  
+
   // Let's build a numpy array
   npy_intp dims[1] = {xsize*ysize};
   PyArrayObject *image_obj = (PyArrayObject *) PyArray_SimpleNewFromData(1, dims, NPY_FLOAT32, image);
   image_obj->flags = NPY_OWNDATA;
-  
+
   return Py_BuildValue("N", image_obj);
-		       
+
 }
 
 static PyMethodDef RenderMethods[] = {
@@ -298,6 +339,3 @@ PyMODINIT_FUNC initrender(void) {
 //{
 //  test_C();
 //}
-
-
-
diff --git a/sphviewer/extensions/scenemodule.c b/sphviewer/extensions/scenemodule.c
index ae6e488..339dfa4 100644
--- a/sphviewer/extensions/scenemodule.c
+++ b/sphviewer/extensions/scenemodule.c
@@ -142,16 +142,6 @@ long int compute_scene(float *x, float *y, float *z, float *hsml,
     ymin = extent[2];
     ymax = extent[3];
 
-
-    for(i=0;i<n;i++){
-      if( (x[i] >= xmin) & (x[i] <= xmax) & 
-	  (y[i] >= ymin) & (y[i] <= ymax) ) {
-	
-	kview[idx] = i;
-	idx += 1;
-      }
-    }
-
 #pragma omp parallel for firstprivate(n,xmin,xmax,ymin,ymax,xsize,ysize)
     for(i=0;i<n;i++){
       x[i] = (x[i] - xmin) / (xmax-xmin) * xsize;
@@ -161,7 +151,7 @@ long int compute_scene(float *x, float *y, float *z, float *hsml,
   }
   // If the camera is not at the infinity, let's put it at a certain 
   // distance from the object.
-  else
+  else if(projection == 1)
     {
       float FOV = 2.0*fabsf(atanf(1.0/zoom));
       xmax = 1.0;
@@ -181,17 +171,6 @@ long int compute_scene(float *x, float *y, float *z, float *hsml,
       extent[3] = yfovmax;
 
       float zpart;
-      for(i=0;i<n;i++){
-	zpart = (z[i]-(-1.0*r));
-	if( (zpart > 0) & 
-	    (fabsf(x[i]) <= fabsf(zpart)*xmax/zoom) &
-	    (fabsf(y[i]) <= fabsf(zpart)*ymax/zoom) ) {
-	  
-	  kview[idx] = i;
-	  idx += 1;
-	}
-      }
-
 #pragma omp parallel for firstprivate(n,xmin,xmax,ymin,ymax,xsize,ysize,lbin,zoom,r,zpart)
       for(i=0;i<n;i++){
 	zpart = (z[i]-(-1.0*r))/zoom;
@@ -200,6 +179,41 @@ long int compute_scene(float *x, float *y, float *z, float *hsml,
 	hsml[i] = (hsml[i]/zpart)/(xmax-xmin) * xsize;
       }
     }          
+  // Dome projection/Fish-eye camera
+  else if(projection == 2)
+    {
+      float xfovmax = 2.0*fabsf(atanf(1.0/zoom))/3.141592;
+      float xfovmin = -xfovmax;
+      float yfovmax = 0.5*(xfovmax-xfovmin)*ysize/xsize;
+      float yfovmin = -yfovmax;
+
+      extent[0] = xfovmin;
+      extent[1] = xfovmax;
+      extent[2] = yfovmin;
+      extent[3] = yfovmax;
+
+      
+      float zpart, radius, polar;
+#pragma omp parallel for firstprivate(n,xfovmin,xfovmax,yfovmin,yfovmax,xsize,ysize,lbin,zoom,r,zpart,radius,polar)
+      for(i=0;i<n;i++){
+	zpart = (z[i]-(-1.0*r));
+	radius = sqrtf(x[i]*x[i]+y[i]*y[i]+zpart*zpart);
+	polar = acosf(zpart/radius)/3.141592/sqrtf(x[i]*x[i]+y[i]*y[i]);
+	x[i] = (polar*x[i]-xfovmin)/(xfovmax-xfovmin)*xsize;
+	y[i] = (polar*y[i]-yfovmin)/(yfovmax-yfovmin)*ysize;
+	hsml[i] = (hsml[i]/radius)/3.141592/(xfovmax-xfovmin)*xsize;
+      }
+    }          
+
+      for(i=0;i<n;i++){
+	if( (hsml[i] > 0) &
+	    (x[i] >= -hsml[i]) & (x[i] <= xsize+hsml[i]) &
+	    (y[i] >= -hsml[i]) & (y[i] <= ysize+hsml[i]) ) { //this considers edge-particles with intersecting kernels
+	  
+	  kview[idx] = i;
+	  idx += 1;
+      }
+    }
   return idx;
 }