/*# ---------------------------------------------------------------------
#
# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
#                        pour la Sant�)
# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
# Previous Authors : Laurent Guigues, Jean-Pierre Roux
# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
#
#  This software is governed by the CeCILL-B license under French law and
#  abiding by the rules of distribution of free software. You can  use,
#  modify and/ or redistribute the software under the terms of the CeCILL-B
#  license as circulated by CEA, CNRS and INRIA at the following URL
#  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
#  or in the file LICENSE.txt.
#
#  As a counterpart to the access to the source code and  rights to copy,
#  modify and redistribute granted by the license, users are provided only
#  with a limited warranty  and the software's author,  the holder of the
#  economic rights,  and the successive licensors  have only  limited
#  liability.
#
#  The fact that you are presently reading this means that you have had
#  knowledge of the CeCILL-B license and that you accept its terms.
# ------------------------------------------------------------------------ */

#include "wxSphereView.h"


wxSphereView::wxSphereView( wxWindow *parent, vtkMPRBaseData *vtkmprbasedata, vtkImageData *imageData )
: wxVtk2DBaseView(parent)
{
	_delta				=	1;
	_vtkmprbasedata		=	vtkmprbasedata;
	_imageDataOriginal	=	imageData;

	_imageSphere		=	vtkImageData::New();
	_imageSphere->SetDimensions (150,150,500);
	_imageSphere->SetScalarTypeToUnsignedShort();
	_imageSphere->AllocateScalars();   
	_imageSphere->Update();   


	vtkBaseData *vtkbasedata = new vtkBaseData();
	vtkbasedata->SetMarImageData( new marImageData(_imageSphere) );
	this->SetVtkBaseData(vtkbasedata);

    _transform			=	vtkTransform::New();
    _transform1			=	vtkTransform::New();
    _transform2			=	vtkTransform::New();
	_transform ->Identity();
	_transform1->Identity();
	_transform2->Identity();

	_radio=25;
}

//-------------------------------------------------------------------

wxSphereView::~wxSphereView()
{
	_transform  -> Delete();
	_transform1 -> Delete();
	_transform2 -> Delete();
	ResetlstId();
}

//----------------------------------------------------------------------------

double wxSphereView::GetRadio()
{
	return _radio;
}

//----------------------------------------------------------------------------

void wxSphereView::SetRadio(double radio)
{
	if (radio<0)
	{
		radio=0;
	}
	_radio=radio;
}

//----------------------------------------------------------------------------

void wxSphereView::Configure()
{
	wxVtk2DBaseView::Configure();

	_vtkinteractorstylesphere = new vtkInteractorStyleSphere();
	((vtkInteractorStyleBaseView*)GetInteractorStyleBaseView())->AddInteractorStyleMaracas( _vtkinteractorstylesphere );
	double points[4][3];

// EED purify 12/sep/2006
	int i,j;
	for (i=0;i<4;i++)
	{
		for (j=0;j<3;j++)
		{
			points[i][j]=0;
		}
	}

	InitSphere(points);
	DefineImageSphere();
}

//----------------------------------------------------------------------------

void wxSphereView::RefreshPoint()
{
	double x	= _vtkmprbasedata->GetX() - _centerX;
	double y	= _vtkmprbasedata->GetY() - _centerY;
	double z	= _vtkmprbasedata->GetZ() - _centerZ;
	double alpha= atan2(x,z);
	double beta = atan2( y , sqrt(z*z+x*x) );

	alpha		= alpha*180/3.1416;
	beta		= beta*180/3.1416;

	_transform1->Identity();
	_transform1->RotateY(alpha);
	_transform1->RotateX(-beta);

	_radio= sqrt(x*x + y*y +z*z);

	RefreshView();
}

//----------------------------------------------------------------------------

void wxSphereView::RefreshView()
{
	DefineImageSphere();
	wxVtk2DBaseView::Refresh();
}

//----------------------------------------------------------------------------

void wxSphereView::RotationEnd()
{
	_transform1->RotateWXYZ(_ang,_vxb,_vyb,0);
	_transform2->Identity();
	SetDeltaVoxel(1);
	ResetlstId();
}

//----------------------------------------------------------------------------

