<commit_before>#ifndef __DeformableModelBasicAdaptator__
#define __DeformableModelBasicAdaptator__

/*!
 * \file DeformableModelBasicAdaptator.h
 *
 * \brief Deformation of triangular mesh to gradient feature in an image.
 *
 * \author <NAME> - NeuroPoly (http://www.neuropoly.info)
 */

#include "Image3D.h"
#include <itkPoint.h>
#include <vector>
#include "Mesh.h"
#include "Vertex.h"
#include "util/Matrix3x3.h"
#include "util/MatrixNxM.h"
#include "SpinalCord.h"
#include <itkImageAlgorithm.h>
#include <itkImageFileReader.h>
#include <itkSingleValuedCostFunction.h>

#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataNormals.h>
#include <vtkPointData.h>
#include <vtkPoints.h>

typedef itk::CovariantVector<double,3> PixelType;
typedef itk::Image< PixelType, 3 > ImageVectorType;
typedef ImageVectorType::IndexType Index;
typedef itk::Point< double, 3 > PointType;

/*!
 * \class FoncteurDeformableBasicLocalAdaptation
 * \brief Equation class for deformable model optimization
 *
 * This class compute the value and the derivatives of the energy equation used for the optimization of the deformable model.
 */
class FoncteurDeformableBasicLocalAdaptation: public itk::SingleValuedCostFunction
{
public:
	typedef itk::SingleValuedCostFunction	Superclass;
	typedef Superclass::ParametersType		ParametersType;
	typedef Superclass::DerivativeType		DerivativeType;

	FoncteurDeformableBasicLocalAdaptation(Image3D* image, Mesh* m, ParametersType &pointsInitiaux, int nbPoints) :
		image_(image), mesh_(m), pointsInitiaux_(pointsInitiaux), nbParametres_(3*nbPoints)
	{
        verbose_ = false;
        
		listeTriangles_ = m->getListTriangles();

		type_image_factor = image->getTypeImageFactor();

		InitParameters();

		trianglesBarycentre_.resize(listeTriangles_.size()/3);

		computeOptimalPoints(pointsInitiaux);

		listeTrianglesContenantPoint = mesh_->getConnectiviteTriangles();
		pointsVoisins = mesh_->getNeighbors();
	}

	virtual void GetDerivative (const ParametersType &parameters, DerivativeType &derivative) const
	{
		//vector<CVector3> der;
		unsigned int nbTrianglesInt = listeTriangles_.size(), nbTriangles = nbTrianglesInt/3, nbPoints = nbParametres_/3;
		vector<CVector3> trianglesBarycentre(nbTriangles);
		for (unsigned int i=0; i<nbTrianglesInt; i+=3)
		{
			trianglesBarycentre[(i+1)/3] = CVector3((parameters[3*listeTriangles_[i]]+parameters[3*listeTriangles_[i+1]]+parameters[3*listeTriangles_[i+2]])/3,
													(parameters[3*listeTriangles_[i]+1]+parameters[3*listeTriangles_[i+1]+1]+parameters[3*listeTriangles_[i+2]+1])/3,
													(parameters[3*listeTriangles_[i]+2]+parameters[3*listeTriangles_[i+1]+2]+parameters[3*listeTriangles_[i+2]+2])/3);
		}

		derivative = DerivativeType(nbParametres_);

		CVector3 index, gradient, laplacien, distancePoint, expect, ci1, c1, ci2, c2, point, voisin;
		for (int i=0; i<nbPoints; i++)
		{
			point(parameters[3*i],parameters[3*i+1],parameters[3*i+2]);
			ci1 = point - transformation_*CVector3(pointsInitiaux_[3*i],pointsInitiaux_[3*i+1],pointsInitiaux_[3*i+2]);
			c1 = CVector3::ZERO; c2 = CVector3::ZERO;
			for (unsigned int j=0; j<pointsVoisins[i].size(); j++) {
				voisin(parameters[3*pointsVoisins[i][j]],parameters[3*pointsVoisins[i][j]+1],parameters[3*pointsVoisins[i][j]+2]);
				c2 += point - voisin;
				c1 += ci1 - voisin + transformation_*CVector3(pointsInitiaux_[3*pointsVoisins[i][j]],pointsInitiaux_[3*pointsVoisins[i][j]+1],pointsInitiaux_[3*pointsVoisins[i][j]+2]);
			}
			derivative[3*i] = 2*alpha*c1[0] + 2*beta*c2[0]; // Internal energy
			derivative[3*i+1] = 2*alpha*c1[1] + 2*beta*c2[1];
			derivative[3*i+2] = 2*alpha*c1[2] + 2*beta*c2[2];
			for (int k=0; k<3; k++)
			{
				for (unsigned int j=0; j<listeTrianglesContenantPoint[i].size(); j++)
				{
					expect = listeXiOpt[listeTrianglesContenantPoint[i][j]];
					if (image_->TransformPhysicalPointToContinuousIndex(expect,index))
					{
						gradient = type_image_factor*image_->GetContinuousPixelVector(index).Normalize();
						distancePoint = listeXiOpt[listeTrianglesContenantPoint[i][j]]-trianglesBarycentre[listeTrianglesContenantPoint[i][j]];
						derivative[3*i+k] += -(2.0/3.0)*listeWi[listeTrianglesContenantPoint[i][j]]*gradient[k]*(gradient*distancePoint);
					}
				}
			}
			//der.push_back(CVector3(derivative[3*indexParam[i]],derivative[3*indexParam[i]+1],derivative[3*indexParam[i]+2]));
			//cout << derivative[i] << " " << derivative[i+1] << " " << derivative[i+2] << " " << endl;
		}
	}
 
