Index: cpp/frams/genetics/fS/fS_conv.h
===================================================================
--- cpp/frams/genetics/fS/fS_conv.h	(revision 1017)
+++ cpp/frams/genetics/fS/fS_conv.h	(revision 1030)
@@ -17,5 +17,5 @@
 	GenoConv_fS0s() : GenoConverter()
 	{
-		name = "Solid encoding";
+		name = "Solids tree-structure encoding";
 
 		in_format = "S";
Index: cpp/frams/genetics/fS/fS_general.cpp
===================================================================
--- cpp/frams/genetics/fS/fS_general.cpp	(revision 1017)
+++ cpp/frams/genetics/fS/fS_general.cpp	(revision 1030)
@@ -15,5 +15,4 @@
 
 int fS_Genotype::precision = 4;
-bool fS_Genotype::TURN_WITH_ROTATION = false;
 std::map<string, double> Node::minValues;
 std::map<string, double> Node::defaultValues;
@@ -33,8 +32,8 @@
 				{RY,        -M_PI},
 				{RZ,        -M_PI},
-				{SIZE,      0.01},
-				{SIZE_X,    Model::getMinPart().scale.x},
-				{SIZE_Y,    Model::getMinPart().scale.y},
-				{SIZE_Z,    Model::getMinPart().scale.z}
+				{SCALE,      0.01},
+				{SCALE_X,    Model::getMinPart().scale.x},
+				{SCALE_Y,    Model::getMinPart().scale.y},
+				{SCALE_Z,    Model::getMinPart().scale.z}
 		};
 	}
@@ -51,8 +50,8 @@
 				{RY,        M_PI},
 				{RZ,        M_PI},
-				{SIZE,      100.0},
-				{SIZE_X,    Model::getMaxPart().scale.x},
-				{SIZE_Y,    Model::getMaxPart().scale.y},
-				{SIZE_Z,    Model::getMaxPart().scale.z}
+				{SCALE,      100.0},
+				{SCALE_X,    Model::getMaxPart().scale.x},
+				{SCALE_Y,    Model::getMaxPart().scale.y},
+				{SCALE_Z,    Model::getMaxPart().scale.z}
 		};
 	}
@@ -68,8 +67,8 @@
 				{RY,        0.0},
 				{RZ,        0.0},
-				{SIZE,      1.0},
-				{SIZE_X,    Model::getDefPart().scale.x},
-				{SIZE_Y,    Model::getDefPart().scale.y},
-				{SIZE_Z,    Model::getDefPart().scale.z}
+				{SCALE,      1.0},
+				{SCALE_X,    Model::getDefPart().scale.x},
+				{SCALE_Y,    Model::getDefPart().scale.y},
+				{SCALE_Z,    Model::getDefPart().scale.z}
 		};
 	}
@@ -111,7 +110,14 @@
 }
 
+void rotateVector(Pt3D &vector, const Pt3D &rotation)
+{
+	Orient rotmatrix = Orient_1;
+	rotmatrix.rotate(rotation);
+	vector = rotmatrix.transform(vector);
+}
+
 void State::rotate(const Pt3D &rotation)
 {
-       fS_Utils::rotateVector(v, rotation);
+       rotateVector(v, rotation);
        v.normalize();
 }
@@ -218,9 +224,9 @@
 void Node::extractModifiers(Substring &restOfGenotype)
 {
-	int partTypePosition = getPartPosition(restOfGenotype);
-	if (partTypePosition == -1)
+	int partShapePosition = getPartPosition(restOfGenotype);
+	if (partShapePosition == -1)
 		throw fS_Exception("Part type missing", restOfGenotype.start);
 
-	for (int i = 0; i < partTypePosition; i++)
+	for (int i = 0; i < partShapePosition; i++)
 	{
 		// Extract modifiers and joint
@@ -233,5 +239,5 @@
 			throw fS_Exception("Invalid modifier", restOfGenotype.start + i);
 	}
-	restOfGenotype.startFrom(partTypePosition);
+	restOfGenotype.startFrom(partShapePosition);
 }
 
@@ -242,5 +248,5 @@
 		throw fS_Exception("Invalid part type", restOfGenotype.start);
 
-	partType = itr->second;
+	partShape = itr->second;
 	restOfGenotype.startFrom(1);
 }
@@ -327,5 +333,5 @@
 		size_t len = length - valueStartIndex;
 		double value = fS_stod(val, restOfGenotype.start + start + valueStartIndex, &len);
-		if((key==SIZE_X || key==SIZE_Y || key==SIZE_Z) && value <= 0.0)
+		if((key==SCALE_X || key==SCALE_Y || key==SCALE_Z) && value <= 0.0)
 			throw fS_Exception("Invalid value of radius parameter", restOfGenotype.start + start + valueStartIndex);
 