void wxSphereView::RotationStart(double vx, double vy, bool ok_v, bool ok_ang)
{
	if (ok_ang==false)
	{
		_ang = -sqrt( vx*vx + vy*vy ) / 1.0;
	}

	if (ok_v==false){
		_vxb=-vy;
		_vyb=vx;
	}

	_transform2->Identity();
	_transform2->RotateWXYZ(_ang,_vxb,_vyb,0);
	SetDeltaVoxel(3);
	ResetlstId();
}

//----------------------------------------------------------------------------

void wxSphereView::GetPointSphere(double p[3],double r1,double angA,double angB)
{
	double in[3],out[3];
	in[0]=0;   
	in[1]=r1;   
	in[2]=0;
	vtkTransform *transform = vtkTransform::New();
	transform->Identity();
	transform->RotateX(angB);
	transform->RotateZ(angA);
	transform->TransformPoint(in,out);
	p[0]=out[0];
	p[1]=out[1];
	p[2]=out[2];
	transform->Delete(); 
}

//----------------------------------------------------------------------------

void wxSphereView::RotatePointOverTheSphere( double pp[3], double p[3],double cc[3])
{

	double out[3];
	_transform->TransformPoint(p,out);
	pp[0] = out[0] + cc[0];
	pp[1] = out[1] + cc[1];
	pp[2] = out[2] + cc[2];

}

//----------------------------------------------------------------------------

void wxSphereView::TransferePoints(double pp1[3],double pp2[3],double AngX,double AngY,vtkImageData *image)
{
	double t;
	double difX = pp2[0]-pp1[0];
	double difY = pp2[1]-pp1[1];
	double difZ = pp2[2]-pp1[2];

	double	max		= 200;

	int dimOrg[3];
	int dimRes[3];
	int z;
	_imageDataOriginal->GetDimensions(dimOrg);		
	image->GetDimensions(dimRes);		

	int i;
	double x1=pp1[0];
	double y1=pp1[1];
	double z1=pp1[2];
	int xx=-1,yy=-1,zz=-1;

	for (i=0;i<max;i++)
	{
		t  = i/max;
		xx = (int) (x1+t*difX);
		yy = (int) (y1+t*difY);
		zz = (int) (z1+t*difZ);

		z=i;
		if ((xx>=0) && (xx<dimOrg[0]) && (yy>=0) && (yy<dimOrg[1]) && (zz>=0) && (zz<dimOrg[2]) &&
			(AngX>=0) && (AngX<dimRes[0]) && (AngY>=0) && (AngY<dimRes[1]) && (z>=0) && (z<dimRes[2]) )
		{
			unsigned short *pOrg=(unsigned short*)_imageDataOriginal->GetScalarPointer (xx,yy,zz); 
			unsigned short *pRes=(unsigned short*)image->GetScalarPointer( (int)AngX , (int)AngY , z ); 
			*pRes=*pOrg;
		}
	}
}

//----------------------------------------------------------------------------

void wxSphereView::ResetlstId()
{
	int i,size=_lstId.size();
	for (i=size-1;i>=0;i--)
	{
		delete _lstId[i];
	}
	_lstId.clear();
}

//----------------------------------------------------------------------------

int wxSphereView::GetIdOfImage(double radio)
{
	int id=0;
	int dim[3];	
	_imageSphere->GetDimensions(dim);	
	int sizeMaxList = dim[2];
	// Search in list >> alpha beta radio
	int i,size=_lstId.size();
	for (i=0; i<size;i++)
	{
		//idAlBeRa *tmp=_lstId[i]; // JPRx
		if ((_lstId[i]->_radio==radio) && (_lstId[i]->_deltavoxel==_delta)) 
		{
			return _lstId[i]->_id;
		}
	}
	if (size>sizeMaxList)
	{
		delete _lstId[size-1];
		_lstId.pop_back(); 
	}
	if (size!=0){
		id=_lstId[0]->_id+1;
		id = id % sizeMaxList;
	} else {
		id = 0;
	}

	FiltreImage(id,radio);
	_lstId.insert(_lstId.begin(),1,new idAlBeRa(id,radio,_delta) ); 

	return id;
}

//----------------------------------------------------------------------------

void wxSphereView::DefineImageSphere()
{
	int id;
	id=GetIdOfImage( _radio );
	GetVtkBaseData()->SetZ( id );
}


//----------------------------------------------------------------------------
void wxSphereView::SetDeltaVoxel(int delta)
{
	_delta=delta;
}

