diff --git a/src/Map3d.cpp b/src/Map3d.cpp
index ec19b149..f93247b3 100644
--- a/src/Map3d.cpp
+++ b/src/Map3d.cpp
@@ -1191,6 +1191,9 @@ bool Map3d::add_las_file(PointFile pointFile) {
   return true;
 }
 
+/**
+ * query rtrees and iterate results to find adjacent features
+ */
 void Map3d::collect_adjacent_features(TopoFeature* f) {
   std::vector<PairIndexed> re;
   Box2 b = f->get_bbox2d();
diff --git a/src/TopoFeature.cpp b/src/TopoFeature.cpp
index 92ac62fd..a27af900 100644
--- a/src/TopoFeature.cpp
+++ b/src/TopoFeature.cpp
@@ -303,6 +303,11 @@ void TopoFeature::get_citygml_attributes(std::wostream& of, const AttributeMap&
   }
 }
 
+/**
+ * create GDAL feature with attributes and geometry
+ * geometry type is a MultiPolygon
+ * geometry contains polygon and vertical walls as triangles
+ */
 bool TopoFeature::get_multipolygon_features(OGRLayer* layer, std::string className, bool writeAttributes, const AttributeMap& extraAttributes) {
   OGRFeatureDefn *featureDefn = layer->GetLayerDefn();
   OGRFeature *feature = OGRFeature::CreateFeature(featureDefn);
@@ -377,6 +382,9 @@ bool TopoFeature::get_multipolygon_features(OGRLayer* layer, std::string classNa
   return true;
 }
 
+/**
+ * create GDAL attribute in feature based on OGR feature definition
+ */
 bool TopoFeature::writeAttribute(OGRFeature* feature, OGRFeatureDefn* featureDefn, std::string name, std::string value) {
   int fi = featureDefn->GetFieldIndex(name.c_str());
   if (fi == -1) {
@@ -817,6 +825,10 @@ Point2 TopoFeature::get_point2(int ringi, int pi) {
     return _p2->inners()[ringi - 1][pi];
 }
 
+/**
+ * return next vertex in the ring
+ * return first vertex of ring when last vertex is supplied
+ */
 Point2 TopoFeature::get_next_point2_in_ring(int ringi, int i, int& pi) {
   Ring2 ring;
   if (ringi == 0)
@@ -949,6 +961,10 @@ bool TopoFeature::point_in_polygon(const Point2& p) {
   return insideOuter;
 }
 
+/**
+ * cleanup elevation information vectors
+ * clean _lidarelevs and _p2z
+ */
 void TopoFeature::cleanup_lidarelevs() {
   _lidarelevs.clear();
   _lidarelevs.shrink_to_fit();
@@ -1044,6 +1060,10 @@ bool TopoFeature::get_attribute(std::string attributeName, std::string& attribut
   return false;
 }
 
+/**
+ * lift polygon to a single height
+ * used for Flat class
+ */
 void TopoFeature::lift_all_boundary_vertices_same_height(int height) {
   int ringi = 0;
   Ring2& oring = _p2->outer();
@@ -1066,6 +1086,10 @@ std::vector<TopoFeature*>* TopoFeature::get_adjacent_features() {
   return _adjFeatures;
 }
 
+/**
+ * lift each vertex to its calculated height
+ * used for Boundary3D and TIN classes
+ */
 void TopoFeature::lift_each_boundary_vertices(float percentile) {
   //-- assign value for each vertex based on percentile
   bool hasHeight = false;
diff --git a/src/geomtools.cpp b/src/geomtools.cpp
index 83ee884b..2dc29a20 100644
--- a/src/geomtools.cpp
+++ b/src/geomtools.cpp
@@ -125,12 +125,15 @@ void mark_domains(CDT& ct,
   }
 }
 
-//explore set of facets connected with non constrained edges,
-//and attribute to each such set a nesting level.
-//We start from facets incident to the infinite vertex, with a nesting
-//level of 0. Then we recursively consider the non-explored facets incident 
-//to constrained edges bounding the former set and increase the nesting level by 1.
-//Facets in the domain are those with an odd nesting level.
+/**
+ * mark the triangles that are inside the original 2D polygon.
+ * explore set of facets connected with non constrained edges,
+ * and attribute to each such set a nesting level.
+ * start from facets incident to the infinite vertex, with a nesting
+ * level of 0. Then recursively consider the non-explored facets incident
+ * to constrained edges bounding the former set and increase the nesting level by 1.
+ * facets in the domain are those with an odd nesting level.
+ */
 void mark_domains(CDT& cdt) {
   for (CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it) {
     it->info().nesting_level = -1;
@@ -147,6 +150,11 @@ void mark_domains(CDT& cdt) {
   }
 }
 
+/**
+ * triangulation the polygon with lifted vertices
+ * add height points to the triangulation in case of Terrain and Forest class
+ * push created vertices and triangles to the output vectors
+ */
 bool getCDT(Polygon2* pgn,
   const std::vector< std::vector<int> > &z,
   std::vector< std::pair<Point3, std::string> > &vertices,