diff --git a/setup.py b/setup.py
index ed83d2f38cc344e426db31e507a05ac55a7f98e9..5abd22a545bbdb88c22af962e931255b4d066045 100644
--- a/setup.py
+++ b/setup.py
@@ -1,9 +1,10 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
-from setuptools import setup, find_packages
+from setuptools import find_packages
+from setuptools import setup
 
-description = "Create figures automatically from LPy"
+description = "Create figures automatically from L-Py."
 readme = open('README.md').read()
 
 # find packages
@@ -19,9 +20,7 @@ setup_kwds = dict(
     url='',
     license='LGPL-3.0',
     zip_safe=False,
-
     packages=pkgs,
-
     package_dir={'': 'src'},
     entry_points={},
     keywords='',
diff --git a/src/lpy_tools/__init__.py b/src/lpy_tools/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..56fafa58b3f43decb7699b93048b8b87e0f695aa 100644
--- a/src/lpy_tools/__init__.py
+++ b/src/lpy_tools/__init__.py
@@ -0,0 +1,2 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
diff --git a/src/lpy_tools/point_sampler/__init__.py b/src/lpy_tools/point_sampler/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..56fafa58b3f43decb7699b93048b8b87e0f695aa 100644
--- a/src/lpy_tools/point_sampler/__init__.py
+++ b/src/lpy_tools/point_sampler/__init__.py
@@ -0,0 +1,2 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
diff --git a/src/lpy_tools/point_sampler/generatePointCloud.py b/src/lpy_tools/point_sampler/generatePointCloud.py
index 8ec0fe1ddf2d2209a810efc18ea7ee24a754ae91..0f2826d312b8b0b75b50bf9266d547fbc0bc31e3 100755
--- a/src/lpy_tools/point_sampler/generatePointCloud.py
+++ b/src/lpy_tools/point_sampler/generatePointCloud.py
@@ -1,21 +1,28 @@
-### Main automating script for point sampling (generating labelled point cloud on virtual model) ###
-### Usage: python generatePointCloud.py LpyModelFile labelDictionaryFile numberOfPoints ###
-### Example: python generatePointCloud.py Arabidopsis.lpy labelDictionary.txt 1000 ###
-### The label dictionary file "labelDictionary.txt" should be updated accordingly for each model ###
-### Author: Ayan Chaudhury (ayanchaudhury.cs@gmail.com) ###
-### INRIA team MOSAIC ###
-### Updated: August 2021 ###
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Main automating script for point sampling (generating labelled point cloud on virtual model)
+Usage: python generatePointCloud.py LpyModelFile labelDictionaryFile numberOfPoints
+Example: python generatePointCloud.py Arabidopsis.lpy labelDictionary.txt 1000
+The label dictionary file "labelDictionary.txt" should be updated accordingly for each model
+
+Author: Ayan Chaudhury (ayanchaudhury.cs@gmail.com)
+INRIA team MOSAIC
+Updated: August 2021
+"""
 
-from openalea import lpy
-from openalea.plantgl.all import *
 import argparse
+
+from openalea import lpy
+
 from scan_utils import *
 
 coloredPointCloudFileName = "pointsWithColor.xyz"
 rawPointCloudFileName = "rawPoints.xyz"
 labelFileName = "rawLabels.txt"
-#dictFileName = "labelDictionary.txt"
-#model_filename = "Arabido.lpy"
+# dictFileName = "labelDictionary.txt"
+# model_filename = "Arabido.lpy"
 
 
 parser = argparse.ArgumentParser()
@@ -27,10 +34,8 @@ model_filename = args.InputFileName
 dictFileName = args.InputLabelDictionary
 pointsToBeSampled = int(args.TotalNumberOfPoints)
 
-
-
 lsys = lpy.Lsystem(model_filename)
 lstring = lsys.derive()
 lscene = lsys.sceneInterpretation(lstring)
-pointSampler(lstring, lscene, dictFileName, pointsToBeSampled, coloredPointCloudFileName, rawPointCloudFileName, labelFileName)
-
+pointSampler(lstring, lscene, dictFileName, pointsToBeSampled, coloredPointCloudFileName, rawPointCloudFileName,
+             labelFileName)
diff --git a/src/lpy_tools/point_sampler/scan_utils.py b/src/lpy_tools/point_sampler/scan_utils.py
index 90949dcb209afe1b7924d9f180bf4f1bba7ce20a..124549aa2722f777ee20913721df0ea8ab778132 100755
--- a/src/lpy_tools/point_sampler/scan_utils.py
+++ b/src/lpy_tools/point_sampler/scan_utils.py
@@ -1,151 +1,182 @@
-### This script contains the functions to tessalate the shape into triangles, and performs the point sampling ###
-### Author: Ayan Chaudhury (ayanchaudhury.cs@gmail.com) ###
-### INRIA team MOSAIC ###
-### Updated: August 2021 ###
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+This script contains the functions to tessellate the shape into triangles, and performs the point sampling.
+
+Author: Ayan Chaudhury (ayanchaudhury.cs@gmail.com)
+INRIA team MOSAIC
+Updated: August 2021
+"""
 
-from openalea.plantgl.all import *
-import numpy as np
-from math import *
-from utilities import *
 import ast
 
+from openalea.plantgl.algo import Tesselator
+from openalea.plantgl.math import Vector3
+from openalea.plantgl.math import cross
+from openalea.plantgl.scenegraph import Cylinder
+from openalea.plantgl.scenegraph import Frustum
+from openalea.plantgl.scenegraph import Oriented
+from openalea.plantgl.scenegraph import Primitive
+from openalea.plantgl.scenegraph import Sphere
+from openalea.plantgl.scenegraph import Translated
+
+from lpy_tools.point_sampler.utilities import addTriangleIdToGlobalList
+from lpy_tools.point_sampler.utilities import addVertexToGlobalList
+from lpy_tools.point_sampler.utilities import createLabelledPointForDisplay
+from lpy_tools.point_sampler.utilities import determineInsideness
+from lpy_tools.point_sampler.utilities import extractCoord
+from lpy_tools.point_sampler.utilities import obtainTriangleIdxList
+from lpy_tools.point_sampler.utilities import obtainVertexList
+from lpy_tools.point_sampler.utilities import sampleRandomPoints
+from lpy_tools.point_sampler.utilities import selectRandomTriangle
+from lpy_tools.point_sampler.utilities import write2File
+
 
-# Given a geometry in the lscene, the function performs tessalation to
+# Given a geometry in the lscene, the function performs tesselation to
 # decompose the shape into triangles. Returns the list of points & the
 # corresponding triangle indices.
 def performTessalation(currentGeometry):
-  t = Tesselator()
-  currentGeometry.apply(t)
-  triangleset = t.result
-  ptList = triangleset.pointList
-  idxList = triangleset.indexList
-  return ptList, idxList
-
+    t = Tesselator()
+    currentGeometry.apply(t)
+    triangleset = t.result
+    ptList = triangleset.pointList
+    idxList = triangleset.indexList
+    return ptList, idxList
 
 
 # This is the main function that performs the point sampling. It takes as argument the lstring, lscene coming
 # from the virtual model, the label dictionary as text file total number of points to be sampled, and the 3
 # output file names 
 def pointSampler(lstring, lscene, dictionaryFile, totalNumberOfResampledPts, pointsWithColor, rawPoints, rawLabels):
-  #totalNumberOfResampledPts = 2048 # desired number of points
-  randomPointLabels = [] # stores the labels of all points
-  globalVertexList = [] #stores all vertices (without repeatition)
-  globalTriangleIdList = [] #stores id of vertices (as in globalVertexList) for the triangles
-  globalTriangleLabelList = [] #stores the label of each triangle of globalTriangleIdList
-  globalGeometryInfoList = [] #stores the geometry of all primitives
-  cyllist = [] # stores all cylinder primitives
-  geomInfoList = [] # stores geometry of all the primitives
+    # totalNumberOfResampledPts = 2048 # desired number of points
+    randomPointLabels = []  # stores the labels of all points
+    globalVertexList = []  # stores all vertices (without repetition)
+    globalTriangleIdList = []  # stores id of vertices (as in globalVertexList) for the triangles
+    globalTriangleLabelList = []  # stores the label of each triangle of globalTriangleIdList
+    globalGeometryInfoList = []  # stores the geometry of all primitives
+    cyllist = []  # stores all cylinder primitives
+    geomInfoList = []  # stores geometry of all the primitives
+
+    # reading the label dictionary from the file
+    with open(dictionaryFile) as f:
+        data = f.read()
+    # reconstructing the data as a dictionary
+    labelTable = ast.literal_eval(data)
 
-  # reading the label dictionary from the file
-  with open(dictionaryFile) as f:
-    data = f.read()
-  # reconstructing the data as a dictionary
-  labelTable = ast.literal_eval(data)
+    # loop over all elements in the lscene
+    for shape in lscene:
+        botpoint = Vector3(0, 0, 0)  # the model is based at the origin
+        heading = Vector3(0, 0, 1)
+        geometry = shape
+        id = shape.id
+        currentModule = lstring[id].name
+        labelCurrentModule = labelTable[currentModule]
+        [ptList, idxList] = performTessalation(shape.geometry)  # get the triangle corresponding to the primitive
+        vertexList = obtainVertexList(ptList)  # store the set of indices of vertices that form triangles (local index)
+        triangleList = obtainTriangleIdxList(idxList)  # store the set of triangle indices
+        ### Now perform the global storing operations ###
+        globalVertexList = addVertexToGlobalList(vertexList, globalVertexList)
+        [globalTriangleIdList, globalTriangleLabelList] = addTriangleIdToGlobalList(vertexList, globalVertexList,
+                                                                                    triangleList, globalTriangleIdList,
+                                                                                    globalTriangleLabelList,
+                                                                                    [labelCurrentModule])
+        # the following loop runs over all the basic primitives we have considered in the model
+        while not isinstance(geometry, Primitive):
+            geometry = geometry.geometry
+            ### When the primitive is a cylinder ###
+            if isinstance(geometry, Cylinder):
+                cyllist.append((botpoint, botpoint + heading * geometry.height, geometry.radius))
+                toppoint = botpoint + heading * geometry.height
+                rad = geometry.radius
+                # print("Cylinder with label:", lstring[id][2], "Bottom point coordinate:", botpoint, "Top point coordinate:", toppoint, "Radius:",rad)
+                # Now we store the label and geometry information of the primitive as the following tuple:
+                # (primitive label, bottom_x, bottom_y, bottom_z, top_x, top_y, top_z, radius)
+                primitiveLabel = [1]  # the label for cylinder is given as '1'
+                geomInfoList.append(primitiveLabel)
+                [bot_x, bot_y, bot_z] = extractCoord(botpoint)  # bottom point
+                l = len(geomInfoList)
+                geomInfoList[l - 1].append(bot_x)
+                geomInfoList[l - 1].append(bot_y)
+                geomInfoList[l - 1].append(bot_z)
+                [top_x, top_y, top_z] = extractCoord(toppoint)  # top point
+                geomInfoList[l - 1].append(top_x)
+                geomInfoList[l - 1].append(top_y)
+                geomInfoList[l - 1].append(top_z)
+                geomInfoList[l - 1].append(rad)  # radius
+            ### When the primitive is a frustum ###
+            elif isinstance(geometry, Frustum):
+                frustumTopPt = botpoint + heading * geometry.height
+                frustumRad = geometry.radius
+                taperFactor = geometry.taper
+                topRad = frustumRad * taperFactor
+                # print("Frustum with label:", lstring[id][2], "Bottom point coordinate:", botpoint, "Top point coordinate:", frustumTopPt, "Base radius:", frustumRad, "Top radius:", topRad)
+                # Now we store the label and geometry information of the primitive as the following tuple:
+                # (primitive label, bottom_x, bottom_y, bottom_z, top_x, top_y, top_z, base radius, top radius)
+                primitiveLabel = [2]  # the label for frustum is given as '2'
+                geomInfoList.append(primitiveLabel)
+                [bot_x, bot_y, bot_z] = extractCoord(botpoint)  # bottom point
+                l = len(geomInfoList)
+                geomInfoList[l - 1].append(bot_x)
+                geomInfoList[l - 1].append(bot_y)
+                geomInfoList[l - 1].append(bot_z)
+                [top_x, top_y, top_z] = extractCoord(frustumTopPt)  # top point
+                geomInfoList[l - 1].append(top_x)
+                geomInfoList[l - 1].append(top_y)
+                geomInfoList[l - 1].append(top_z)
+                geomInfoList[l - 1].append(frustumRad)  # base radius
+                geomInfoList[l - 1].append(topRad)  # top radius
+            ### When the primitive is a sphere ###
+            elif isinstance(geometry, Sphere):
+                sphereRad = geometry.radius
+                # print("Sphere with label:", lstring[id][1], "Radius:", sphereRad)
+                # Now we store the label and geometry information of the primitive as the following tuple:
+                # (primitive label, sphere radius)
+                # So the tuple contains: (primitive label, sphere radius)
+                primitiveLabel = [3]  # the label for sphere is given as '3'
+                geomInfoList.append(primitiveLabel)
+                l = len(geomInfoList)
+                geomInfoList[l - 1].append(sphereRad)
+            elif isinstance(geometry, Translated):
+                botpoint = geometry.translation
+            elif isinstance(geometry, Oriented):
+                heading = cross(geometry.primary, geometry.secondary)
+    # print(vertexList), print(triangleList), print(geomInfoList), print(globalVertexList), print(globalTriangleIdList), print(globalTriangleLabelList)
+    ##### So at this point, we are done with storing all the labelled geometry primitives and the labelled vertices of the coarse mesh model #####
+    ##### Next, we will do the point sampling to generate points on the surface of the model. #####
 
-  # loop over all elements in the lscene
-  for shape in lscene:
-    botpoint = Vector3(0,0,0) # the model is based at the origin
-    heading = Vector3(0,0,1)
-    geometry = shape
-    id = shape.id
-    currentModule = lstring[id].name
-    labelCurrentModule = labelTable[currentModule]
-    [ptList, idxList] = performTessalation(shape.geometry) # get the triangle corresponding to the primitive  
-    vertexList = obtainVertexList(ptList) # store the set of indices of vertices that form triangles (local index)
-    triangleList = obtainTriangleIdxList(idxList) # store the set of triangle indices
-    ### Now perform the global storing operations ###
-    globalVertexList = addVertexToGlobalList(vertexList, globalVertexList)
-    [globalTriangleIdList, globalTriangleLabelList] = addTriangleIdToGlobalList(vertexList, globalVertexList, triangleList, globalTriangleIdList, globalTriangleLabelList, [labelCurrentModule])
-    # the following loop runs over all the basic primitives we have considered in the model
-    while not isinstance(geometry,Primitive):
-      geometry = geometry.geometry
-      ### When the primitive is a cylinder ###
-      if isinstance(geometry,Cylinder):
-        cyllist.append((botpoint,botpoint+heading*geometry.height,geometry.radius))
-        toppoint = botpoint+heading*geometry.height
-        rad = geometry.radius
-        #print("Cylinder with label:", lstring[id][2], "Bottom point coordinate:", botpoint, "Top point coordinate:", toppoint, "Radius:",rad)
-        # Now we store the label and geometry information of the primitive as the following tuple:
-        # (primitive label, bottom_x, bottom_y, bottom_z, top_x, top_y, top_z, radius)
-        primitiveLabel = [1] # the label for cylinder is given as '1'
-        geomInfoList.append(primitiveLabel)
-        [bot_x,bot_y,bot_z] = extractCoord(botpoint) #bottom point
-        l = len(geomInfoList)
-        geomInfoList[l-1].append(bot_x)
-        geomInfoList[l-1].append(bot_y)
-        geomInfoList[l-1].append(bot_z)
-        [top_x,top_y,top_z] = extractCoord(toppoint) #top point
-        geomInfoList[l-1].append(top_x)
-        geomInfoList[l-1].append(top_y)
-        geomInfoList[l-1].append(top_z)
-        geomInfoList[l-1].append(rad) #radius
-      ### When the primitive is a frustum ###
-      elif isinstance(geometry,Frustum):
-        frustumTopPt = botpoint+heading*geometry.height
-        frustumRad = geometry.radius
-        taperFactor = geometry.taper
-        topRad = frustumRad * taperFactor
-        #print("Frustum with label:", lstring[id][2], "Bottom point coordinate:", botpoint, "Top point coordinate:", frustumTopPt, "Base radius:", frustumRad, "Top radius:", topRad)
-        # Now we store the label and geometry information of the primitive as the following tuple:
-        # (primitive label, bottom_x, bottom_y, bottom_z, top_x, top_y, top_z, base radius, top radius)
-        primitiveLabel = [2] # the label for frustum is given as '2'
-        geomInfoList.append(primitiveLabel)
-        [bot_x,bot_y,bot_z] = extractCoord(botpoint) #bottom point
-        l = len(geomInfoList)
-        geomInfoList[l-1].append(bot_x)
-        geomInfoList[l-1].append(bot_y)
-        geomInfoList[l-1].append(bot_z)
-        [top_x,top_y,top_z] = extractCoord(frustumTopPt) #top point
-        geomInfoList[l-1].append(top_x)
-        geomInfoList[l-1].append(top_y)
-        geomInfoList[l-1].append(top_z)
-        geomInfoList[l-1].append(frustumRad) #base radius
-        geomInfoList[l-1].append(topRad) #top radius
-      ### When the primitive is a sphere ###
-      elif isinstance(geometry,Sphere):
-        sphereRad = geometry.radius
-        #print("Sphere with label:", lstring[id][1], "Radius:", sphereRad)
-        # Now we store the label and geometry information of the primitive as the following tuple:
-        # (primitive label, sphere radius)
-        #So the tuple contains: (primitive label, sphere radius)
-        primitiveLabel = [3] # the label for sphere is given as '3'
-        geomInfoList.append(primitiveLabel)
-        l = len(geomInfoList)
-        geomInfoList[l-1].append(sphereRad)
-      elif isinstance(geometry,Translated):
-        botpoint = geometry.translation
-      elif isinstance(geometry,Oriented):
-        heading = cross(geometry.primary, geometry.secondary)
-  #print(vertexList), print(triangleList), print(geomInfoList), print(globalVertexList), print(globalTriangleIdList), print(globalTriangleLabelList)
-  ##### So at this point, we are done with storing all the labelled geometry primitives and the labelled vertices of the coarse mesh model #####
-  ##### Next, we will do the point sampling to generate points on the surface of the model. #####
-  
-  ######################## Basic point sampler ########################
+    ######################## Basic point sampler ########################
 
-  ### Generate exactly totalNumberOfResampledPts number of points ###
-  finalSampledPoints = []
-  finalRandomPointLabels = []
-  totalSampledPointsSoFar = 0
-  initialPointRequirement = totalNumberOfResampledPts
-  while(totalNumberOfResampledPts != totalSampledPointsSoFar):
-    randomPointLabels = []
-    randomTriangleList = selectRandomTriangle(globalTriangleIdList, globalVertexList, initialPointRequirement) # stores a list of random triangles with 'totalNumberOfResampledPts' number of elements. Next, we will sample one point per triangle to generate the point cloud
-    [newSampledPoints, randomPointLabels] = sampleRandomPoints(randomTriangleList, globalTriangleIdList, globalVertexList, initialPointRequirement, globalTriangleLabelList, randomPointLabels) # samples a random point on each triangle in the list 'randomTriangleList'
-    [ultimateSampledPoints, ultimateRandomPointLabels] = determineInsideness(geomInfoList, newSampledPoints, randomPointLabels) # perform insideness testing of the generated points & discard the points which are inside any primitive
-    currentNumberOfPoints = len(ultimateSampledPoints)
-    totalSampledPointsSoFar = totalSampledPointsSoFar + currentNumberOfPoints
-    pointDifference = totalNumberOfResampledPts - totalSampledPointsSoFar
-    initialPointRequirement = pointDifference
-    #finalSampledPoints = ultimateSampledPoints
-    #finalRandomPointLabels = ultimateRandomPointLabels
-    newPtLen = len(ultimateSampledPoints)
-    for l in range(newPtLen):
-      finalSampledPoints.append(ultimateSampledPoints[l])
-      finalRandomPointLabels.append(ultimateRandomPointLabels[l])
-    if(totalSampledPointsSoFar > totalNumberOfResampledPts):
-      break 
+    ### Generate exactly totalNumberOfResampledPts number of points ###
+    finalSampledPoints = []
+    finalRandomPointLabels = []
+    totalSampledPointsSoFar = 0
+    initialPointRequirement = totalNumberOfResampledPts
+    while (totalNumberOfResampledPts != totalSampledPointsSoFar):
+        randomPointLabels = []
+        randomTriangleList = selectRandomTriangle(globalTriangleIdList, globalVertexList,
+                                                  initialPointRequirement)  # stores a list of random triangles with 'totalNumberOfResampledPts' number of elements. Next, we will sample one point per triangle to generate the point cloud
+        [newSampledPoints, randomPointLabels] = sampleRandomPoints(randomTriangleList, globalTriangleIdList,
+                                                                   globalVertexList, initialPointRequirement,
+                                                                   globalTriangleLabelList,
+                                                                   randomPointLabels)  # samples a random point on each triangle in the list 'randomTriangleList'
+        [ultimateSampledPoints, ultimateRandomPointLabels] = determineInsideness(geomInfoList, newSampledPoints,
+                                                                                 randomPointLabels)  # perform insideness testing of the generated points & discard the points which are inside any primitive
+        currentNumberOfPoints = len(ultimateSampledPoints)
+        totalSampledPointsSoFar = totalSampledPointsSoFar + currentNumberOfPoints
+        pointDifference = totalNumberOfResampledPts - totalSampledPointsSoFar
+        initialPointRequirement = pointDifference
+        # finalSampledPoints = ultimateSampledPoints
+        # finalRandomPointLabels = ultimateRandomPointLabels
+        newPtLen = len(ultimateSampledPoints)
+        for l in range(newPtLen):
+            finalSampledPoints.append(ultimateSampledPoints[l])
+            finalRandomPointLabels.append(ultimateRandomPointLabels[l])
+        if (totalSampledPointsSoFar > totalNumberOfResampledPts):
+            break
 
-  [onlyPoints, onlyLabels,finalLabelledPoints] = createLabelledPointForDisplay(finalSampledPoints, finalRandomPointLabels) #this is after insidenss testing
-  write2File(pointsWithColor, finalLabelledPoints) # the data format is: x, y, z, r, g, b (each r,g,b combination is unique for a specific label)
-  write2File(rawPoints, onlyPoints) # the data format is: x, y, z
-  write2File(rawLabels, onlyLabels) # the data format is: label
+    [onlyPoints, onlyLabels, finalLabelledPoints] = createLabelledPointForDisplay(finalSampledPoints,
+                                                                                  finalRandomPointLabels)  # this is after insidenss testing
+    write2File(pointsWithColor,
+               finalLabelledPoints)  # the data format is: x, y, z, r, g, b (each r,g,b combination is unique for a specific label)
+    write2File(rawPoints, onlyPoints)  # the data format is: x, y, z
+    write2File(rawLabels, onlyLabels)  # the data format is: label
diff --git a/src/lpy_tools/point_sampler/utilities.py b/src/lpy_tools/point_sampler/utilities.py
index 2d32cd07cc726d06122f4041f267ac4987a78034..117b204a3fad5c7553ff2d3e4ccd0c18d0f5ef9c 100755
--- a/src/lpy_tools/point_sampler/utilities.py
+++ b/src/lpy_tools/point_sampler/utilities.py
@@ -1,348 +1,356 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
 ### This file contains helper functions for the point sampler ###
 ### Author: Ayan Chaudhury (ayanchaudhury.cs@gmail.com) ###
 ### INRIA team MOSAIC ###
 ### Updated: August 2021 ###
 
-from random import seed, random, randint, uniform, gauss
-from numpy import arange
-from math import *
-import numpy as np 
+import numpy as np
 
 
 # Extract floating point numbers from Vec3 format
 def extractCoord(vec3coord):
-  x = vec3coord[0]
-  y = vec3coord[1]
-  z = vec3coord[2]
-  return x,y,z
+    x = vec3coord[0]
+    y = vec3coord[1]
+    z = vec3coord[2]
+    return x, y, z
+
 
 # Computes the set of indices of vertices that form the triangles
 def obtainVertexList(pointlist):
-  vertexlist = []
-  totalVertices = len(pointlist) # get the no. of triangles in the primitive
-  for i in range(totalVertices):
-    coord = pointlist[i]
-    [x,y,z] = extractCoord(coord)
-    #Store the vertex coordinates to vertexlist
-    vertexlist.append([x])
-    l = len(vertexlist)
-    vertexlist[l-1].append(y)
-    vertexlist[l-1].append(z)
-  return vertexlist
+    vertexlist = []
+    totalVertices = len(pointlist)  # get the no. of triangles in the primitive
+    for i in range(totalVertices):
+        coord = pointlist[i]
+        [x, y, z] = extractCoord(coord)
+        # Store the vertex coordinates to vertexlist
+        vertexlist.append([x])
+        l = len(vertexlist)
+        vertexlist[l - 1].append(y)
+        vertexlist[l - 1].append(z)
+    return vertexlist
 
 
 # Extracts the vertices corresponding to each triangle
 def obtainTriangleIdxList(idxList):
-  triangleList = []
-  totalTriangles = len(idxList)
-  for i in range(totalTriangles):
-    triIndices = idxList[i]
-    [v1,v2,v3] = extractCoord(triIndices)
-    triangleList.append([v1])
-    l = len(triangleList)
-    triangleList[l-1].append(v2)
-    triangleList[l-1].append(v3)
-  return triangleList
+    triangleList = []
+    totalTriangles = len(idxList)
+    for i in range(totalTriangles):
+        triIndices = idxList[i]
+        [v1, v2, v3] = extractCoord(triIndices)
+        triangleList.append([v1])
+        l = len(triangleList)
+        triangleList[l - 1].append(v2)
+        triangleList[l - 1].append(v3)
+    return triangleList
 
 
 # Local triangle vertices are added to the global list
 def addVertexToGlobalList(localVertexList, globalVertexList):
-  localLen = len(localVertexList)
-  for i in range(localLen):
-    vertexExists = localVertexList[i] in globalVertexList
-    if (vertexExists == False):
-      globalVertexList.append(localVertexList[i])
-  return globalVertexList
+    localLen = len(localVertexList)
+    for i in range(localLen):
+        vertexExists = localVertexList[i] in globalVertexList
+        if (vertexExists == False):
+            globalVertexList.append(localVertexList[i])
+    return globalVertexList
 
 
 # Local triangle id's are  added to the global list 
-def addTriangleIdToGlobalList(localVertexList, globalVertexList, localTriangleIdList, globalTriangleIdList, globalTriangleLabelList, currentLabel):
-  totalTriangles = len(localTriangleIdList)
-  totalLocalVertex = len(localVertexList)
-  totalGlobalVertex = len(globalVertexList)
-  for i in range(totalTriangles):
-    currentTriangle = localTriangleIdList[i] #is a list
-    v1_local_id = currentTriangle[0] #is a number
-    v2_local_id = currentTriangle[1]
-    v3_local_id = currentTriangle[2]
-    v1_coord_local = localVertexList[v1_local_id] #is a list
-    v2_coord_local = localVertexList[v2_local_id]
-    v3_coord_local = localVertexList[v3_local_id]
-    v1_idx_Global = globalVertexList.index(v1_coord_local)
-    v2_idx_Global = globalVertexList.index(v2_coord_local)
-    v3_idx_Global = globalVertexList.index(v3_coord_local)
-    currentTriangleGlobalId = [v1_idx_Global, v2_idx_Global, v3_idx_Global]
-    globalTriangleIdList.append(currentTriangleGlobalId)
-    globalTriangleLabelList.append(currentLabel) # add the triangle label
-  return globalTriangleIdList, globalTriangleLabelList
-  
+def addTriangleIdToGlobalList(localVertexList, globalVertexList, localTriangleIdList, globalTriangleIdList,
+                              globalTriangleLabelList, currentLabel):
+    totalTriangles = len(localTriangleIdList)
+    totalLocalVertex = len(localVertexList)
+    totalGlobalVertex = len(globalVertexList)
+    for i in range(totalTriangles):
+        currentTriangle = localTriangleIdList[i]  # is a list
+        v1_local_id = currentTriangle[0]  # is a number
+        v2_local_id = currentTriangle[1]
+        v3_local_id = currentTriangle[2]
+        v1_coord_local = localVertexList[v1_local_id]  # is a list
+        v2_coord_local = localVertexList[v2_local_id]
+        v3_coord_local = localVertexList[v3_local_id]
+        v1_idx_Global = globalVertexList.index(v1_coord_local)
+        v2_idx_Global = globalVertexList.index(v2_coord_local)
+        v3_idx_Global = globalVertexList.index(v3_coord_local)
+        currentTriangleGlobalId = [v1_idx_Global, v2_idx_Global, v3_idx_Global]
+        globalTriangleIdList.append(currentTriangleGlobalId)
+        globalTriangleLabelList.append(currentLabel)  # add the triangle label
+    return globalTriangleIdList, globalTriangleLabelList
+
 
 # Write the given list to a file
 def write2File(myfilename, mylist):
-  import csv
-  with open(myfilename, 'w') as f:
-    csv_writer = csv.writer(f)
-    csv_writer.writerows(mylist)
-    
+    import csv
+    with open(myfilename, 'w') as f:
+        csv_writer = csv.writer(f)
+        csv_writer.writerows(mylist)
+
 
 # Based on the camera position, the distance factor is computed for each triangle. The center of mass is
 # chosen as the representative point inside each triangle (although that is not a good idea for big triangles,
 # and other types of strategies to have a gradient effect inside each triangle will be better). The function
-# returns the distances in the form of probabilties.
+# returns the distances in the form of probabilities.
 def computeDistanceFactor(vertices_1, vertices_2, vertices_3, cameraPosition):
-  #cameraPosition = np.asarray([73.0, -69.0, 31.5])
-  #cameraPosition = np.asarray([0.0, 0.0, 0.0]) # use custom camera position according to the need of the application
-  midpoints = (vertices_1 + vertices_2 + vertices_3) / 3
-  distances = np.linalg.norm((cameraPosition - midpoints), axis = 1)
-  maxdist = max(distances)
-  mindist = min(distances)
-  distances = (distances - mindist) / (maxdist - mindist)
-  distances = 1 - distances ### small distances have high score, long distances have low score ###
-  return distances
-    
+    # cameraPosition = np.asarray([73.0, -69.0, 31.5])
+    # cameraPosition = np.asarray([0.0, 0.0, 0.0]) # use custom camera position according to the need of the application
+    midpoints = (vertices_1 + vertices_2 + vertices_3) / 3
+    distances = np.linalg.norm((cameraPosition - midpoints), axis=1)
+    maxdist = max(distances)
+    mindist = min(distances)
+    distances = (distances - mindist) / (maxdist - mindist)
+    distances = 1 - distances  ### small distances have high score, long distances have low score ###
+    return distances
 
 
 # This function computes the number of points to be sampled in each triangle. Initially the number of points per
 # square unit area is computed, which is multiplied by the area of the triangle under consideration. So that gives
 # the number of points to be sampled for uniform sampling case. Since we aim at having a distance factor, the factor
 # is weighted by a distance probability - small distances yield large weightage, and big distances yield small weightage.
-# The ceil() calculations are done because in many cases there are many tiny triabgles and the total number of estimated
+# The ceil() calculations are done because in many cases there are many tiny triangles and the total number of estimated
 # points inside the triangles are computed as < 1.0. Performing ceil() operation prevents us in having empty triangles.
 # PS- the number of sampled points should be specified as large in these cases, then only we can see the distance effect.
 def distanceSampledTriangles(triangleIdList, distanceList, areaList, totalArea, meanPoints):
-  from math import floor, ceil, exp
-  newTriangleList = []
-  perSqAreaPoints_uniform = ceil(meanPoints / totalArea)
-  totalTriangles = len(triangleIdList)
-  for i in range(totalTriangles):
-    currentArea = areaList[i]
-    currentDistance = distanceList[i]
-    perSqAreaPoints_new = int(ceil(perSqAreaPoints_uniform * currentArea * currentDistance * currentDistance)) #square distance function
-    #perSqAreaPoints_new = int(ceil(perSqAreaPoints_uniform * currentArea * exp(-currentDistance))) # exponential distance function
-    for j in range(perSqAreaPoints_new):
-      newTriangleList.append(i)
-  npNewTriangleList = np.asarray(newTriangleList)
-  return npNewTriangleList
-
+    from math import ceil
+    newTriangleList = []
+    perSqAreaPoints_uniform = ceil(meanPoints / totalArea)
+    totalTriangles = len(triangleIdList)
+    for i in range(totalTriangles):
+        currentArea = areaList[i]
+        currentDistance = distanceList[i]
+        perSqAreaPoints_new = int(
+            ceil(perSqAreaPoints_uniform * currentArea * currentDistance * currentDistance))  # square distance function
+        # perSqAreaPoints_new = int(ceil(perSqAreaPoints_uniform * currentArea * exp(-currentDistance))) # exponential distance function
+        for j in range(perSqAreaPoints_new):
+            newTriangleList.append(i)
+    npNewTriangleList = np.asarray(newTriangleList)
+    return npNewTriangleList
 
 
 # Triangle selection in advanced point sampler where the distance effect is modelled. Triangles near the camera are more
 # frequently selected than the triangles which are far away. 
-def selectRandomTriangleWithDistanceEffect(triangleIdList, triangleVertexList, totalNumberOfResampledPts, cameraPosition):
-  v1 = []
-  v2 = []
-  v3 = []
-  totalTriangles = len(triangleIdList)
-  for i in range(totalTriangles):
-    currentTriangle = triangleIdList[i]
-    v1_local_id = currentTriangle[0]
-    v2_local_id = currentTriangle[1]
-    v3_local_id = currentTriangle[2]
-    v1_coordinate = triangleVertexList[v1_local_id]
-    v2_coordinate = triangleVertexList[v2_local_id]
-    v3_coordinate = triangleVertexList[v3_local_id]
-    v1.append(v1_coordinate)
-    v2.append(v2_coordinate)
-    v3.append(v3_coordinate)
-  v1np = np.asarray(v1) # just a data type conversion for the ease of calculations
-  v2np = np.asarray(v2)
-  v3np = np.asarray(v3)
-  allTriangleAreas = 0.5 * np.linalg.norm(np.cross(v2np - v1np, v3np - v1np), axis = 1)
-  totalArea = allTriangleAreas.sum() # total area of all the triangles (total surface area)
-  distances = computeDistanceFactor(v1np, v2np, v3np, cameraPosition) # close points to the camera will have low value
-  newSampledTriangles = distanceSampledTriangles(triangleIdList, distances, allTriangleAreas, totalArea, totalNumberOfResampledPts)
-  return newSampledTriangles
+def selectRandomTriangleWithDistanceEffect(triangleIdList, triangleVertexList, totalNumberOfResampledPts,
+                                           cameraPosition):
+    v1 = []
+    v2 = []
+    v3 = []
+    totalTriangles = len(triangleIdList)
+    for i in range(totalTriangles):
+        currentTriangle = triangleIdList[i]
+        v1_local_id = currentTriangle[0]
+        v2_local_id = currentTriangle[1]
+        v3_local_id = currentTriangle[2]
+        v1_coordinate = triangleVertexList[v1_local_id]
+        v2_coordinate = triangleVertexList[v2_local_id]
+        v3_coordinate = triangleVertexList[v3_local_id]
+        v1.append(v1_coordinate)
+        v2.append(v2_coordinate)
+        v3.append(v3_coordinate)
+    v1np = np.asarray(v1)  # just a data type conversion for the ease of calculations
+    v2np = np.asarray(v2)
+    v3np = np.asarray(v3)
+    allTriangleAreas = 0.5 * np.linalg.norm(np.cross(v2np - v1np, v3np - v1np), axis=1)
+    totalArea = allTriangleAreas.sum()  # total area of all the triangles (total surface area)
+    distances = computeDistanceFactor(v1np, v2np, v3np,
+                                      cameraPosition)  # close points to the camera will have low value
+    newSampledTriangles = distanceSampledTriangles(triangleIdList, distances, allTriangleAreas, totalArea,
+                                                   totalNumberOfResampledPts)
+    return newSampledTriangles
 
 
 # Random triangle selection in basic point sampler. Among a list of triangles, random selection is performed based on the
 # probability of triangle area. This allows us not to over-sample from areas having large number of small triangles.
 def selectRandomTriangle(triangleIdList, triangleVertexList, totalNumberOfResampledPts):
-  v1 = []
-  v2 = []
-  v3 = []
-  totalTriangles = len(triangleIdList)
-  for i in range(totalTriangles):
-    currentTriangle = triangleIdList[i]
-    v1_local_id = currentTriangle[0]
-    v2_local_id = currentTriangle[1]
-    v3_local_id = currentTriangle[2]
-    v1_coordinate = triangleVertexList[v1_local_id]
-    v2_coordinate = triangleVertexList[v2_local_id]
-    v3_coordinate = triangleVertexList[v3_local_id]
-    v1.append(v1_coordinate)
-    v2.append(v2_coordinate)
-    v3.append(v3_coordinate)
-  v1np = np.asarray(v1) # just a data type conversion for the ease of calculations
-  v2np = np.asarray(v2)
-  v3np = np.asarray(v3)
-  allTriangleAreas = 0.5 * np.linalg.norm(np.cross(v2np - v1np, v3np - v1np), axis = 1)
-  weightedAreas = allTriangleAreas / allTriangleAreas.sum()
-  randomTriangleList = np.random.choice(range(totalTriangles), size = totalNumberOfResampledPts, p = weightedAreas)
-  return randomTriangleList
-  
+    v1 = []
+    v2 = []
+    v3 = []
+    totalTriangles = len(triangleIdList)
+    for i in range(totalTriangles):
+        currentTriangle = triangleIdList[i]
+        v1_local_id = currentTriangle[0]
+        v2_local_id = currentTriangle[1]
+        v3_local_id = currentTriangle[2]
+        v1_coordinate = triangleVertexList[v1_local_id]
+        v2_coordinate = triangleVertexList[v2_local_id]
+        v3_coordinate = triangleVertexList[v3_local_id]
+        v1.append(v1_coordinate)
+        v2.append(v2_coordinate)
+        v3.append(v3_coordinate)
+    v1np = np.asarray(v1)  # just a data type conversion for the ease of calculations
+    v2np = np.asarray(v2)
+    v3np = np.asarray(v3)
+    allTriangleAreas = 0.5 * np.linalg.norm(np.cross(v2np - v1np, v3np - v1np), axis=1)
+    weightedAreas = allTriangleAreas / allTriangleAreas.sum()
+    randomTriangleList = np.random.choice(range(totalTriangles), size=totalNumberOfResampledPts, p=weightedAreas)
+    return randomTriangleList
+
 
 # Randomly chooses 3 variables in the range [0,1]
 # such that their sum equals to 1 
 def generateAlphaBetaGamma():
-  alpha = np.random.rand(1)
-  beta = np.random.rand(1)
-  gamma = np.random.rand(1)
-  totalSum = alpha + beta + gamma
-  newAlpha = alpha / totalSum
-  newBeta = beta / totalSum
-  newGamma = gamma / totalSum
-  return newAlpha, newBeta, newGamma
-  
+    alpha = np.random.rand(1)
+    beta = np.random.rand(1)
+    gamma = np.random.rand(1)
+    totalSum = alpha + beta + gamma
+    newAlpha = alpha / totalSum
+    newBeta = beta / totalSum
+    newGamma = gamma / totalSum
+    return newAlpha, newBeta, newGamma
+
 
 # Given a list of triangles, this function generates one point inside each triangle in a random position,
 # and returns sampled points along with its label obtained from the label of the triangles in which the points are sampled from
-def sampleRandomPoints(randomTriangleIdList, globalTriangleIdList, globalTriangleVertexList, totalNumberOfResampledPts, globalTriangleLabelList, randomPointLabels):
-  allSampledPoints = []
-  totalNumberOfResampledPts = randomTriangleIdList.shape[0]
-  for i in range(totalNumberOfResampledPts):
-    currentRandomTriangle = randomTriangleIdList[i]
-    currentLabel = globalTriangleLabelList[currentRandomTriangle][0]
-    currentTriangleVertexIds = globalTriangleIdList[currentRandomTriangle] #a list
-    v1_id = currentTriangleVertexIds[0]
-    v2_id = currentTriangleVertexIds[1]
-    v3_id = currentTriangleVertexIds[2]
-    v1_corrd = globalTriangleVertexList[v1_id] #a list
-    v2_corrd = globalTriangleVertexList[v2_id]
-    v3_corrd = globalTriangleVertexList[v3_id]
-    [alpha, beta, gamma] = generateAlphaBetaGamma()
-    sampledPoint = [sum(x) for x in zip([i * alpha for i in v1_corrd],[i * beta for i in v2_corrd],[i * gamma for i in v3_corrd])]
-    sampledPoint_x = sampledPoint[0][0] #because these are in the form of list of list
-    sampledPoint_y = sampledPoint[1][0]
-    sampledPoint_z = sampledPoint[2][0]
-    newPoint = [sampledPoint_x, sampledPoint_y, sampledPoint_z]
-    allSampledPoints.append(newPoint)
-    randomPointLabels.append(currentLabel)
-  return allSampledPoints, randomPointLabels
+def sampleRandomPoints(randomTriangleIdList, globalTriangleIdList, globalTriangleVertexList, totalNumberOfResampledPts,
+                       globalTriangleLabelList, randomPointLabels):
+    allSampledPoints = []
+    totalNumberOfResampledPts = randomTriangleIdList.shape[0]
+    for i in range(totalNumberOfResampledPts):
+        currentRandomTriangle = randomTriangleIdList[i]
+        currentLabel = globalTriangleLabelList[currentRandomTriangle][0]
+        currentTriangleVertexIds = globalTriangleIdList[currentRandomTriangle]  # a list
+        v1_id = currentTriangleVertexIds[0]
+        v2_id = currentTriangleVertexIds[1]
+        v3_id = currentTriangleVertexIds[2]
+        v1_corrd = globalTriangleVertexList[v1_id]  # a list
+        v2_corrd = globalTriangleVertexList[v2_id]
+        v3_corrd = globalTriangleVertexList[v3_id]
+        [alpha, beta, gamma] = generateAlphaBetaGamma()
+        sampledPoint = [sum(x) for x in
+                        zip([i * alpha for i in v1_corrd], [i * beta for i in v2_corrd], [i * gamma for i in v3_corrd])]
+        sampledPoint_x = sampledPoint[0][0]  # because these are in the form of list of list
+        sampledPoint_y = sampledPoint[1][0]
+        sampledPoint_z = sampledPoint[2][0]
+        newPoint = [sampledPoint_x, sampledPoint_y, sampledPoint_z]
+        allSampledPoints.append(newPoint)
+        randomPointLabels.append(currentLabel)
+    return allSampledPoints, randomPointLabels
 
 
 # Returns unique elements of a list
-def distinctElements(list1): 
-  distinct_list = []
-  for x in list1:
-    if x not in distinct_list:
-      distinct_list.append(x)
-  return distinct_list
+def distinctElements(list1):
+    distinct_list = []
+    for x in list1:
+        if x not in distinct_list:
+            distinct_list.append(x)
+    return distinct_list
 
 
 # Given a list of labels, this function performs relabelling sequentially starting
 # from 1. So if a given list of label is [30, 31, 37, 37, 37, 42, 42, 51, 51, 51],
 # then the new labelling is obtained as, [1, 2, 3, 3, 3, 4, 4, 5, 5, 5]
 def refineLabels(a):
-  newLabelList = []
-  b = distinctElements(a)
-  b.sort()
-  totalLabels = len(b)
-  labelSeq = np.linspace(1, totalLabels, totalLabels)
-  labelSeq = [int(i) for i in labelSeq]
-  for i in range(len(a)):
-    for j in range(len(b)):
-        if (a[i] == b[j]):
-            newLabelList.append(labelSeq[j])
-  return newLabelList
+    newLabelList = []
+    b = distinctElements(a)
+    b.sort()
+    totalLabels = len(b)
+    labelSeq = np.linspace(1, totalLabels, totalLabels)
+    labelSeq = [int(i) for i in labelSeq]
+    for i in range(len(a)):
+        for j in range(len(b)):
+            if (a[i] == b[j]):
+                newLabelList.append(labelSeq[j])
+    return newLabelList
 
 
 # This function takes as argument the list of points and labels, and returns
 # 3 types of formats: (x,y,z,r,g,b), (x,y,z), (labels). The r,g,b triplets are
 # generated uniquely for different labels
 def createLabelledPointForDisplay(newSampledPoints, randomPointLabels):
-  refinedLabels = refineLabels(randomPointLabels)
-  allLabelledPoints = []
-  onlyPoints = []
-  onlyLabels = []
-  totalPoints = len(newSampledPoints)
-  #We keep the first 3 label colours as r, g, b & the rest are random colours
-  color1 = [255, 0, 0] + list(np.random.choice(range(256), size=2000)) #label max up to 2000
-  color2 = [0, 255, 0] + list(np.random.choice(range(256), size=2000))
-  color3 = [0, 0, 255] + list(np.random.choice(range(256), size=2000))
-  for i in range(totalPoints):
-    currentPoint = newSampledPoints[i]
-    currentLabel = refinedLabels[i]
-    color = [color1[currentLabel-1], color2[currentLabel-1], color3[currentLabel-1]] # labels start from '1', array idx from '0'
-    currentLabelledPoint = currentPoint + color
-    allLabelledPoints.append(currentLabelledPoint)
-    onlyPoints.append(currentPoint)
-    onlyLabels.append([currentLabel])
-  return onlyPoints, onlyLabels, allLabelledPoints
-  
+    refinedLabels = refineLabels(randomPointLabels)
+    allLabelledPoints = []
+    onlyPoints = []
+    onlyLabels = []
+    totalPoints = len(newSampledPoints)
+    # We keep the first 3 label colours as r, g, b & the rest are random colours
+    color1 = [255, 0, 0] + list(np.random.choice(range(256), size=2000))  # label max up to 2000
+    color2 = [0, 255, 0] + list(np.random.choice(range(256), size=2000))
+    color3 = [0, 0, 255] + list(np.random.choice(range(256), size=2000))
+    for i in range(totalPoints):
+        currentPoint = newSampledPoints[i]
+        currentLabel = refinedLabels[i]
+        color = [color1[currentLabel - 1], color2[currentLabel - 1],
+                 color3[currentLabel - 1]]  # labels start from '1', array idx from '0'
+        currentLabelledPoint = currentPoint + color
+        allLabelledPoints.append(currentLabelledPoint)
+        onlyPoints.append(currentPoint)
+        onlyLabels.append([currentLabel])
+    return onlyPoints, onlyLabels, allLabelledPoints
 
 
 # This is a simple trick used for the ease of floating point calculations,
 # especially to handle rounding off problems
 def applyFloatingPointTrick(n1, n2, n):
-  from math import floor, ceil
-  if (0.1 <= n < 10.0):
-    n1 = floor(n1 * 10) / 10.0
-    n2 = ceil(n2 * 10) / 10.0
-  elif(0.01 <= n < 0.1):
-    n1 = floor(n1 * 100) / 100.0
-    n2 = ceil(n2 * 100) / 100.0
-  elif(0.001 <= n < 0.01):
-    n1 = floor(n1 * 1000) / 1000.0
-    n2 = ceil(n2 * 1000) / 1000.0
-  return n1, n2
+    from math import floor, ceil
+    if (0.1 <= n < 10.0):
+        n1 = floor(n1 * 10) / 10.0
+        n2 = ceil(n2 * 10) / 10.0
+    elif (0.01 <= n < 0.1):
+        n1 = floor(n1 * 100) / 100.0
+        n2 = ceil(n2 * 100) / 100.0
+    elif (0.001 <= n < 0.01):
+        n1 = floor(n1 * 1000) / 1000.0
+        n2 = ceil(n2 * 1000) / 1000.0
+    return n1, n2
 
 
 # Given a list of geometric primitives in a model, a list of points and their corresponding labels,
 # this function checks which points are falling inside any of the primitives. These points are then
 # discarded and rest of the points along with their labels are returned
 def determineInsideness(geomInfoList, allPointsList, allPointLabels):
-  from math import floor, ceil
-  selectedPointsList = []
-  selectedPointsLabels = []
-  totalPoints = len(allPointsList)
-  totalPrimitives = len(geomInfoList)
-  for i in range(totalPoints):
-    currentPoint = np.asarray(allPointsList[i])
-    insidenessFlag = False
-    for j in range(totalPrimitives):
-      currentPrimitive = geomInfoList[j]
-      #cylinder primitive (we assigned it as label '1' in the End() function)
-      if(currentPrimitive[0] == 1):
-        bottomPoint = np.asarray([currentPrimitive[1], currentPrimitive[2], currentPrimitive[3]])
-        topPoint = np.asarray([currentPrimitive[4], currentPrimitive[5], currentPrimitive[6]])
-        rad = currentPrimitive[7]
-        vec = topPoint - bottomPoint
-        c = rad * np.linalg.norm(vec) 
-        condition1 = np.dot(currentPoint - bottomPoint, vec) >= 0
-        condition2 = np.dot(currentPoint - topPoint, vec) <= 0
-        condition3 = np.linalg.norm(np.cross(currentPoint - bottomPoint, vec))
-        minimum_number = min(c, condition3)
-        [c, condition3] = applyFloatingPointTrick(c, condition3, minimum_number)
-        z = (condition1 and condition2 and condition3 < c)
-        #if z returns true, that means the point is inside the cylinder
-        if(z == True):
-          insidenessFlag = True
-      #frustum primitive (we assigned it as label '2' in the End() function)
-      elif(currentPrimitive[0] == 2):
-        bottomPoint = np.asarray([currentPrimitive[1], currentPrimitive[2], currentPrimitive[3]])
-        topPoint = np.asarray([currentPrimitive[4], currentPrimitive[5], currentPrimitive[6]])
-        bottomRad = currentPrimitive[7]
-        topRad = currentPrimitive[8]
-        pointToAxisVector = currentPoint - bottomPoint
-        axisVector = topPoint - bottomPoint
-        condition1 = np.dot(currentPoint - bottomPoint, axisVector) >= 0
-        condition2 = np.dot(currentPoint - topPoint, axisVector) <= 0
-        projectedPoint = bottomPoint + np.dot(pointToAxisVector, axisVector) / np.dot(axisVector, axisVector) * axisVector
-        factor1 = np.linalg.norm(projectedPoint - bottomPoint) / np.linalg.norm(bottomPoint - topPoint) * bottomRad
-        factor2 = np.linalg.norm(projectedPoint - topPoint) / np.linalg.norm(bottomPoint - topPoint) * topRad
-        currentRadius = factor1 + factor2
-        currentDistance = np.linalg.norm(projectedPoint - currentPoint)
-        minimum_number = min(currentRadius, currentDistance)
-        [currentRadius, currentDistance] = applyFloatingPointTrick(currentRadius, currentDistance, minimum_number)
-        #if the distance is < frustum radius at that point, then it is inside
-        if (condition1 and condition2 and currentDistance < currentRadius):
-          insidenessFlag = True
-      # If there is sphere or other primitives in the model, then put it here as another condition(s)
-    if (insidenessFlag == False):
-      selectedPointsList.append(allPointsList[i])
-      selectedPointsLabels.append(allPointLabels[i])
-  return selectedPointsList, selectedPointsLabels
-  
+    selectedPointsList = []
+    selectedPointsLabels = []
+    totalPoints = len(allPointsList)
+    totalPrimitives = len(geomInfoList)
+    for i in range(totalPoints):
+        currentPoint = np.asarray(allPointsList[i])
+        insidenessFlag = False
+        for j in range(totalPrimitives):
+            currentPrimitive = geomInfoList[j]
+            # cylinder primitive (we assigned it as label '1' in the End() function)
+            if (currentPrimitive[0] == 1):
+                bottomPoint = np.asarray([currentPrimitive[1], currentPrimitive[2], currentPrimitive[3]])
+                topPoint = np.asarray([currentPrimitive[4], currentPrimitive[5], currentPrimitive[6]])
+                rad = currentPrimitive[7]
+                vec = topPoint - bottomPoint
+                c = rad * np.linalg.norm(vec)
+                condition1 = np.dot(currentPoint - bottomPoint, vec) >= 0
+                condition2 = np.dot(currentPoint - topPoint, vec) <= 0
+                condition3 = np.linalg.norm(np.cross(currentPoint - bottomPoint, vec))
+                minimum_number = min(c, condition3)
+                [c, condition3] = applyFloatingPointTrick(c, condition3, minimum_number)
+                z = (condition1 and condition2 and condition3 < c)
+                # if z returns true, that means the point is inside the cylinder
+                if (z == True):
+                    insidenessFlag = True
+            # frustum primitive (we assigned it as label '2' in the End() function)
+            elif (currentPrimitive[0] == 2):
+                bottomPoint = np.asarray([currentPrimitive[1], currentPrimitive[2], currentPrimitive[3]])
+                topPoint = np.asarray([currentPrimitive[4], currentPrimitive[5], currentPrimitive[6]])
+                bottomRad = currentPrimitive[7]
+                topRad = currentPrimitive[8]
+                pointToAxisVector = currentPoint - bottomPoint
+                axisVector = topPoint - bottomPoint
+                condition1 = np.dot(currentPoint - bottomPoint, axisVector) >= 0
+                condition2 = np.dot(currentPoint - topPoint, axisVector) <= 0
+                projectedPoint = bottomPoint + np.dot(pointToAxisVector, axisVector) / np.dot(axisVector,
+                                                                                              axisVector) * axisVector
+                factor1 = np.linalg.norm(projectedPoint - bottomPoint) / np.linalg.norm(
+                    bottomPoint - topPoint) * bottomRad
+                factor2 = np.linalg.norm(projectedPoint - topPoint) / np.linalg.norm(bottomPoint - topPoint) * topRad
+                currentRadius = factor1 + factor2
+                currentDistance = np.linalg.norm(projectedPoint - currentPoint)
+                minimum_number = min(currentRadius, currentDistance)
+                [currentRadius, currentDistance] = applyFloatingPointTrick(currentRadius, currentDistance,
+                                                                           minimum_number)
+                # if the distance is < frustum radius at that point, then it is inside
+                if (condition1 and condition2 and currentDistance < currentRadius):
+                    insidenessFlag = True
+            # If there is sphere or other primitives in the model, then put it here as another condition(s)
+        if (insidenessFlag == False):
+            selectedPointsList.append(allPointsList[i])
+            selectedPointsLabels.append(allPointLabels[i])
+    return selectedPointsList, selectedPointsLabels
+
 
 # Modelling plain noise. Given a list of points, corresponding labels, and all the geometric
 # primitives of the model, this function applies random normal distribution to the point position.
@@ -352,51 +360,51 @@ def determineInsideness(geomInfoList, allPointsList, allPointLabels):
 # works fine. However, other alternatives can also be modelled (or simply the number of trials can be
 # increased to have higher chances). The function returns noisy points and their labels.
 def convertToNoisyData(geomInfo, points, labels):
-  noisyPoints = []
-  noisyPointLabels = []
-  sigma = 0.07
-  n = len(points)
-  for i in range(n):
-    for j in range(5):
-      currentPoint = np.asarray(points[i])
-      x = np.random.normal(currentPoint[0], sigma)
-      y = np.random.normal(currentPoint[1], sigma)
-      z = np.random.normal(currentPoint[2], sigma)
-      noisyPoint = [[x, y, z]]
-      currentLabel = [labels[i]]
-      [insideTestedPts, insideTestedLabels] = determineInsideness(geomInfo, noisyPoint, currentLabel)
-      if(len(insideTestedPts) > 0):
-        noisyPoints.append(insideTestedPts[0])
-        noisyPointLabels.append(insideTestedLabels[0])
-        break
-  return noisyPoints, noisyPointLabels
+    noisyPoints = []
+    noisyPointLabels = []
+    sigma = 0.07
+    n = len(points)
+    for i in range(n):
+        for j in range(5):
+            currentPoint = np.asarray(points[i])
+            x = np.random.normal(currentPoint[0], sigma)
+            y = np.random.normal(currentPoint[1], sigma)
+            z = np.random.normal(currentPoint[2], sigma)
+            noisyPoint = [[x, y, z]]
+            currentLabel = [labels[i]]
+            [insideTestedPts, insideTestedLabels] = determineInsideness(geomInfo, noisyPoint, currentLabel)
+            if (len(insideTestedPts) > 0):
+                noisyPoints.append(insideTestedPts[0])
+                noisyPointLabels.append(insideTestedLabels[0])
+                break
+    return noisyPoints, noisyPointLabels
 
 
 # Modelling sensor noise. As a point gets further away from the camera, the noise increases. The insideness
 # is tested following the similar strategy as done in plain noise modelling. The function returns noisy points
 # and their labels.
 def convertToSensorNoisyData(geomInfo, points, labels, cameraPosition):
-  cameraPosition = np.asarray([0.0, 0.0, 0.0]) # use custom camera position here
-  noisyPoints = []
-  noisyPointLabels = []
-  maxSigma = 0.1 # maximum sigma occus to the furthest point(s)
-  distances = np.linalg.norm((cameraPosition - points), axis = 1)
-  maxdist = max(distances)
-  mindist = min(distances)
-  distances = (distances - mindist) / (maxdist - mindist) # normalize distance in the range [0,1]
-  n = len(points)
-  for i in range(n):
-    for j in range(5):
-      currentPoint = np.asarray(points[i])
-      currentSigma = maxSigma * (distances[i]) # small distances have lower effect than large distances
-      x = np.random.normal(currentPoint[0], currentSigma)
-      y = np.random.normal(currentPoint[1], currentSigma)
-      z = np.random.normal(currentPoint[2], currentSigma)
-      noisyPoint = [[x, y, z]]
-      currentLabel = [labels[i]]
-      [insideTestedPts, insideTestedLabels] = determineInsideness(geomInfo, noisyPoint, currentLabel)
-      if(len(insideTestedPts) > 0):
-        noisyPoints.append(insideTestedPts[0])
-        noisyPointLabels.append(insideTestedLabels[0])
-        break
-  return noisyPoints, noisyPointLabels
+    cameraPosition = np.asarray([0.0, 0.0, 0.0])  # use custom camera position here
+    noisyPoints = []
+    noisyPointLabels = []
+    maxSigma = 0.1  # maximum sigma occurs to the furthest point(s)
+    distances = np.linalg.norm((cameraPosition - points), axis=1)
+    maxdist = max(distances)
+    mindist = min(distances)
+    distances = (distances - mindist) / (maxdist - mindist)  # normalize distance in the range [0,1]
+    n = len(points)
+    for i in range(n):
+        for j in range(5):
+            currentPoint = np.asarray(points[i])
+            currentSigma = maxSigma * (distances[i])  # small distances have lower effect than large distances
+            x = np.random.normal(currentPoint[0], currentSigma)
+            y = np.random.normal(currentPoint[1], currentSigma)
+            z = np.random.normal(currentPoint[2], currentSigma)
+            noisyPoint = [[x, y, z]]
+            currentLabel = [labels[i]]
+            [insideTestedPts, insideTestedLabels] = determineInsideness(geomInfo, noisyPoint, currentLabel)
+            if (len(insideTestedPts) > 0):
+                noisyPoints.append(insideTestedPts[0])
+                noisyPointLabels.append(insideTestedLabels[0])
+                break
+    return noisyPoints, noisyPointLabels
diff --git a/src/lpy_tools/snapshots/__init__.py b/src/lpy_tools/snapshots/__init__.py
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..56fafa58b3f43decb7699b93048b8b87e0f695aa 100644
--- a/src/lpy_tools/snapshots/__init__.py
+++ b/src/lpy_tools/snapshots/__init__.py
@@ -0,0 +1,2 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
diff --git a/src/lpy_tools/snapshots/lscene_snapshots.py b/src/lpy_tools/snapshots/lscene_snapshots.py
index 7be60d242184bda2a22cac73721abc7a83ae6fb6..1c81d9442989d826d12d361aca56a80a60d9c59e 100644
--- a/src/lpy_tools/snapshots/lscene_snapshots.py
+++ b/src/lpy_tools/snapshots/lscene_snapshots.py
@@ -1,4 +1,4 @@
-# -*- python -*-
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
 #       File author(s):
@@ -13,37 +13,41 @@
 #
 # -----------------------------------------------------------------------
 
-'''
-    Library of tools to take snapshots of the result of a l-system simulations
+"""
+Library of tools to take snapshots of the result of a l-system simulations
 
-    Authors: Christophe Godin, Inria, 2021
-    Licence: Open source LGPL
+Authors: Christophe Godin, Inria, 2021
+Licence: Open source LGPL
 
-    The lscene_snapshot library offers 2 main functions:
-      simulate_and_shoot(model_filename, variables, suffix, cam_settings = camdict)
-      grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed_variables_dict, cam_settings = cam_setting_fun1, short_names = variable_short_names)
+The lscene_snapshot library offers 2 main functions:
+  simulate_and_shoot(model_filename, variables, suffix, cam_settings = camdict)
+  grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed_variables_dict, cam_settings = cam_setting_fun1, short_names = variable_short_names)
 
-'''
+"""
 
 import os
-from numpy import arange
-from openalea.plantgl.all import *
-from openalea import lpy
 from math import pi
 
-from matplotlib.offsetbox import TextArea, DrawingArea, OffsetImage, AnnotationBbox
-
-import matplotlib.pyplot as plt
 import matplotlib.image as mpimg
+import matplotlib.pyplot as plt
+from numpy import arange
+from openalea import lpy
+from openalea.plantgl.algo import ZBufferEngine
+from openalea.plantgl.algo import eColorBased
+from openalea.plantgl.math import Matrix3
+from openalea.plantgl.math import Vector3
+from openalea.plantgl.math import norm
+from openalea.plantgl.scenegraph import BoundingBox
+from openalea.plantgl.scenegraph import Color3
 
 
 def read_camera(filename):
-    ''' Camera file contains
+    """ Camera file contains
         e.g.:
         104 10 4 -0.583715 -0.430251 35.4873 149.352 -0.430251 35.4873 92 -2.2 -1.8 40.0189 480.226 0.600283 30 79.7913 1 1
         interpreted as:
         azimuth, Elevation, ?4, Center(3-uple), Eye(3-uple), Zoom, Translation(2-uple),?40.0189, zfar, znear, DefaultViewAngle(deg), CurrentViewAngle(deg), ?1, ?1
-    '''
+    """
     with open(filename, 'r') as file:
         data = file.read().replace('\n', '')
     tokens = data.split()
@@ -63,10 +67,10 @@ def read_camera(filename):
         'znear': float(tokens[14]),
         'DefaultViewAngle': (float(tokens[15])),
         'CurrentViewAngle': (float(tokens[16])),
-        'projection_mode': int(tokens[17]),  # Flag othographic, pers
+        'projection_mode': int(tokens[17]),  # Flag orthographic, pers
         'geomsys': int(tokens[18]),  # Flag PlanGL frame (Z up) / OpenGL (Y up)
     }
-    '''
+    """
     target_point = None             # looks a bit above z = 0. Value is a 3-uple (x,y,z)
     zoomcoef = 1.                   # increase to zoom out
     camdist = 1.                    # increase to widen the observation window
@@ -82,15 +86,15 @@ def read_camera(filename):
     ambiant = (255, 255, 255)
     diffuse = (0, 0, 0)
     specular = (0, 0, 0))
-    '''
+    """
 
     return camdict
 
 
 def get_bb(viewer):
-    '''
+    """
     gets the bounding box of a scene contained in a LPy viewer.
-    '''
+    """
     sc = viewer.getCurrentScene()
     return BoundingBox(sc)
 
@@ -98,7 +102,7 @@ def get_bb(viewer):
 # No longer used (prefer offline equivalent function below)
 def take_snapshot(viewer, filename_prefix=".", suffix="",
                   cam_settings={'camdist': 1., 'zoomcoef': 1., 'bb': None, 'target_point': None}):
-    '''
+    """
     take a snapshot of a computed lpy_scene, with camera attributes defined in cam_settings:
     - target_point is the point at which the camera is looking (Vector3)
     - camdist fixes the distance of the camera to the target point (in screen units)
@@ -108,7 +112,7 @@ def take_snapshot(viewer, filename_prefix=".", suffix="",
 
     filename_prefix: A prefix directory may be defined (default is current directory '.')
     suffix: name of the image
-    '''
+    """
 
     if 'camdist' in cam_settings:
         camdist = cam_settings['camdist']
@@ -173,7 +177,7 @@ def take_snapshot(viewer, filename_prefix=".", suffix="",
     # Set beckground Color
     viewer.frameGL.setBgColor(Color3(255, 255, 255))
     # Sets the size of the PlantGL window in pixels
-    viewer.frameGL.setSize((400, 400))
+    viewer.frameGL.setSize(400, 400)
 
     # Defines the new camera position and orientation as
     # shooting direction = c-np
@@ -187,8 +191,8 @@ def take_snapshot(viewer, filename_prefix=".", suffix="",
 
 def take_offline_snapshot(lscene, filename_prefix=".", suffix="",
                           cam_settings={'camdist': 1., 'zoomcoef': 1., 'bb': None, 'target_point': None}):
-    """ take a snapshot of a computed lpy_scene without the need of running LPy.
-
+    """
+    take a snapshot of a computed lpy_scene without the need of running LPy.
     Camera attributes defined in cam_settings:
     - target_point is the point at which the camera is looking (Vector3)
     - camdist fixes the distance of the camera to the target point (in screen units)
@@ -218,7 +222,6 @@ def take_offline_snapshot(lscene, filename_prefix=".", suffix="",
     #  for key, val in d.items():
     #    exec(key + '=val')
 
-
     if 'camdist' in cam_settings:
         camdist = cam_settings['camdist']
     else:
@@ -280,7 +283,7 @@ def take_offline_snapshot(lscene, filename_prefix=".", suffix="",
     else:
         specular = (0, 0, 0)
 
-        # Determine the point to look at
+    # Determine the point to look at
     if target_point == None:
         # define bounding box used to take the picture
         if bb == None:
@@ -293,9 +296,9 @@ def take_offline_snapshot(lscene, filename_prefix=".", suffix="",
     else:
         c = target_point  # target point to look at is given as an argument
 
-        # DISPLAY WINDOW
+    # DISPLAY WINDOW
 
-        # rendering engine (alternative offline to the viewer
+    # rendering engine (alternative offline to the viewer
     zb = ZBufferEngine(width, height, (255, 255, 255), renderingStyle=eColorBased)
     zb.multithreaded = False
 
@@ -341,11 +344,11 @@ def take_offline_snapshot(lscene, filename_prefix=".", suffix="",
 # !!!!!!!!!!! NEW VERSION NOT YET TESTED --> TO BE TESTED
 def take_circular_snapshots(viewer, delta_a=36, filename_prefix=".", suffix="",
                             cam_settings={'camdist': 1., 'zoomcoef': 1., 'bb': None, 'target_point': None}):
-    '''
+    """
     take a snapshot of a computed lpy_scene contained in the viewer.
 
     - delta_a is the increment angle in degrees
-    '''
+    """
 
     if 'camdist' in cam_settings:
         camdist = cam_settings['camdist']
@@ -453,10 +456,10 @@ def simulate_and_shoot(model_filename, variables, suffix, cam_settings):
 
 
 def build_suffix(var, short_names=None):
-    '''
+    """
     takes a dict of variables and constructs a str suffix identifying uniquely this set of parameter
     (for further defining a name corresponding to this set of parameters for storage on disk of corresponding results)
-    '''
+    """
     sep = ' '  # could also be ' ', '_' or '#'
     stg = ''
     index = 0
@@ -478,7 +481,7 @@ def build_suffix(var, short_names=None):
 
 
 def build_variables_dict(x, y, free_variable_list, fixed_variables):
-    '''
+    """
     builds a dictionary of variables by blending free variables (made of two variables)
 
     - x,y are the values of the two free free_variables
@@ -486,7 +489,7 @@ def build_variables_dict(x, y, free_variable_list, fixed_variables):
     - fixed_variables is a dict of variables with already defined values (i.e. fixed)
 
     Precondition: len(free_variable_list) >= 2
-    '''
+    """
     assert len(free_variable_list) >= 2  # (at list names for x and y should be defined)
 
     vardict = fixed_variables.copy()
@@ -499,7 +502,7 @@ def build_variables_dict(x, y, free_variable_list, fixed_variables):
 def grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed_variables,
                             cam_settings={'camdist': 1., 'zoomcoef': 1., 'bb': None, 'target_point': None},
                             short_names=None):
-    '''
+    """
     launches simulations of the lpy model model_filename for free parameters values defined in simpoints,
     - simpoints is a dict whose keys are the y values and the values are lists of x-values for each y key
     - free variable names are given in free_variable_list and
@@ -510,7 +513,7 @@ def grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed
       or a function returning a dict for args = (x,y) = pair of values of the free parameters
         cam_settings(x,y) --> {'camdist':val1, 'zoomcoef':val2, 'bb':val3, 'target_point':val4}
 
-    '''
+    """
 
     # create list of variable names
     fixed_variable_list = fixed_variables.keys()
@@ -560,17 +563,18 @@ def grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed
         index_list.append(xi * XMAX + yi + 1)  # these indexes must be shifted by 1 as they start at 1
     print("array dim = ", YMAX, XMAX)
     print("index_list(", len(index_list), ") = ", index_list)
+
     plot_images(image_name_list, index_list, dir=".", size=(YMAX, XMAX))
 
 
 def plot_images(image_name_list, index_list=None, dir='.', size=(1, 1)):
-    '''
+    """
     Plots the grid of png images @ grid positions defined by index_list
 
     - index_list = 1 dimensional position = flattened 2D position in a grid of size[0]xsize[1]
     - dir specifies the directory where to print
     - size is the dimensions of the underlying 2D grid
-    '''
+    """
 
     # A figure is a canvas that contains axes (= places on which to plot)
     # by default the figure has one axis
@@ -582,12 +586,12 @@ def plot_images(image_name_list, index_list=None, dir='.', size=(1, 1)):
     dim1 = min(size[1] * 3, 15)
     print("Image: grid size= ", size[1], size[0], " image --> width,heigth = ", dim1, dim2)
     fig = plt.figure(figsize=(dim1, dim2))
+    print("Coucou")
 
     i = 1
-
     for iname in image_name_list:
-
-        image = mpimg.imread(dir + '/' + iname)
+        # image = mpimg.imread(dir + '/' + iname)
+        image = mpimg.imread(iname)
 
         if not index_list == None:
             k = index_list[i - 1]
diff --git a/tests/test-snapshots.py b/tests/test-snapshots.py
new file mode 100644
index 0000000000000000000000000000000000000000..85231423b0d58be2278f225210f1e376b5939a17
--- /dev/null
+++ b/tests/test-snapshots.py
@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Take automaticaly pictures of scenes generated by a Lsystem
+
+Authors: Christophe Godin, Inria, 2021
+Licence: Open source LGPL
+
+See README.md in lpy_snapshot module for documentation
+
+Launch with:
+shell> ipython --gui=qt test-snapshots.py
+
+This should display a graphical window with a figure containing 4 subfigures
+"""
+
+import matplotlib
+
+matplotlib.use('tkagg')
+print(f"Using {matplotlib.get_backend()} backend for Matplotlib!")
+
+from openalea.plantgl.all import Vector3
+
+from lpy_tools.snapshots.lscene_snapshots import grid_simulate_and_shoot
+
+model_filename = 'test.lpy'
+
+##############################
+# Input variables
+##############################
+
+Narray = [1, 2, 3, 4]
+
+###########################################
+# Simulate the plant and generate the scans
+###########################################
+
+# def build_suffix(var):
+#     return "O"+str(var["MAX_ORDER"])+"-D"+str(var["NB_DAYS"])+"-MD"+str(var["MD"])+'+' if str(var["SHOW_MERISTEMS"]) else '-'
+
+# Definition of shortnames for variables (used for generating small image names)
+variable_short_names = {}
+
+# Set fixed variables input to the L-system for all the simulations here:
+
+# Set variables x,y that you want to vary in the L-system to build the figure
+# x = [x1,x2,.. ]
+# this is made using a dictionary providing for each md the list of x values to simulate
+# y = dict key,
+# x list = dict values giving for each key the list of x for which a
+# simulation sim(x,y) must be computed
+
+fixed_variables_dict = {}
+free_variable_list = ['N', 'cst']  # Should contain the two variables which will vary (simpoints) to make the grid
+simpoints = {0: Narray}
+
+# Building of the image
+# setting camera options
+target_point = Vector3(0, 0, 0.)  # looks a bit above z = 0
+zoomcoef = 0.5  # increase to zoom out
+camdist = 100.  # increase to widen the observation window
+
+
+# camdict = {'camdist':camdist, 'zoomcoef':zoomcoef, 'bb':None, 'target_point':target_point}
+
+def cam_setting_fun1(x, y):
+    """Defines the camera setting values for each pair of parameter values (x,y).
+
+    Parameters
+    ----------
+    x represents time
+    y represents meristem delay
+    """
+    t = target_point
+    z = zoomcoef
+    c = camdist
+
+    return {'camdist': c, 'zoomcoef': z, 'bb': None, 'target_point': t, 'elevation': 0.0}
+
+
+# Test for reusing a camera file recorded from L-py
+# camdict_save = read_camera('camA')
+# print(camdict_save)
+# grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed_variables_dict, cam_settings = camdict_save) #short_names = variable_short_names)
+
+grid_simulate_and_shoot(simpoints, model_filename, free_variable_list, fixed_variables_dict,
+                        cam_settings=cam_setting_fun1)  # short_names = variable_short_names)
diff --git a/tests/test.lpy b/tests/test.lpy
new file mode 100644
index 0000000000000000000000000000000000000000..243518d89da280aef54d5b2bac1c1c4e73042ba1
--- /dev/null
+++ b/tests/test.lpy
@@ -0,0 +1,14 @@
+extern(N=2)
+extern(cst=0)
+
+Axiom: _(0.1)-(90)G(10)
+
+derivation length: N
+production:
+G(x) --> G(x/3)+G(x/3)--G(x/3)+G(x/3)
+
+interpretation:
+
+G(x) :
+  nproduce F(x)
+endlsystem