//----------------------------------------------------------------------------
void wxSphereView::SetVoxel(double i, double j, int delta,double id,  unsigned short gris)
{
	int ii,jj,delta2;
	unsigned short *pRes;
	int dimRes[3];
	_imageSphere->GetDimensions(dimRes);

	delta2=delta-1;
	for ( ii=(int)(i-delta2) ; ii<=(int)(i+delta2) ; ii++ )
	{
		for ( jj=(int)(j-delta2) ; jj<=(int)(j+delta2) ; jj++ )
		{
			if ( (ii>=0)&&(ii<dimRes[0]) &&  
				 (jj>=0)&&(jj<dimRes[1]) )
			{
				pRes = (unsigned short*)_imageSphere->GetScalarPointer( ii , jj , (int)id );
				*pRes=gris;
			}
		}
	}

}

//----------------------------------------------------------------------------

void wxSphereView::SetXYZtoParent(double i, double j)
{

	double factor = 0.75;
	double radio2	= _radio*_radio;
	double pxx,pyy,d2x,d2y;
	double cc[3],p[3],pp[3];
	cc[0]=_centerX;
	cc[1]=_centerY;
	cc[2]=_centerZ;
	double aa;
	int dimRes[3],dimOrig[3];
	_imageSphere->GetDimensions(dimRes);
	d2x=dimRes[0]/2;
	d2y=dimRes[1]/2;
	_imageDataOriginal->GetDimensions(dimOrig);

	p[0]  = (i - d2x)*factor;
	pxx=p[0]*p[0];
	p[1]  = (j - d2y)*factor;
	pyy=p[1]*p[1];
	aa = pxx + pyy;
	if (radio2>aa){
		aa=radio2-aa;
		p[2]  = sqrt(aa);
		RotatePointOverTheSphere( pp, p,cc);
		if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) && 
			 (pp[1]>=0) && (pp[1]<dimOrig[1]) && 
			 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
		{
			if (_vtkmprbasedata){
				_vtkmprbasedata->SetX(pp[0]);
				_vtkmprbasedata->SetY(pp[1]);
				_vtkmprbasedata->SetZ(pp[2]);
			}
		}
	}
}


//----------------------------------------------------------------------------

void wxSphereView::FiltreImageB(int id, double radio, bool ok,int deltaTMP)
{	
	double factor = 0.75;
	double radioB	= radio/3;
	double radio2	= radio*radio;
	double pxx,pyy,d2x,d2y;
	double cc[3],p[3],pp[3];
	cc[0]=_centerX;
	cc[1]=_centerY;
	cc[2]=_centerZ;
	double aa;
	unsigned short *pOrig;
	int dimRes[3],dimOrig[3];
	double i,j;
	int ext[6];
	_imageSphere->GetExtent(ext);
	_imageSphere->GetDimensions(dimRes);
	//JCP 24 - 04 -09
	//_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
	_imageSphere->SetExtent(ext);
	//JCP 24 - 04 -09
	d2x=dimRes[0]/2;
	d2y=dimRes[1]/2;
//	double deltaTMP=_delta;
	_imageDataOriginal->GetDimensions(dimOrig);

	int start,end;
	int limitA,limitB;
	limitA	= (int) ( (-radioB/factor)+d2x );
	limitB	= (int) ( (radioB/factor)+d2x );
	if (ok==true){
		start	= limitA;
		end		= limitB;
	} else {
		start=0;
		end=dimRes[0];
	}

	for ( i=start ; i<end ; i=i+deltaTMP )
	{
		p[0]  = (i - d2x)*factor;
		pxx=p[0]*p[0];
		for (j=start;j<end;j=j+deltaTMP)
		{
			p[1]  = (j - d2y)*factor;
			pyy=p[1]*p[1];
			aa = pxx + pyy;

			if  (( ((ok==false) && (!((i>limitA) && (i<limitB) && (j>limitA) && (j<limitB)))) )
				    ||
					(ok==true))
			{
				if (radio2>aa){
					aa=radio2-aa;
					p[2]  = sqrt(aa);
					RotatePointOverTheSphere( pp, p,cc);
					if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) && 
						 (pp[1]>=0) && (pp[1]<dimOrig[1]) && 
						 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
					{
						pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( (int)(pp[0]) , (int)(pp[1]) , (int)(pp[2]) ); 
						SetVoxel(i,j,deltaTMP,id,*pOrig);
					} else {
						SetVoxel(i,j,deltaTMP,id,2000);
					}
				} else {
					SetVoxel(i,j,deltaTMP,id,0);
				}
			}
		}
	}

	_imageSphere->Modified();  
	_imageSphere->Update();
}