	virtual MeasureType GetValue (const ParametersType &parameters) const
	{
		double result = 0.0, interne1 = 0.0, interne2 = 0.0, externe = 0.0;
		
		unsigned int nbTrianglesInt = listeTriangles_.size(), nbTriangles = nbTrianglesInt/3, nbPoints = nbParametres_/3;
		vector<CVector3> trianglesBarycentre(nbTriangles);
		for (unsigned int i=0; i<nbTrianglesInt; i+=3)
		{
			trianglesBarycentre[(i+1)/3] = CVector3((parameters[3*listeTriangles_[i]]+parameters[3*listeTriangles_[i+1]]+parameters[3*listeTriangles_[i+2]])/3,
													(parameters[3*listeTriangles_[i]+1]+parameters[3*listeTriangles_[i+1]+1]+parameters[3*listeTriangles_[i+2]+1])/3,
													(parameters[3*listeTriangles_[i]+2]+parameters[3*listeTriangles_[i+1]+2]+parameters[3*listeTriangles_[i+2]+2])/3);
		}

		CVector3 index, gradient, expect;
		for (unsigned int i=0; i<nbTriangles; i++)
		{
			expect = listeXiOpt[i];
			if (image_->TransformPhysicalPointToContinuousIndex(expect,index))
			{
				gradient = type_image_factor*image_->GetContinuousPixelVector(index).Normalize();
				externe += listeWi[i]*pow(gradient*(listeXiOpt[i]-trianglesBarycentre[i]),2);
			}
		}
		
		CVector3 c1, c2;

		for (int i=0; i<nbPoints; i++)
		{
			for (unsigned int j=0; j<pointsVoisins[i].size(); j++)
			{
				c2 = CVector3(parameters[3*i]-parameters[3*pointsVoisins[i][j]],parameters[3*i+1]-parameters[3*pointsVoisins[i][j]+1],parameters[3*i+2]-parameters[3*pointsVoisins[i][j]+2]);//point - voisin;
				c1 = c2 - transformation_*CVector3(pointsInitiaux_[3*i]-pointsInitiaux_[3*pointsVoisins[i][j]],pointsInitiaux_[3*i+1]-pointsInitiaux_[3*pointsVoisins[i][j]+1],pointsInitiaux_[3*i+2]-pointsInitiaux_[3*pointsVoisins[i][j]+2]);
				interne1 += pow(c1[0],2)+pow(c1[1],2)+pow(c1[2],2);
				interne2 += pow(c2[0],2)+pow(c2[1],2)+pow(c2[2],2);
			}
		}

		result = externe + alpha*interne1 + beta*interne2;
		//cout << "Energie Externe = " << externe << endl << "Energie Interne = " << interne1 << endl << "Energie Totale : " << result << endl;
		return result;
	}

	virtual unsigned int GetNumberOfParameters (void) const { return nbParametres_; }

	MeasureType GetInitialValue ()
	{
		return GetValue(pointsInitiaux_);
	}

	double getInitialNormDerivatives ()
	{
		DerivativeType derivative;
		GetDerivative (pointsInitiaux_, derivative);
		double normeGradient = 0.0;
		for (int i=0; i<nbParametres_; i++)
			normeGradient += pow(derivative[i],2);
		return sqrt(normeGradient);
	}
	
	void setInitialParameters(const ParametersType &pointsInitiaux)
	{
		listeXiOpt.clear();
		listeWi.clear();

		computeOptimalPoints(pointsInitiaux);
	}

