wxSphereView.cxx

Go to the documentation of this file.

#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;
}