//----------------------------------------------------------------------------

void wxSphereView::FiltreImage(int id, double radio)
{

	_transform -> Identity();
	_transform -> Concatenate(_transform1);
	_transform -> Concatenate(_transform2);

	FiltreImageB(id,radio,false, _delta);
	FiltreImageB(id,radio,true, 1);
}


//----------------------------------------------------------------------------

/*
void wxSphereView::FiltreImage(int id, double radio)
{
	double radio2	= radio*radio;
	double radio2TMP= (radio/2)*(radio/2);
	double pxx,pyy,d2x,d2y;
	double cc[3],p[3],pp[3];
	cc[0]=_centerX;
	cc[1]=_centerY;
	cc[2]=_centerZ;
	double aa;
	unsigned short *pOrig;
	int dimRes[3],dimOrig[3];
	double i,j;
	_imageSphere->GetDimensions(dimRes);
	_imageSphere->SetExtent(0,dimRes[0]-1,0,dimRes[1]-1,0,dimRes[2]-1);
	d2x=dimRes[0]/2;
	d2y=dimRes[1]/2;
	double deltaTMP=_delta;
	_imageDataOriginal->GetDimensions(dimOrig);

	for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
	{
		p[0]  = (i - d2x)*0.75;
		pxx=p[0]*p[0];
		for (j=0;j<dimRes[1];j=j+deltaTMP)
		{
			p[1]  = (j - d2y)*0.75;
			pyy=p[1]*p[1];
			aa = pxx + pyy;

			if (aa>radio2TMP)
			{
				if (radio2>aa){
					aa=radio2-aa;
					p[2]  = sqrt(aa);
					RotatePointOverTheSphere( pp, p,cc);
					if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) && 
						 (pp[1]>=0) && (pp[1]<dimOrig[1]) && 
						 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
					{
						pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] ); 
						SetVoxel(i,j,deltaTMP,id,*pOrig);
					} else {
						SetVoxel(i,j,deltaTMP,id,2000);
					}
				} else {
					SetVoxel(i,j,deltaTMP,id,0);
				}
			}
		}
	}


	deltaTMP=1;
	for ( i=0 ; i<dimRes[0] ; i=i+deltaTMP )
	{
		p[0]  = (i - d2x)*0.75;
		pxx=p[0]*p[0];
		for (j=0;j<dimRes[1];j=j+deltaTMP)
		{
			p[1]  = (j - d2y)*0.75;
			pyy=p[1]*p[1];
			aa = pxx + pyy;
			if (aa<=radio2TMP)
			{
				if (radio2>aa){
					aa=radio2-aa;
					p[2]  = sqrt(aa);
					RotatePointOverTheSphere( pp, p,cc);
					if ( (pp[0]>=0) && (pp[0]<dimOrig[0]) && 
						 (pp[1]>=0) && (pp[1]<dimOrig[1]) && 
						 (pp[2]>=0) && (pp[2]<dimOrig[2]) )
					{
						pOrig=(unsigned short*)_imageDataOriginal->GetScalarPointer( pp[0] , pp[1] , pp[2] ); 
						SetVoxel(i,j,deltaTMP,id,*pOrig);
					} else {
						SetVoxel(i,j,deltaTMP,id,2000);
					}
				} else {
					SetVoxel(i,j,deltaTMP,id,0);
				}
			}
		}
	}

	_imageSphere->Modified();  
	_imageSphere->Update();
}
*/
/*
void wxSphereView::FiltreImage(vtkImageData *imageSphere)
{
	int dim[3],i,j,k;
	imageSphere->GetDimensions(dim);
	for (i=0;i<dim[0];i++)
	{
		for (j=0;j<dim[1];j++)
		{
			for (k=0;k<dim[2];k++)
			{
				unsigned short *pRes=(unsigned short*)imageSphere->GetScalarPointer (i,j,k); 
				*pRes=0;
			}
		}
	}

	double deltaA=90;
	double cc[3],p1[3],p2[3],pp1[3],pp2[3];
	cc[0]=_centerX;
	cc[1]=_centerY;
	cc[2]=_centerZ;
	double r1	= _sl_radio->GetValue() - _sl_thickness->GetValue()/2;
	double r2	= _sl_radio->GetValue() + _sl_thickness->GetValue()/2;
	if (r1<10)
	{
		r1=10;
	}
	double alpha= _sl_alpha->GetValue();
	double beta	= _sl_beta->GetValue();

	double angA,angB;
	for (angA=-deltaA;angA<deltaA;angA++)
	{
		for (angB=-deltaA;angB<deltaA;angB++)
		{
			GetPointSphere(p1,r1,angA,angB);
			GetPointSphere(p2,r2,angA,angB);
			RotatePointOverTheSphere( pp1, alpha, beta, p1 ,cc );
			RotatePointOverTheSphere( pp2, alpha, beta, p2 ,cc );
			TransferePoints(pp1,pp2,angA+alpha+180,angB+beta+90,imageSphere);
		}
	}
}
*/