@@ -342,5 +348,9 @@
 	if (item != params.end())
 		return item->second;
-	return defaultValues.at(key);
+
+	auto defaultItem = defaultValues.find(key);
+	if(defaultItem == defaultValues.end())
+		throw fS_Exception("Default value missing", 0);
+	return defaultItem->second;
 }
 
@@ -382,5 +392,5 @@
 		state->rotate(getVectorRotation());
 
-		double distance = getDistance();
+		double distance = calculateDistanceFromParent();
 		state->addVector(distance);
 	}
@@ -433,10 +443,10 @@
 }
 
-void Node::calculateSize(Pt3D &scale)
-{
-	double sizeMultiplier = getParam(SIZE) * state->s;
-	scale.x = getParam(SIZE_X) * sizeMultiplier;
-	scale.y = getParam(SIZE_Y) * sizeMultiplier;
-	scale.z = getParam(SIZE_Z) * sizeMultiplier;
+void Node::calculateScale(Pt3D &scale)
+{
+	double scaleMultiplier = getParam(SCALE) * state->s;
+	scale.x = getParam(SCALE_X) * scaleMultiplier;
+	scale.y = getParam(SCALE_Y) * scaleMultiplier;
+	scale.z = getParam(SCALE_Z) * scaleMultiplier;
 }
 
@@ -444,8 +454,8 @@
 {
 	double result;
-	Pt3D size;
-	calculateSize(size);
-	double radiiProduct = size.x * size.y * size.z;
-	switch (partType)
+	Pt3D scale;
+	calculateScale(scale);
+	double radiiProduct = scale.x * scale.y * scale.z;
+	switch (partShape)
 	{
 		case Part::Shape::SHAPE_CUBOID:
@@ -464,8 +474,8 @@
 }
 
-bool Node::isPartSizeValid()
-{
-	Pt3D size;
-	calculateSize(size);
+bool Node::isPartScaleValid()
+{
+	Pt3D scale;
+	calculateScale(scale);
 	double volume = calculateVolume();
 	Part_MinMaxDef minP = Model::getMinPart();
@@ -474,13 +484,13 @@
 	if (volume > maxP.volume || minP.volume > volume)
 		return false;
-	if (size.x < minP.scale.x || size.y < minP.scale.y || size.z < minP.scale.z)
+	if (scale.x < minP.scale.x || scale.y < minP.scale.y || scale.z < minP.scale.z)
 		return false;
-	if (size.x > maxP.scale.x || size.y > maxP.scale.y || size.z > maxP.scale.z)
+	if (scale.x > maxP.scale.x || scale.y > maxP.scale.y || scale.z > maxP.scale.z)
 		return false;
 
-	if (partType == Part::Shape::SHAPE_ELLIPSOID && fS_Utils::max3(size) != fS_Utils::min3(size))
+	if (partShape == Part::Shape::SHAPE_ELLIPSOID && scale.maxComponentValue() != scale.minComponentValue())
 		// When not all radii have different values
 		return false;
-	if (partType == Part::Shape::SHAPE_CYLINDER && size.y != size.z)
+	if (partShape == Part::Shape::SHAPE_CYLINDER && scale.y != scale.z)
 		// If base radii have different values
 		return false;
@@ -488,9 +498,4 @@
 }
 
-bool Node::hasPartSizeParam()
-{
-	return params.count(SIZE_X) > 0 || params.count(SIZE_Y) > 0 || params.count(SIZE_Z) > 0;
-}
-
 Pt3D Node::getVectorRotation()
 {
@@ -501,5 +506,5 @@
 {
 	Pt3D rotation = Pt3D(getParam(RX, 0.0), getParam(RY, 0.0), getParam(RZ, 0.0));
-	if(fS_Genotype::TURN_WITH_ROTATION)
+	if(genotypeParams.turnWithRotation)
 		rotation += getVectorRotation();
 	return rotation;
@@ -537,10 +542,10 @@
 void Node::createPart()
 {
-	part = new Part(partType);
+	part = new Part(partShape);
 	part->p = Pt3D(state->location);
 
 	part->friction = getParam(FRICTION) * state->fr;
 	part->ingest = getParam(INGESTION) * state->ing;
-	calculateSize(part->scale);
+	calculateScale(part->scale);
 	part->setRot(getRotation());
 }
@@ -581,5 +586,5 @@
 		result += std::string(count, mod).c_str();
 	}
-	result += SHAPE_TO_GENE.at(partType);
+	result += SHAPE_TO_GENE.at(partShape);
 
 	if (!neurons.empty())
@@ -658,10 +663,17 @@
 }
 
-fS_Genotype::fS_Genotype(const string &geno)
-{
+fS_Genotype::fS_Genotype(const string &g)
+{
+	string geno(g);
+	geno.erase(remove(geno.begin(), geno.end(), ' '), geno.end());
+	geno.erase(remove(geno.begin(), geno.end(), '\n'), geno.end());
 	try
 	{
 		GenotypeParams genotypeParams;
 		genotypeParams.modifierMultiplier = 1.1;
+		genotypeParams.distanceTolerance = 0.1;
+		genotypeParams.relativeDensity = 10.0;
+		genotypeParams.turnWithRotation = false;
+		genotypeParams.paramMutationStrength = 0.4;
 
 		size_t modeSeparatorIndex = geno.find(MODE_SEPARATOR);
@@ -669,5 +681,21 @@
 			throw fS_Exception("Genotype parameters missing", 0);
 
-		genotypeParams.modifierMultiplier = fS_stod(geno, 0, &modeSeparatorIndex);
+		std::vector<SString> paramStrings;
+		strSplit(SString(geno.c_str(), modeSeparatorIndex), ',', false, paramStrings);
+
+		if(paramStrings.size() >= 1 && paramStrings[0] != "")
+		{
+			size_t len0 = paramStrings[0].length();
+			genotypeParams.modifierMultiplier = fS_stod(paramStrings[0].c_str(), 0, &len0);
+		}
+		if(paramStrings.size() >= 2 && paramStrings[1] != "")
+		{
+			genotypeParams.turnWithRotation = bool(atoi(paramStrings[1].c_str()));
+		}
+		if(paramStrings.size() >= 3 && paramStrings[2] != "")
+		{
+			size_t len2 = paramStrings[2].length();
+			genotypeParams.paramMutationStrength = fS_stod(paramStrings[2].c_str(), 0, &len2);
+		}
 
 		int genoStart = modeSeparatorIndex + 1;
@@ -680,4 +708,9 @@
 		delete startNode;
 		throw e;
+	}
+	catch(...)
+	{
+		delete startNode;
+		throw fS_Exception("Unknown exception in fS", 0);
 	}
 }
@@ -753,5 +786,10 @@
 	geno.reserve(100);
 
-	geno += doubleToString(startNode->genotypeParams.modifierMultiplier, fS_Genotype::precision).c_str();
+	GenotypeParams gp = startNode->genotypeParams;
+	geno += doubleToString(gp.modifierMultiplier, precision).c_str();
+	geno += ",";
+	geno += doubleToString(gp.turnWithRotation, precision).c_str();
+	geno += ",";
+	geno += doubleToString(gp.paramMutationStrength, precision).c_str();
 	geno += MODE_SEPARATOR;
 
@@ -850,5 +888,5 @@
 	for (int i = 0; i < int(nodes.size()); i++)
 	{
-		if (!nodes[i]->isPartSizeValid())
+		if (!nodes[i]->isPartScaleValid())
 		{
 			return 1 + nodes[i]->partDescription->start;
@@ -884,2 +922,26 @@
 }
 
+double Node::calculateDistanceFromParent()
+{
+	Pt3D scale;
+	calculateScale(scale);
+	Pt3D parentScale;
+	parent->calculateScale(parentScale);    // Here we are sure that parent is not nullptr
+	Part *tmpPart = PartDistanceEstimator::buildTemporaryPart(partShape, scale, getRotation());
+	Part *parentTmpPart = PartDistanceEstimator::buildTemporaryPart(parent->partShape, parentScale, parent->getRotation());
+
+	double result;
+	try
+	{
+		tmpPart->p = state->v;
+		result = PartDistanceEstimator::calculateDistance(*tmpPart, *parentTmpPart, genotypeParams.distanceTolerance, genotypeParams.relativeDensity);
+	}
+	catch (...)
+	{
+		throw fS_Exception("Exception thrown while calculating distance from parent", 0);
+	}
+
+	delete tmpPart;
+	delete parentTmpPart;
+	return result;
+}
Index: cpp/frams/genetics/fS/fS_general.h
===================================================================
--- cpp/frams/genetics/fS/fS_general.h	(revision 1017)
+++ cpp/frams/genetics/fS/fS_general.h	(revision 1030)
@@ -43,8 +43,8 @@
 #define FRICTION "f"
 #define STIFFNESS "st"
-#define SIZE "s"
-#define SIZE_X "x"
-#define SIZE_Y "y"
-#define SIZE_Z "z"
+#define SCALE "s"
+#define SCALE_X "x"
+#define SCALE_Y "y"
+#define SCALE_Z "z"
 #define ROT_X "tx"
 #define ROT_Y "ty"
@@ -83,7 +83,7 @@
 const int JOINT_COUNT = JOINTS.length();
 const string MODIFIERS = "IFS";
-const char SIZE_MODIFIER = 's';
-const vector<string> PARAMS {INGESTION, FRICTION, ROT_X, ROT_Y, ROT_Z, RX, RY, RZ, SIZE, SIZE_X, SIZE_Y, SIZE_Z};
-const vector<string> SIZE_PARAMS {SIZE, SIZE_X, SIZE_Y, SIZE_Z};
+const char SCALE_MODIFIER = 's';
+const vector<string> PARAMS {INGESTION, FRICTION, ROT_X, ROT_Y, ROT_Z, RX, RY, RZ, SCALE, SCALE_X, SCALE_Y, SCALE_Z};
+const vector<string> SCALE_PARAMS {SCALE, SCALE_X, SCALE_Y, SCALE_Z};
 
 /** @name Default values of node parameters*/
@@ -247,4 +247,13 @@
 struct GenotypeParams{
 	double modifierMultiplier;	// Every modifier changes the underlying value by this multiplier
+	/// When calculating the distance between parts, the internal result is a range of numbers
+	/// distanceTolerance is the maximal allowed size of this range
+	double distanceTolerance;
+	/// Used for deriving density for MeshBuilder
+	double relativeDensity;
+	///
+	bool turnWithRotation;
+	///
+	double paramMutationStrength;
 };
 
@@ -265,5 +274,7 @@
 	Part *part;     /// A part object built from node. Used in building the Model
 	int partCodeLen; /// The length of substring that directly describes the corresponding part
-	GenotypeParams genotypeParams;
+	static std::map<string, double> minValues;	/// Min parameter values
+	static std::map<string, double> defaultValues;	/// Default parameter values
+	static std::map<string, double> maxValues;	/// Max parameter values
 
 	vector<Node *> children;    /// Vector of all direct children
@@ -273,15 +284,13 @@
 	void prepareParams();
 
-	double getDistance();
-
 	void cleanUp();
 
+	/// Get part's  rotation
 	Pt3D getRotation();
 
+	/// Get the rotation of vector
 	Pt3D getVectorRotation();
 
-	bool isPartSizeValid();
-
-	bool hasPartSizeParam();
+	bool isPartScaleValid();
 
 	/**
@@ -362,11 +371,9 @@
 
 public:
-	static std::map<string, double> minValues;
-	static std::map<string, double> defaultValues;
-	static std::map<string, double> maxValues;
 	char joint = DEFAULT_JOINT;           /// Set of all joints
-	Part::Shape partType;  /// The type of the part
+	Part::Shape partShape;  /// The type of the part
 	State *state = nullptr; /// The phenotypic state that inherits from ancestors
 	std::map<string, double> params; /// The map of all the node params
+	GenotypeParams genotypeParams; /// Parameters that affect the whole genotype
 
 	Node(Substring &genotype, Node *parent, GenotypeParams genotypeParams);
@@ -381,8 +388,8 @@
 
 	/**
-	 * Calculate the effective size of the part (after applying all multipliers and params)
-	 * @return The effective size
-	 */
-	void calculateSize(Pt3D &scale);
+	 * Calculate the effective scale of the part (after applying all multipliers and params)
+	 * @return The effective scales
+	 */
+	void calculateScale(Pt3D &scale);
 
 	/**
@@ -404,4 +411,7 @@
 	double getParam(const string &key);
 	double getParam(const string &key, double defaultValue);
+
+	/// Calculate distance between the part its parent
+	double calculateDistanceFromParent();
 };
 
@@ -439,6 +449,5 @@
 
 
-	static int precision;
-	static bool TURN_WITH_ROTATION;
+	static int precision; /// Number of decimal places for numbers in genotype
 
 	/**
@@ -450,4 +459,5 @@
 	~fS_Genotype();
 
+	/// Calculate the State field for all the nodes
 	void getState(bool calculateLocation);
 
@@ -495,5 +505,5 @@
 
 	/**
-	 * Check if sizes of all parts in genotype are valid
+	 * Check if scales of all parts in genotype are valid
 	\retval error_position 1-based
 	\retval 0 when all part sizes are valid
Index: cpp/frams/genetics/fS/fS_oper.cpp
===================================================================
--- cpp/frams/genetics/fS/fS_oper.cpp	(revision 1017)
+++ cpp/frams/genetics/fS/fS_oper.cpp	(revision 1030)
@@ -51,5 +51,5 @@
 		if (errorPosition != 0)
 		{
-			logPrintf("GenoOper_fS", "checkValidity", LOG_WARN, "Invalid part size");
+			logPrintf("GenoOper_fS", "checkValidity", LOG_WARN, "Invalid part scale");
 			return errorPosition;
 		}
@@ -226,5 +226,5 @@
 const char *GenoOper_fS::getSimplest()
 {
-	return "1.1:C{x=0.80599;y=0.80599;z=0.80599}";
+	return "1.1,0,0.4:C{x=0.80599;y=0.80599;z=0.80599}";
 }
 
@@ -297,7 +297,7 @@
 	geno.getState(false);
 	Node *node = geno.chooseNode();
-	char partType = SHAPE_TO_GENE.at(availablePartShapes[rndUint(availablePartShapes.size())]);
-
-	Substring substring(&partType, 0, 1);
+	char partShape = SHAPE_TO_GENE.at(availablePartShapes[rndUint(availablePartShapes.size())]);
+
+	Substring substring(&partShape, 0, 1);
 	Node *newNode = new Node(substring, node, node->genotypeParams);
 	// Add random rotation
@@ -322,6 +322,6 @@
 		newNode->params[selectedParam] = RndGen.Uni(-M_PI / 2, M_PI / 2);
 	}
-	// Assign part size to default value
-	double volumeMultiplier = pow(node->getParam(SIZE) * node->state->s, 3);
+	// Assign part scale to default value
+	double volumeMultiplier = pow(node->getParam(SCALE) * node->state->s, 3);
 	double minVolume = Model::getMinPart().volume;
 	double defVolume = Model::getDefPart().volume * volumeMultiplier;    // Default value after applying modifiers
@@ -330,8 +330,8 @@
 	double relativeVolume = volume / volumeMultiplier;    // Volume without applying modifiers
 
-	double newRadius = std::cbrt(relativeVolume / volumeMultipliers.at(newNode->partType));
-	newNode->params[SIZE_X] = newRadius;
-	newNode->params[SIZE_Y] = newRadius;
-	newNode->params[SIZE_Z] = newRadius;
+	double newRadius = std::cbrt(relativeVolume / volumeMultipliers.at(newNode->partShape));
+	newNode->params[SCALE_X] = newRadius;
+	newNode->params[SCALE_Y] = newRadius;
+	newNode->params[SCALE_Z] = newRadius;
 	node->children.push_back(newNode);
 
@@ -339,7 +339,7 @@
 	{
 		geno.getState(false);
-		mutateSizeParam(newNode, SIZE_X, true);
-		mutateSizeParam(newNode, SIZE_Y, true);
-		mutateSizeParam(newNode, SIZE_Z, true);
+		mutateScaleParam(newNode, SCALE_X, true);
+		mutateScaleParam(newNode, SCALE_Y, true);
+		mutateScaleParam(newNode, SCALE_Z, true);
 	}
 	return true;
@@ -350,5 +350,5 @@
 	Node *randomNode, *selectedChild;
 	// Choose a parent with children
-	// It may be difficult to choose a eligible node, so the number of tries should be high
+	// It may be difficult to choose an eligible node, so the number of tries should be high
 	for (int i = 0; i < 10 * mutationTries; i++)
 	{
@@ -380,40 +380,40 @@
 		Node *randomNode = geno.chooseNode();
 		int index = rndUint(availShapesLen);
-		if (availablePartShapes[index] == randomNode->partType)
+		if (availablePartShapes[index] == randomNode->partShape)
 			index = (index + 1 + rndUint(availShapesLen - 1)) % availShapesLen;
 		Part::Shape newType = availablePartShapes[index];
 
 #ifdef _DEBUG
-		if(newType == randomNode->partType)
+		if(newType == randomNode->partShape)
 			throw fS_Exception("Internal error: invalid part type chosen in mutation.", 1);
 #endif
 
 		geno.getState(false);
-		double sizeMultiplier = randomNode->getParam(SIZE) * randomNode->state->s;
-		double relativeVolume = randomNode->calculateVolume() / pow(sizeMultiplier, 3.0);
-
-		if (!ensureCircleSection || newType == Part::Shape::SHAPE_CUBOID || (randomNode->partType == Part::Shape::SHAPE_ELLIPSOID && newType == Part::Shape::SHAPE_CYLINDER))
-		{
-			double radiusQuotient = std::cbrt(volumeMultipliers.at(randomNode->partType) / volumeMultipliers.at(newType));
-			randomNode->params[SIZE_X] = randomNode->getParam(SIZE_X) * radiusQuotient;
-			randomNode->params[SIZE_Y] = randomNode->getParam(SIZE_Y) * radiusQuotient;
-			randomNode->params[SIZE_Z] = randomNode->getParam(SIZE_Z) * radiusQuotient;
-		} else if (randomNode->partType == Part::Shape::SHAPE_CUBOID && newType == Part::Shape::SHAPE_CYLINDER)
-		{
-			double newRadius = 0.5 * (randomNode->getParam(SIZE_X) + randomNode->getParam(SIZE_Y));
-			randomNode->params[SIZE_X] = 0.5 * relativeVolume / (M_PI * newRadius * newRadius);
-			randomNode->params[SIZE_Y] = newRadius;
-			randomNode->params[SIZE_Z] = newRadius;
+		double scaleMultiplier = randomNode->getParam(SCALE) * randomNode->state->s;
+		double relativeVolume = randomNode->calculateVolume() / pow(scaleMultiplier, 3.0);
+
+		if (!ensureCircleSection || newType == Part::Shape::SHAPE_CUBOID || (randomNode->partShape == Part::Shape::SHAPE_ELLIPSOID && newType == Part::Shape::SHAPE_CYLINDER))
+		{
+			double radiusQuotient = std::cbrt(volumeMultipliers.at(randomNode->partShape) / volumeMultipliers.at(newType));
+			randomNode->params[SCALE_X] = randomNode->getParam(SCALE_X) * radiusQuotient;
+			randomNode->params[SCALE_Y] = randomNode->getParam(SCALE_Y) * radiusQuotient;
+			randomNode->params[SCALE_Z] = randomNode->getParam(SCALE_Z) * radiusQuotient;
+		} else if (randomNode->partShape == Part::Shape::SHAPE_CUBOID && newType == Part::Shape::SHAPE_CYLINDER)
+		{
+			double newRadius = 0.5 * (randomNode->getParam(SCALE_X) + randomNode->getParam(SCALE_Y));
+			randomNode->params[SCALE_X] = 0.5 * relativeVolume / (M_PI * newRadius * newRadius);
+			randomNode->params[SCALE_Y] = newRadius;
+			randomNode->params[SCALE_Z] = newRadius;
 		} else if (newType == Part::Shape::SHAPE_ELLIPSOID)
 		{
 			double newRelativeRadius = cbrt(relativeVolume / volumeMultipliers.at(newType));
-			randomNode->params[SIZE_X] = newRelativeRadius;
-			randomNode->params[SIZE_Y] = newRelativeRadius;
-			randomNode->params[SIZE_Z] = newRelativeRadius;
+			randomNode->params[SCALE_X] = newRelativeRadius;
+			randomNode->params[SCALE_Y] = newRelativeRadius;
+			randomNode->params[SCALE_Z] = newRelativeRadius;
 		} else
 		{
 			throw fS_Exception("Invalid part type", 1);
 		}
-		randomNode->partType = newType;
+		randomNode->partShape = newType;
 		return true;
 	}
@@ -446,11 +446,11 @@
 		return false;
 	// Do not allow invalid changes in part size
-	bool isRadiusOfBase = key == SIZE_Y || key == SIZE_Z;
-	bool isRadius = isRadiusOfBase || key == SIZE_X;
+	bool isRadiusOfBase = key == SCALE_Y || key == SCALE_Z;
+	bool isRadius = isRadiusOfBase || key == SCALE_X;
 	if (ensureCircleSection && isRadius)
 	{
-		if (randomNode->partType == Part::Shape::SHAPE_ELLIPSOID)
+		if (randomNode->partShape == Part::Shape::SHAPE_ELLIPSOID)
 			return false;
-		if (randomNode->partType == Part::Shape::SHAPE_CYLINDER && isRadiusOfBase)
+		if (randomNode->partShape == Part::Shape::SHAPE_CYLINDER && isRadiusOfBase)
 			return false;
 	}
@@ -490,11 +490,14 @@
 bool GenoOper_fS::mutateParamValue(Node *node, string key)
 {
-	// Do not allow invalid changes in part size
-	if (std::find(SIZE_PARAMS.begin(), SIZE_PARAMS.end(), key) == SIZE_PARAMS.end())
-	{
-		node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), Node::minValues.at(key), Node::maxValues.at(key), true);
+	// Do not allow invalid changes in part scale
+	if (std::find(SCALE_PARAMS.begin(), SCALE_PARAMS.end(), key) == SCALE_PARAMS.end())
+	{
+		double max = Node::maxValues.at(key);
+		double min = Node::minValues.at(key);
+		double stddev = (max - min) * node->genotypeParams.paramMutationStrength;
+		node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), min, max, stddev, true);
 		return true;
 	} else
-		return mutateSizeParam(node, key, ensureCircleSection);
+		return mutateScaleParam(node, key, ensureCircleSection);
 }
 
@@ -524,5 +527,5 @@
 	randomNode->modifiers[randomModifier] += rndUint(2) == 0 ? 1 : -1;
 
-	bool isSizeMod = tolower(randomModifier) == SIZE_MODIFIER;
+	bool isSizeMod = tolower(randomModifier) == SCALE_MODIFIER;
 	if (isSizeMod && geno.checkValidityOfPartSizes() != 0)
 	{
@@ -678,10 +681,10 @@
 }
 
-bool GenoOper_fS::mutateSizeParam(Node *node, string key, bool ensureCircleSection)
+bool GenoOper_fS::mutateScaleParam(Node *node, string key, bool ensureCircleSection)
 {
 	double oldValue = node->getParam(key);
 	double volume = node->calculateVolume();
 	double valueAtMinVolume, valueAtMaxVolume;
-	if(key == SIZE)
+	if(key == SCALE)
 	{
 		valueAtMinVolume = oldValue * std::cbrt(Model::getMinPart().volume / volume);
@@ -696,8 +699,9 @@
 	double min = std::max(Node::minValues.at(key), valueAtMinVolume);
 	double max = std::min(Node::maxValues.at(key), valueAtMaxVolume);
-
-	node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), min, max, true);
-
-	if (!ensureCircleSection || node->isPartSizeValid())
+	double stdev = (max - min) * node->genotypeParams.paramMutationStrength;
+
+	node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), min, max, stdev, true);
+
+	if (!ensureCircleSection || node->isPartScaleValid())
 		return true;
 	else
Index: cpp/frams/genetics/fS/fS_oper.h
===================================================================
--- cpp/frams/genetics/fS/fS_oper.h	(revision 1017)
+++ cpp/frams/genetics/fS/fS_oper.h	(revision 1030)
@@ -142,6 +142,6 @@
 
 	/**
-	 * Change the value of the size parameter by given multiplier
-	 * Do not change the value if any of the size restrictions is not satisfied
+	 * Change the value of the scale parameter by given multiplier
+	 * Do not change the value if any of the scale restrictions is not satisfied
 	 * @param paramKey
 	 * @param multiplier
@@ -149,5 +149,5 @@
 	 * @return True if the parameter value was change, false otherwise
 	 */
-	bool mutateSizeParam(Node *node, string key, bool ensureCircleSection);
+	bool mutateScaleParam(Node *node, string key, bool ensureCircleSection);
 };
 
Index: cpp/frams/genetics/fS/part_distance_estimator.h
===================================================================
--- cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1017)
+++ cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1030)
@@ -8,120 +8,84 @@
 #include "frams/model/geometry/meshbuilder.h"
 
-class fS_Utils
+class PartDistanceEstimator
 {
 public:
-	static void rotateVector(Pt3D &vector, const Pt3D &rotation)
+
+	static Part *buildTemporaryPart(Part::Shape shape, const Pt3D &scale, const Pt3D &rotation)
 	{
-		Orient rotmatrix = Orient_1;
-		rotmatrix.rotate(rotation);
-		vector = rotmatrix.transform(vector);
+		Part *tmpPart1 = new Part(shape);
+		tmpPart1->scale = scale;
+		tmpPart1->setRot(rotation);
+		return tmpPart1;
 	}
 
-	static double avg(double a, double b)
-	{
-		return 0.5 * (a + b);
-	}
-
-	static double min3(const Pt3D &p)
-	{
-		double tmp = p.x;
-		if (p.y < tmp)
-			tmp = p.y;
-		if (p.z < tmp)
-			tmp = p.z;
-		return tmp;
-	}
-
-	static double max3(const Pt3D &p)
-	{
-		double tmp = p.x;
-		if (p.y > tmp)
-			tmp = p.y;
-		if (p.z > tmp)
-			tmp = p.z;
-		return tmp;
-	}
-};
-
-class PartDistanceEstimator
-{
-
-public:
-	static constexpr double PRECISION = 0.05;
-	static constexpr double RELATIVE_DENSITY = 5.0;
-
-
-	static Part *buildTemporaryPart(Part::Shape shape, const Pt3D &scale, const Pt3D &rotations)
-	{
-		Part *tmpPart = new Part(shape);
-		tmpPart->scale = scale;
-		tmpPart->setRot(rotations);
-		return tmpPart;
-	}
-
-	static vector <Pt3D> findSphereCenters(Part *part)
+	/// Get some of the points from the surface of the part
+	static vector <Pt3D> findSurfacePoints(Part *part, double  relativeDensity)
 	{
 		// Divide by maximal radius to avoid long computations
-		MeshBuilder::PartSurface surface(RELATIVE_DENSITY / fS_Utils::max3(part->scale));
+		MeshBuilder::PartSurface surface(relativeDensity / part->scale.maxComponentValue());
 		surface.initialize(part);
 
-		vector <Pt3D> centers;
+		vector <Pt3D> points;
 		Pt3D point;
 		while (surface.tryGetNext(point))
 		{
-			centers.push_back(point);
+			points.push_back(point);
 		}
-		return centers;
+		return points;
 	}
 
-	static bool isCollision(Part *parentPart, vector <Pt3D> &centers, Pt3D &vectorBetweenParts)
+	/// Check if there is a collision between the parts
+	static bool isCollision(Part *part, vector <Pt3D> &points, Pt3D &vectorBetweenParts)
 	{
 		static double CBRT_3 = std::cbrt(3);
-		double maxParentReachSq = pow(CBRT_3 * fS_Utils::max3(parentPart->scale), 2);
-		for (int i = 0; i < int(centers.size()); i++)
+		double maxPartReachSq = pow(CBRT_3 * part->scale.maxComponentValue(), 2);
+		for (int i = 0; i < int(points.size()); i++)
 		{
-			Pt3D shifted = centers[i] + vectorBetweenParts;
-			double distanceToCenterSq = shifted.x * shifted.x + shifted.y * shifted.y + shifted.z * shifted.z;
-			if (distanceToCenterSq <= maxParentReachSq && GeometryUtils::isPointInsidePart(shifted, parentPart))
+			Pt3D shifted = points[i] + vectorBetweenParts;
+			double distanceToPointSq = shifted.x * shifted.x + shifted.y * shifted.y + shifted.z * shifted.z;
+			if (distanceToPointSq <= maxPartReachSq && GeometryUtils::isPointInsidePart(shifted, part))
 				return true;
 		}
 		return false;
 	}
+
+
+	static double calculateDistance(Part tmpPart1, Part tmpPart2, double distanceTolerance, double relativeDensity)
+	{
+		/// tmpPart1 and tmpPart2 are copied for purpose and should not be passed as reference
+		/// This function can change some of the properties of those parts
+		Pt3D directionVersor = tmpPart1.p - tmpPart2.p;
+		directionVersor.normalize();
+
+		tmpPart1.p = Pt3D(0);
+		tmpPart2.p = Pt3D(0);
+
+		static double CBRT_3 = std::cbrt(3);
+		vector <Pt3D> points = PartDistanceEstimator::findSurfacePoints(&tmpPart1, relativeDensity);
+
+		double minDistance = tmpPart2.scale.minComponentValue() + tmpPart1.scale.minComponentValue();
+		double maxDistance = CBRT_3 * (tmpPart2.scale.maxComponentValue() + tmpPart1.scale.maxComponentValue());
+		double currentDistance = 0.5 * (maxDistance + minDistance);
+		int collisionDetected = false;
+		while (maxDistance - minDistance > distanceTolerance)
+		{
+			Pt3D vectorBetweenParts = directionVersor * currentDistance;
+			collisionDetected = PartDistanceEstimator::isCollision(&tmpPart2, points, vectorBetweenParts);
+
+			if (collisionDetected)
+			{
+				minDistance = currentDistance;
+				currentDistance = 0.5 * (maxDistance + currentDistance);
+			} else
+			{
+				maxDistance = currentDistance;
+				currentDistance = 0.5 * (currentDistance + minDistance);
+			}
+		}
+		return currentDistance;
+	}
 };
 
-double Node::getDistance()
-{
-	Pt3D size;
-	calculateSize(size);
-	Pt3D parentSize;
-	parent->calculateSize(parentSize);    // Here we are sure that parent is not nullptr
-	Part *tmpPart = PartDistanceEstimator::buildTemporaryPart(partType, size, getRotation());
-	Part *parentTmpPart = PartDistanceEstimator::buildTemporaryPart(parent->partType, parentSize, parent->getRotation());
-
-	vector <Pt3D> centers = PartDistanceEstimator::findSphereCenters(tmpPart);
-
-	double minDistance = 0.0;
-	double maxDistance = 2 * (fS_Utils::max3(parentSize) + fS_Utils::max3(size));
-	double currentDistance = fS_Utils::avg(maxDistance, minDistance);
-	int collisionDetected = false;
-	while (maxDistance - minDistance > PartDistanceEstimator::PRECISION)
-	{
-		Pt3D vectorBetweenParts = state->v * currentDistance;
-		collisionDetected = PartDistanceEstimator::isCollision(parentTmpPart, centers, vectorBetweenParts);
-
-		if (collisionDetected)
-		{
-			minDistance = currentDistance;
-			currentDistance = fS_Utils::avg(maxDistance, currentDistance);
-		} else
-		{
-			maxDistance = currentDistance;
-			currentDistance = fS_Utils::avg(currentDistance, minDistance);
-		}
-	}
-	delete tmpPart;
-	delete parentTmpPart;
-	return currentDistance;
-}
 
 #endif //_PART_DISTANCE_ESTIMATOR_H_