	void computeOptimalPoints(const ParametersType &pointsInitiaux)
	{
		listeXiOpt.clear();
		listeWi.clear();

		unsigned int sizeBary = trianglesBarycentre_.size(), nbTrianglesInt = listeTriangles_.size();
		CVector3 point1, normale1, point2, normale2, point3, normale3;
		for (unsigned int i=0; i<nbTrianglesInt; i+=3) {
			point1(pointsInitiaux[3*listeTriangles_[i]],pointsInitiaux[3*listeTriangles_[i]+1],pointsInitiaux[3*listeTriangles_[i]+2]);
			point2(pointsInitiaux[3*listeTriangles_[i+1]],pointsInitiaux[3*listeTriangles_[i+1]+1],pointsInitiaux[3*listeTriangles_[i+1]+2]);
			point3(pointsInitiaux[3*listeTriangles_[i+2]],pointsInitiaux[3*listeTriangles_[i+2]+1],pointsInitiaux[3*listeTriangles_[i+2]+2]);
			trianglesBarycentre_[(i+1)/3] = Vertex((point1+point2+point3)/3,-((point1-point2)^(point1-point3)).Normalize());
		}
		
        // define startPos as meanRadius of the mesh
        int startPos = 0;
        
        
		double resultCkMax = 0.0, resultCk;
		CVector3 xi, ni, ck, ci, index, gradient;
		int k;
        
        Matrice imageDistance = Matrice(sizeBary,2*line_search+1);

		listeXiOpt.resize(sizeBary);
		listeWi.resize(sizeBary);
		vector<double> listeDistancePointsOpt(sizeBary);
        
		for (unsigned int i=0; i<sizeBary; i++)
		{
			xi = trianglesBarycentre_[i].getPosition();
			ni = trianglesBarycentre_[i].getNormal();
			k = 0;
			resultCkMax = 0.0;
			for (int j=-startPos; j<=line_search; j++) {
				if (image_->TransformPhysicalPointToContinuousIndex(xi + j*deltaNormale*ni,index)) {
					resultCk = type_image_factor*ni*image_->GetContinuousPixelVector(index) - tradeOff*deltaNormale*deltaNormale*j*j;
                    //imageDistance(i,j+startPos) = resultCk;<commit_msg>BUG: fixed bug on sct_propseg<commit_after>#ifndef __DeformableModelBasicAdaptator__
#define __DeformableModelBasicAdaptator__

/*!
 * \file DeformableModelBasicAdaptator.h
 *
 * \brief Deformation of triangular mesh to gradient feature in an image.
 *
 * \author <NAME> - NeuroPoly (http://www.neuropoly.info)
 */

#include "Image3D.h"
#include <itkPoint.h>
#include <vector>
#include "Mesh.h"
#include "Vertex.h"
#include "util/Matrix3x3.h"
#include "util/MatrixNxM.h"
#include "SpinalCord.h"
#include <itkImageAlgorithm.h>
#include <itkImageFileReader.h>
#include <itkSingleValuedCostFunction.h>

#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataNormals.h>
#include <vtkPointData.h>
#include <vtkPoints.h>

typedef itk::CovariantVector<double,3> PixelType;
typedef itk::Image< PixelType, 3 > ImageVectorType;
typedef ImageVectorType::IndexType Index;
typedef itk::Point< double, 3 > PointType;

/*!
 * \class FoncteurDeformableBasicLocalAdaptation
 * \brief Equation class for deformable model optimization
 *
 * This class compute the value and the derivatives of the energy equation used for the optimization of the deformable model.
 */
class FoncteurDeformableBasicLocalAdaptation: public itk::SingleValuedCostFunction
{
public:
	typedef itk::SingleValuedCostFunction	Superclass;
	typedef Superclass::ParametersType		ParametersType;
	typedef Superclass::DerivativeType		DerivativeType;

	FoncteurDeformableBasicLocalAdaptation(Image3D* image, Mesh* m, ParametersType &pointsInitiaux, int nbPoints) :
		image_(image), mesh_(m), pointsInitiaux_(pointsInitiaux), nbParametres_(3*nbPoints)
	{
        verbose_ = false;
        
		listeTriangles_ = m->getListTriangles();

		type_image_factor = image->getTypeImageFactor();

		InitParameters();

		trianglesBarycentre_.resize(listeTriangles_.size()/3);

		computeOptimalPoints(pointsInitiaux);

		listeTrianglesContenantPoint = mesh_->getConnectiviteTriangles();
		pointsVoisins = mesh_->getNeighbors();
	}

	virtual void GetDerivative (const ParametersType &parameters, DerivativeType &derivative) const
	{
		//vector<CVector3> der;
		unsigned int nbTrianglesInt = listeTriangles_.size(), nbTriangles = nbTrianglesInt/3, nbPoints = nbParametres_/3;
		vector<CVector3> trianglesBarycentre(nbTriangles);
		for (unsigned int i=0; i<nbTrianglesInt; i+=3)
		{
			trianglesBarycentre[(i+1)/3] = CVector3((parameters[3*listeTriangles_[i]]+parameters[3*listeTriangles_[i+1]]+parameters[3*listeTriangles_[i+2]])/3,
													(parameters[3*listeTriangles_[i]+1]+parameters[3*listeTriangles_[i+1]+1]+parameters[3*listeTriangles_[i+2]+1])/3,
													(parameters[3*listeTriangles_[i]+2]+parameters[3*listeTriangles_[i+1]+2]+parameters[3*listeTriangles_[i+2]+2])/3);
		}

		derivative = DerivativeType(nbParametres_);

		CVector3 index, gradient, laplacien, distancePoint, expect, ci1, c1, ci2, c2, point, voisin;
		for (int i=0; i<nbPoints; i++)
		{
			point(parameters[3*i],parameters[3*i+1],parameters[3*i+2]);
			ci1 = point - transformation_*CVector3(pointsInitiaux_[3*i],pointsInitiaux_[3*i+1],pointsInitiaux_[3*i+2]);
			c1 = CVector3::ZERO; c2 = CVector3::ZERO;
			for (unsigned int j=0; j<pointsVoisins[i].size(); j++) {
				voisin(parameters[3*pointsVoisins[i][j]],parameters[3*pointsVoisins[i][j]+1],parameters[3*pointsVoisins[i][j]+2]);
				c2 += point - voisin;
				c1 += ci1 - voisin + transformation_*CVector3(pointsInitiaux_[3*pointsVoisins[i][j]],pointsInitiaux_[3*pointsVoisins[i][j]+1],pointsInitiaux_[3*pointsVoisins[i][j]+2]);
			}
			derivative[3*i] = 2*alpha*c1[0] + 2*beta*c2[0]; // Internal energy
			derivative[3*i+1] = 2*alpha*c1[1] + 2*beta*c2[1];
			derivative[3*i+2] = 2*alpha*c1[2] + 2*beta*c2[2];
			for (int k=0; k<3; k++)
			{
				for (unsigned int j=0; j<listeTrianglesContenantPoint[i].size(); j++)
				{
					expect = listeXiOpt[listeTrianglesContenantPoint[i][j]];
					if (image_->TransformPhysicalPointToContinuousIndex(expect,index))
					{
						gradient = type_image_factor*image_->GetContinuousPixelVector(index).Normalize();
						distancePoint = listeXiOpt[listeTrianglesContenantPoint[i][j]]-trianglesBarycentre[listeTrianglesContenantPoint[i][j]];
						derivative[3*i+k] += -(2.0/3.0)*listeWi[listeTrianglesContenantPoint[i][j]]*gradient[k]*(gradient*distancePoint);
					}
				}
			}
			//der.push_back(CVector3(derivative[3*indexParam[i]],derivative[3*indexParam[i]+1],derivative[3*indexParam[i]+2]));
			//cout << derivative[i] << " " << derivative[i+1] << " " << derivative[i+2] << " " << endl;
		}
	}
 
	virtual MeasureType GetValue (const ParametersType &parameters) const
	{
		double result = 0.0, interne1 = 0.0, interne2 = 0.0, externe = 0.0;
		
		unsigned int nbTrianglesInt = listeTriangles_.size(), nbTriangles = nbTrianglesInt/3, nbPoints = nbParametres_/3;
		vector<CVector3> trianglesBarycentre(nbTriangles);
		for (unsigned int i=0; i<nbTrianglesInt; i+=3)
		{
			trianglesBarycentre[(i+1)/3] = CVector3((parameters[3*listeTriangles_[i]]+parameters[3*listeTriangles_[i+1]]+parameters[3*listeTriangles_[i+2]])/3,
													(parameters[3*listeTriangles_[i]+1]+parameters[3*listeTriangles_[i+1]+1]+parameters[3*listeTriangles_[i+2]+1])/3,
													(parameters[3*listeTriangles_[i]+2]+parameters[3*listeTriangles_[i+1]+2]+parameters[3*listeTriangles_[i+2]+2])/3);
		}

		CVector3 index, gradient, expect;
		for (unsigned int i=0; i<nbTriangles; i++)
		{
			expect = listeXiOpt[i];
			if (image_->TransformPhysicalPointToContinuousIndex(expect,index))
			{
				gradient = type_image_factor*image_->GetContinuousPixelVector(index).Normalize();
				externe += listeWi[i]*pow(gradient*(listeXiOpt[i]-trianglesBarycentre[i]),2);
			}
		}
		
		CVector3 c1, c2;

		for (int i=0; i<nbPoints; i++)
		{
			for (unsigned int j=0; j<pointsVoisins[i].size(); j++)
			{
				c2 = CVector3(parameters[3*i]-parameters[3*pointsVoisins[i][j]],parameters[3*i+1]-parameters[3*pointsVoisins[i][j]+1],parameters[3*i+2]-parameters[3*pointsVoisins[i][j]+2]);//point - voisin;
				c1 = c2 - transformation_*CVector3(pointsInitiaux_[3*i]-pointsInitiaux_[3*pointsVoisins[i][j]],pointsInitiaux_[3*i+1]-pointsInitiaux_[3*pointsVoisins[i][j]+1],pointsInitiaux_[3*i+2]-pointsInitiaux_[3*pointsVoisins[i][j]+2]);
				interne1 += pow(c1[0],2)+pow(c1[1],2)+pow(c1[2],2);
				interne2 += pow(c2[0],2)+pow(c2[1],2)+pow(c2[2],2);
			}
		}

		result = externe + alpha*interne1 + beta*interne2;
		//cout << "Energie Externe = " << externe << endl << "Energie Interne = " << interne1 << endl << "Energie Totale : " << result << endl;
		return result;
	}

	virtual unsigned int GetNumberOfParameters (void) const { return nbParametres_; }

	MeasureType GetInitialValue ()
	{
		return GetValue(pointsInitiaux_);
	}

	double getInitialNormDerivatives ()
	{
		DerivativeType derivative;
		GetDerivative (pointsInitiaux_, derivative);
		double normeGradient = 0.0;
		for (int i=0; i<nbParametres_; i++)
			normeGradient += pow(derivative[i],2);
		return sqrt(normeGradient);
	}
	
	void setInitialParameters(const ParametersType &pointsInitiaux)
	{
		listeXiOpt.clear();
		listeWi.clear();

		computeOptimalPoints(pointsInitiaux);
	}

	void computeOptimalPoints(const ParametersType &pointsInitiaux)
	{
		listeXiOpt.clear();
		listeWi.clear();

		unsigned int sizeBary = trianglesBarycentre_.size(), nbTrianglesInt = listeTriangles_.size();
		CVector3 point1, normale1, point2, normale2, point3, normale3;
		for (unsigned int i=0; i<nbTrianglesInt; i+=3) {
			point1(pointsInitiaux[3*listeTriangles_[i]],pointsInitiaux[3*listeTriangles_[i]+1],pointsInitiaux[3*listeTriangles_[i]+2]);
			point2(pointsInitiaux[3*listeTriangles_[i+1]],pointsInitiaux[3*listeTriangles_[i+1]+1],pointsInitiaux[3*listeTriangles_[i+1]+2]);
			point3(pointsInitiaux[3*listeTriangles_[i+2]],pointsInitiaux[3*listeTriangles_[i+2]+1],pointsInitiaux[3*listeTriangles_[i+2]+2]);
			trianglesBarycentre_[(i+1)/3] = Vertex((point1+point2+point3)/3,-((point1-point2)^(point1-point3)).Normalize());
		}
		
        // define startPos as meanRadius of the mesh
        int startPos = line_search;
        
		double resultCkMax = 0.0, resultCk;
		CVector3 xi, ni, ck, ci, index, gradient;
		int k;
        
        Matrice imageDistance = Matrice(sizeBary,2*line_search+1);

		listeXiOpt.resize(sizeBary);
		listeWi.resize(sizeBary);
		vector<double> listeDistancePointsOpt(sizeBary);
        
		for (unsigned int i=0; i<sizeBary; i++)
		{
			xi = trianglesBarycentre_[i].getPosition();
			ni = trianglesBarycentre_[i].getNormal();
			k = 0;
			resultCkMax = 0.0;
			for (int j=-startPos; j<=line_search; j++) {
				if (image_->TransformPhysicalPointToContinuousIndex(xi + j*deltaNormale*ni,index)) {
					resultCk = type_image_factor*ni*image_->GetContinuousPixelVector(index) - tradeOff*deltaNormale*deltaNormale*j*j;
                    //imageDistance(i,j+startPos) = resultCk;
{'fraction_non_alphanumeric': 0.15327725399424152, 'fraction_numerical': 0.030768362219838537, 'mean_word_length': 5.903351519875292, 'pattern_counts': {'":': 0, '<': 131, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3513797', 'n_tokens_mistral': 7327, 'n_tokens_neox': 6586, 'n_words': 945}