//----------------------------------------------------------------------------

void wxSphereView::InitSphere(double points[4][3])
{
	double cc[3];
    double r = SphereFindCenter(points,cc); // 4-points , center
    if (r > 0)
    {
        _centerX	= (int)(cc[0]);
        _centerY	= (int)(cc[1]);
        _centerZ	= (int)(cc[2]);
    } else {
		int dim[3];
		_imageDataOriginal->GetDimensions(dim);
        _centerX	= (int)(dim[0]/2);
        _centerY	= (int)(dim[1]/2);
        _centerZ	= (int)(dim[2]/2);
	}
}

//----------------------------------------------------------------------------

// Calculate center and radius of sphere given four points
// http://home.att.net/~srschmitt/script_sphere_solver.html
// source code HTML <script language=JavaScript>
double wxSphereView::SphereFindCenter(double P[4][3], double cc[3])
{
    int i;
    double r, m11, m12, m13, m14, m15;
	double a[4][4];

    for (i = 0; i < 4; i++)                    // find minor 11
    {
        a[i][0] = P[i][0];
        a[i][1] = P[i][1];
        a[i][2] = P[i][2];
        a[i][3] = 1;
    }
    m11 = determinant( a, 4 );

    for (i = 0; i < 4; i++)                    // find minor 12 
    {
        a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
        a[i][1] = P[i][1];
        a[i][2] = P[i][2];
        a[i][3] = 1;
    }
    m12 = determinant( a, 4 );

    for (i = 0; i < 4; i++)                    // find minor 13
    {
        a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
        a[i][1] = P[i][0];
        a[i][2] = P[i][2];
        a[i][3] = 1;
    }
    m13 = determinant( a, 4 );

    for (i = 0; i < 4; i++)                    // find minor 14
    {
        a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
        a[i][1] = P[i][0];
        a[i][2] = P[i][1];
        a[i][3] = 1;
    }
    m14 = determinant( a, 4 );


    for (i = 0; i < 4; i++)                    // find minor 15
    {
        a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1] + P[i][2]*P[i][2];
        a[i][1] = P[i][0];
        a[i][2] = P[i][1];
        a[i][3] = P[i][2];
    }
    m15 = determinant( a, 4 );

    if (m11 == 0)
    {
        r = 0;
    }
    else
    {
		// center of sphere
        cc[0] =  0.5*m12/m11;	//cx                  
        cc[1] = -0.5*m13/m11;	//cy
        cc[2] =  0.5*m14/m11;	//cz
		// Sphere radio 
        r  = sqrt( cc[0]*cc[0] + cc[1]*cc[1] + cc[2]*cc[2] - m15/m11 );
    }

    return r;                                  // the radius
}
//----------------------------------------------------------------------------

//  Recursive definition of determinate using expansion by minors.
double wxSphereView::determinant(double a[4][4], int n)
{
    int i, j, j1, j2;
    double d;
	double m[4][4];

	for (i=0;i<4;i++)
	{
		for (j=0;j<4;j++)
		{
			m[i][j]=0;
		}
	}

    if (n == 2)                                // terminate recursion
    {
        d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
    }
    else 
    {
        d = 0;
        for (j1 = 0; j1 < n; j1++ )            // do each column
        {
            for (i = 1; i < n; i++)            // create minor
            {
                j2 = 0;
                for (j = 0; j < n; j++)
                {
                    if (j == j1) continue;
                    m[i-1][j2] = a[i][j];
                    j2++;
                }
            }
            
            // sum (+/-)cofactor * minor  
            d = d + pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );
        }
    }

    return d;
}