wxPathologyWidget_01.cxx

Go to the documentation of this file.

#include "vtkActor.h"
#include "vtkAppendFilter.h"
#include "vtkCamera.h"
#include "vtkCallbackCommand.h"
#include "vtkCellArray.h"
#include "vtkCylinderSource.h"
#include "vtkClosePolyData.h"
#include "vtkCommand.h"
#include "vtkDoubleArray.h"
#include "vtkExtractVOI.h"
#include "vtkImageReader.h"
#include "vtkImageViewer.h"
#include "vtkImageViewer2.h"
#include "vtkImageToStructuredPoints.h"
#include "vtkImageThreshold.h" 
#include "vtkImageWriter.h"
#include "vtkImageClip.h"
#include "vtkImageResample.h"
#include "vtkImageReslice.h"
#include "vtkImageResample.h"
#include "vtkImageThreshold.h"
#include "vtkImageCast.h"
#include "vtkImageData.h"
#include "vtkMath.h"
#include "vtkMarchingCubes.h"
#include "vtkObjectFactory.h"
#include "vtkPolyDataMapper.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkPolyDataConnectivityFilter.h"
#include "vtkProperty.h"
#include "vtkPriorityQueue.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkPolyDataWriter.h"
#include "vtkPolyDataReader.h"
#include "vtkPointData.h"
#include "vtkSphereSource.h"
#include "vtkSTLWriter.h"
#include "vtkStripper.h"
#include "vtkTriangleFilter.h"
#include "vtkTransform.h"

#include "wxPathologyWidget_01.h"
#include "kernel/vtkOtsuSphereSource.h"
#include <wx/splitter.h>




//-------------------------------------------------------------------
//-------------------------------------------------------------------
//-------------------------------------------------------------------
wxPathologyWidget_01::wxPathologyWidget_01(wxWindow *parent, marInterface* mar)
: wxPanel( parent, -1) 
{
        

        _mar                                                            = mar;
        wxBoxSizer                      *sizer                  = new wxBoxSizer(wxVERTICAL  );
    wxSplitterWindow    *pnlSplitter    = new wxSplitterWindow( this , -1);
        wxPanel                         *viewPanel              = CreateViewPanel(pnlSplitter);
        wxPanel                         *controlPanel   = CreateControlPanel(pnlSplitter);

        sizer           -> Add( pnlSplitter ,1,wxGROW  ,0);
        pnlSplitter     -> SetMinimumPaneSize( 150 );
    pnlSplitter -> SplitVertically( viewPanel, controlPanel, 600 );
        this            -> SetSizer(sizer);



        //DHC STL SURFACES
        stlInterna = NULL;
        stlExterna = NULL;
        stlImageData = NULL;
        workingDir = NULL;
        
        //PATHOLOGY
        p1SphereSource = NULL;
        p1Mapper = NULL;
        p1Actor = NULL;

        p2SphereSource = NULL;
        p2Mapper = NULL;
        p2Actor = NULL;

        patSphereSource = NULL;
        patMapper = NULL;
        patActor = NULL;


        pathologyFrame  = NULL;

        patImageData    = NULL;
        voiFilter               = NULL;
        cubesFilter             = NULL;
        dataMapper              = NULL;
        dataActor               = NULL;

        outlineFilter   = NULL;
        outlineMapper   = NULL;
        outlineActor    = NULL;


        p3SphereSource  = NULL;
        p3Mapper                = NULL;
        p3Actor                 = NULL;
        
        
}
//-------------------------------------------------------------------
wxPathologyWidget_01::~wxPathologyWidget_01(){
        
}
//-------------------------------------------------------------------
wxPanel* wxPathologyWidget_01::CreateViewPanel(wxWindow *parent)
{
        wxPanel *panel          =       new wxPanel(parent,-1);
    wxBoxSizer *sizer   =       new wxBoxSizer(wxVERTICAL);
        _maracasSurfaceWidget = new wxSurfaceWidget(panel);
        sizer->Add(_maracasSurfaceWidget        , 1, wxEXPAND, 0);
        panel->SetSizer(sizer);
        panel->SetAutoLayout(true);
        panel->SetSize(400,400);
        panel->Layout();
        return panel;
}
//-------------------------------------------------------------------
wxPanel* wxPathologyWidget_01::CreateControlPanel(wxWindow *parent)
{
        wxPanel *panel          = new wxPanel(parent,-1);
        wxFlexGridSizer *sizer = new wxFlexGridSizer(2);


        //-- PATHOLOGY WIDGETS
        
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
    sizer->Add(new wxStaticText(panel,-1,_T("  - - -  PATHOLOGY EXTRACTION - - - ")));
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        wxButton *btnP1 = new wxButton(panel,-1,_T("Set P1"));
        wxButton *btnP2 = new wxButton(panel,-1,_T("Set P2"));
        wxButton *btnPat = new wxButton(panel,-1,_T("Set Region"));
        wxButton *btnExtract = new wxButton(panel,-1,_T("Extract Region"));
        wxButton *btnAxis = new wxButton(panel,-1, _T("Axis"));
        sizer->Add(btnP1);
        sizer->Add(btnP2);
        sizer->Add(btnPat);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        sizer->Add(btnExtract);sizer->Add(new wxStaticText(panel,-1,_T("  ")));

        patSliderOpacity        = new wxSlider( panel, -1, 50 , 0, 100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
        sizer->Add(new wxStaticText(panel, -1,_T("Opacity Pathological VOI")));
        sizer->Add(patSliderOpacity);
        patSliderMarchingCubes  = new wxSlider( panel, -1, 254 , 0, 1000, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
        sizer->Add(new wxStaticText(panel, -1,_T(" Marching Cubes Level Pathological VOI")));
        sizer->Add(patSliderMarchingCubes);
        sizer->Add(btnAxis);sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        //-- PATHOLOGY WIDGETS


        
        //-- STL WIDGETS
        
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
    sizer->Add(new wxStaticText(panel,-1,_T("  - - -  STL Diego - - - ")));
    sizer->Add(new wxStaticText(panel,-1,_T("  ")));
        stlSliderDeltaGauss     = new wxSlider( panel, -1, 100 , 0, 300 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
        stlSliderMarchingCubes= new wxSlider( panel, -1, 128 , 0, 256 , wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS );
        
        sizer->Add(new wxStaticText(panel,-1,_T(" Delta Gauss")));
        sizer->Add(stlSliderDeltaGauss);

        sizer->Add(new wxStaticText(panel,-1,_T(" Marching Cubes Level")));
        sizer->Add(stlSliderMarchingCubes);

        stlSliderOpacityInternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);
        stlSliderOpacityExternal = new wxSlider(panel, -1, 100,0,100, wxDefaultPosition, wxDefaultSize, wxSL_HORIZONTAL | wxSL_LABELS);

        sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL Internal")));
        sizer->Add(stlSliderOpacityInternal);

        sizer->Add(new wxStaticText(panel, -1,_T(" Opacity STL External")));
        sizer->Add(stlSliderOpacityExternal);

        wxButton *btnFileSTL = new wxButton(panel,-1,_T("Generate STL files"));
        sizer->Add(btnFileSTL);
        sizer->Add(new wxStaticText(panel,-1,_T("  ")));
         //-- STL WIDGETS
        

        
        
        
 



        panel->SetSizer(sizer);
        panel->SetAutoLayout(true);
        panel->SetSize(400,600);
        panel->Layout();


        // -- PATHOLOGY CONNECT WIDGETS
        Connect(btnP1->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP1);
        Connect(btnP2->GetId(),wxEVT_COMMAND_BUTTON_CLICKED, (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetP2);
        Connect(btnPat->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnSetPat);
        Connect(btnExtract->GetId(),wxEVT_COMMAND_BUTTON_CLICKED,(wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractPat);
        Connect(patSliderOpacity->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatOpacity);
        Connect(patSliderMarchingCubes->GetId(), wxEVT_COMMAND_SLIDER_UPDATED   , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangePatMarchingCubes);
        Connect(btnAxis->GetId(), wxEVT_COMMAND_BUTTON_CLICKED  , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnExtractAxis);
                

        // -- PATHOLOGY CONNECT WIDGETS
        
        
        // -- STL CONNECT WIDGETS
        Connect(btnFileSTL->GetId()                             , wxEVT_COMMAND_BUTTON_CLICKED  , (wxObjectEventFunction) &wxPathologyWidget_01::OnBtnFileSTL           );
        Connect(stlSliderDeltaGauss->GetId()    , wxEVT_SCROLL_THUMBRELEASE             , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLGaussLevel  );
        Connect(stlSliderMarchingCubes->GetId() , wxEVT_SCROLL_THUMBRELEASE             , (wxObjectEventFunction) &wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel);
        Connect(stlSliderOpacityInternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLInternal);
        Connect(stlSliderOpacityExternal->GetId(), wxEVT_COMMAND_SLIDER_UPDATED , (wxObjectEventFunction) &wxPathologyWidget_01::OnOpacitySTLExternal);
        // -- STL CONNECT WIDGETS

        

        

        return panel;
}
//------------------------------------------------------------------------
void wxPathologyWidget_01::Refresh()
{
        _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->Render();
}

//------------------------------------------------------------------------
void wxPathologyWidget_01::ConfigureVTK()
{
    wxBusyCursor wait;

        vtkImageData    *imagedata              = _mar->_experiment->getDynData( )->getVolume( )->castVtk();
        //-Maracas-
        _maracasSurfaceWidget->ShowMARACASData( _mar );

        //CONFIGURACION ADICIONAL
    this->ConfigureSTL();
        this->ConfigurePathologyExtraction();
        
}


void wxPathologyWidget_01::ConfigureVTK(vtkImageData *imagedata, int x, int y, int z, double param){

        this->px = x;
        this->py = y;
        this->pz = z;

        wxBusyCursor wait;
        _maracasSurfaceWidget->ShowMARACASData( _mar );


        //CONFIGURACION ADICIONAL
        this->ConfigureMPR();
        this->ConfigureSTL();
        this->ConfigurePathologyExtraction();
}





// ==================================================================================================
// CODIGO NUEVO ADICIONADO POR DIEGO CANTOR
// ==================================================================================================


// ------------------------------------------------------------------------
// START MPR FUNCTIONS - DHC
// ------------------------------------------------------------------------

void wxPathologyWidget_01::ConfigureMPR(){

        mprSphereSource = vtkSphereSource::New();
    mprMapper = vtkPolyDataMapper::New();
    mprActor = vtkActor::New();

        mprMapper->SetInput(mprSphereSource->GetOutput());
        mprActor->SetMapper(mprMapper);
        mprActor->PickableOff( );

        mprSphereSource->SetCenter(this->px,this->py,this->pz);
        mprSphereSource->SetRadius(2.5);

        mprActor->GetProperty()->SetColor( 1.0, 0.0, 0.0 );
    vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
        ren->AddActor(mprActor);

}

// ------------------------------------------------------------------------
// END MPR FUNCTIONS - DHC
// ------------------------------------------------------------------------



// ------------------------------------------------------------------------
// START STL FUNCTIONS - DHC
// ------------------------------------------------------------------------

void wxPathologyWidget_01::ConfigureSTL()
{
        stlExterna = vtkPolyData::New();
        stlInterna = vtkPolyData::New();

        dsm1 = vtkPolyDataMapper ::New();
    dsm1->SetInput (stlInterna); 
    dsm1->ScalarVisibilityOff();

    actorInternal = vtkActor::New();
    actorInternal->SetMapper (dsm1);
    actorInternal->GetProperty()->SetColor (0,1,0);

    dsm2 = vtkPolyDataMapper ::New();
    dsm2->SetInput (stlExterna);
    dsm2->ScalarVisibilityOff();

    actorExternal= vtkActor::New();
    actorExternal->SetMapper (dsm2);
    actorExternal->GetProperty()->SetRepresentationToWireframe();
    vtkRenderer *ren = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();
    ren->AddActor(actorInternal);
    ren->AddActor(actorExternal);

        stlExtractor = new vtkSTLExtractor();  
}

void wxPathologyWidget_01::generateSTLSurfaces()
{
        stlExtractor->setVolume(stlImageData);
        stlExtractor->setSigmaLevel(stlDeltaGaussLevel);
        stlExtractor->setMarchingCubesLevel(stlMarchingCubesLevel);
        stlExtractor->calculate();
        
        stlInterna->DeepCopy(stlExtractor->getInnerSurface());
        stlExterna->DeepCopy(stlExtractor->getOuterSurface());
}


void wxPathologyWidget_01::OnOpacitySTLExternal(wxScrollEvent& event){
        double value = ((double)stlSliderOpacityExternal->GetValue())/100;
    actorExternal->GetProperty( )->SetOpacity( value );
        Refresh();
}


void wxPathologyWidget_01::OnOpacitySTLInternal(wxScrollEvent& event){
        double value = ((double)stlSliderOpacityInternal->GetValue())/100;
    actorInternal->GetProperty( )->SetOpacity( value );
        Refresh();
}

void wxPathologyWidget_01::OnBtnFileSTL(wxCommandEvent& event)
{
        wxBusyCursor wait;
        wxString dirSTL = _mar->_parameters->getStringParam( 
    marParameters::e_installation_directory ); 
        dirSTL = ( dirSTL == _T("NO_DIRECTORY") ) ? wxGetHomeDir( ) : dirSTL;
        wxDirDialog dialog( this, _T("Choose a directory..."), ( !dirSTL.IsEmpty( ) )?
        dirSTL: wxGetHomeDir( ) );

        if( dialog.ShowModal( ) == wxID_OK ) 
        {
           
        
                // ------------------------------------------------------------------------
                //  1.  GENERATE STL FILES
                // ------------------------------------------------------------------------
                const char* fileprefix = "c:\\Creatis\\";
                std::string prefix = fileprefix;
                

                // 1.1. Se hace un filtro triangular puesto que el stl writer solo recibe poligonos triangulares.

        vtkTriangleFilter *filtro = vtkTriangleFilter::New();
                filtro->SetInput(stlInterna);
                vtkPolyDataConnectivityFilter *pdcf = vtkPolyDataConnectivityFilter::New();
        pdcf->SetInput( filtro->GetOutput() );
        vtkClosePolyData *cpd = vtkClosePolyData::New();
        cpd->SetInput( pdcf->GetOutput() );

                // 1.2 se escribe a disco el archivo stl de la superficie interna
        cpd->Update();
        vtkSTLWriter *writerI = vtkSTLWriter::New();
        writerI->SetInput( cpd->GetOutput() );
        prefix = fileprefix;
        prefix += "internal.stl";
        writerI->SetFileName(prefix.c_str());
        writerI->SetFileTypeToASCII();
        writerI->Write();
        writerI->Delete();

                // 1.3 se escribe a disco el archivo stl de la superficie externa
                filtro->SetInput(stlExterna);
        cpd->Update();
        vtkSTLWriter *writerE = vtkSTLWriter::New();
        writerE->SetInput( cpd->GetOutput() );
        prefix = fileprefix;
        prefix += "external.stl";
        writerE->SetFileName( prefix.c_str() );
        writerE->SetFileTypeToASCII();
        writerE->Write();
        writerE->Delete();
   
        filtro->Delete();
        cpd->Delete();
        pdcf->Delete();
        }

        //By default *always* update e_installation_directory:
        _mar->_parameters->setStringParam( marParameters::e_installation_directory, dialog.GetPath( ) ); 
        _mar->saveParameters( );
}


void wxPathologyWidget_01::OnChangeSTLGaussLevel(wxScrollEvent& event)
{
        stlDeltaGaussLevel  = ((double)stlSliderDeltaGauss->GetValue())/100;
        generateSTLSurfaces();
        Refresh();
}


void wxPathologyWidget_01::OnChangeSTLMarchingCubesLevel(wxScrollEvent& event)
{
        stlMarchingCubesLevel = ((double)stlSliderMarchingCubes->GetValue());
        generateSTLSurfaces();
        Refresh();
        
}


// ------------------------------------------------------------------------
// END STL FUNCTIONS - DHC
// ------------------------------------------------------------------------


// ------------------------------------------------------------------------
// START PATHOLOGY FUNCTIONS - DHC
// ------------------------------------------------------------------------

double GetDistanciaEuclideana(double* a, double* b){
        return sqrt(((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2])));
}


double GetValorPromedio(vtkImageData *data){
        int ext[6];
        int i, j, k;
        int numVoxels = 0;
        double promedio = 0;
        double valor;
        
        
        data->GetExtent(ext);

        
        for(i=ext[0];i<=ext[1];i++){
                for(j=ext[2];j<=ext[3];j++){
                        for(k=ext[4];k<=ext[5];k++){
                                valor =data->GetScalarComponentAsDouble(i,j,k,0);
                                promedio += valor;
                                numVoxels++;
                        }
                }
        }
        return (double)promedio/(double)numVoxels;
}


void wxPathologyWidget_01::ConfigurePathologyExtraction()
{
}

void wxPathologyWidget_01::generatePathologySurface(){
        // 1. Inicializacion
}


void wxPathologyWidget_01::OnBtnSetP1(wxCommandEvent &event){
        vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();
        
        if(p1Actor){
                sw->GetRenderer()->RemoveActor(p1Actor);
        }
        
        
        sw->GetSphereCenter(p1Center);

        p1SphereSource = vtkSphereSource::New();
    p1Mapper = vtkPolyDataMapper::New();
    p1Actor = vtkActor::New();

        p1Mapper->SetInput(p1SphereSource->GetOutput());
        p1Actor->SetMapper(p1Mapper);
        p1Actor->PickableOff( );

        p1SphereSource->SetCenter(p1Center[0],p1Center[1],p1Center[2]);
        p1SphereSource->SetRadius(0.5);

        p1Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);

        sw->GetRenderer()->AddActor(p1Actor);
        Refresh();

}

void wxPathologyWidget_01::OnBtnSetP2(wxCommandEvent &event){

        vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();


        if(p2Actor){
                sw->GetRenderer()->RemoveActor(p1Actor);
        }

        sw->GetSphereCenter(p2Center);

        p2SphereSource = vtkSphereSource::New();
    p2Mapper = vtkPolyDataMapper::New();
    p2Actor = vtkActor::New();

        p2Mapper->SetInput(p2SphereSource->GetOutput());
        p2Actor->SetMapper(p2Mapper);
        p2Actor->PickableOff( );

        p2SphereSource->SetCenter(p2Center[0],p2Center[1],p2Center[2]);
        p2SphereSource->SetRadius(0.5);

        p2Actor->GetProperty()->SetColor( 0.0, 1.0, 0.0);

        sw->GetRenderer()->AddActor(p2Actor);
        Refresh();
}




void wxPathologyWidget_01::OnBtnSetPat(wxCommandEvent &event){

        if (p2SphereSource == NULL || p1SphereSource == NULL){
                return;
        }

        vtk3DSurfaceWidget *sw = _maracasSurfaceWidget->GetVtk3DSurfaceWidget();

        if(patActor){
                sw->GetRenderer()->RemoveActor(patActor);
        }

        
        double centro[3];

        
        centro[0] = (p1Center[0] + p2Center[0]) / 2;
        centro[1] = (p1Center[1] + p2Center[1]) / 2;
        centro[2] = (p1Center[2] + p2Center[2]) / 2;

        patSphereSource = vtkSphereSource::New();
    patMapper = vtkPolyDataMapper::New();
    patActor = vtkActor::New();

        patMapper->SetInput(patSphereSource->GetOutput());
        patActor->SetMapper(patMapper);
        

        patSphereSource->SetCenter(centro[0],centro[1],centro[2]);
        patSphereSource->SetRadius(GetDistanciaEuclideana(p1Center,p2Center));

        patActor->GetProperty()->SetColor( 0.15, 0.75, 0.15 );
        patActor->PickableOff( );
        patActor->GetProperty( )->SetOpacity(0.15);

        sw->GetRenderer()->AddActor(patActor);
        Refresh();
}


        

void wxPathologyWidget_01::OnBtnExtractPat(wxCommandEvent &event){

        //RENDERER
        //----------------------------------------------------------------------------------
        vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();

        //IMAGEN ORIGINAL
        //----------------------------------------------------------------------------------
        vtkImageData* originalData = vtkImageData::New();
        originalData->DeepCopy(_mar->_experiment->getDynData()->getVolume()->castVtk());
        
        double spacing[3];
        originalData->GetSpacing(spacing);
        
        int extend[6];
        originalData->GetExtent(extend);        

        

        //==================================================================================
        // REGION DE INTERES EN LA IMAGEN ORIGINAL: patImageData
        //==================================================================================

        double p1crem[3], p2crem[3];

        p1crem[0] = p1Center[0] / spacing[0];
        p1crem[1] = p1Center[1] / spacing[1];
        p1crem[2] = p1Center[2] / spacing[2];

        p2crem[0] = p2Center[0] / spacing[0];
        p2crem[1] = p2Center[1] / spacing[1];
        p2crem[2] = p2Center[2] / spacing[2];

        

        

    //COORDENADAS DE LA IMAGEN 

        double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
        double centro[3];
        
        centro[0] = (p1crem[0] + p2crem[0]) / 2;
        centro[1] = (p1crem[1] + p2crem[1]) / 2;
        centro[2] = (p1crem[2] + p2crem[2]) / 2;

        

        //Obtener los ejes sobre la esfera para delimitar el area.

        double lowX     = centro[0] - (2*radio);
        double highX    = centro[0] + (2*radio);
        double lowY     = centro[1] - (2*radio);
        double highY    = centro[1] + (2*radio);
        double lowZ     = centro[2] - (2*radio);
        double highZ    = centro[2] + (2*radio);


        
        //Extraer la region de interes:patImageData
        //----------------------------------------------------------------------------------

        voiFilter = vtkExtractVOI::New();
        voiFilter->SetInput(originalData);
        voiFilter->SetVOI((int)lowX, (int)highX, (int)lowY, (int)highY, (int)lowZ, (int)highZ);
        voiFilter->Update();
        patImageData = vtkImageData::New();
        patImageData->DeepCopy(voiFilter->GetOutput());
//      patImageData->SetSpacing(spacing);

        //==================================================================================
        // REGION DE INTERES ESPECIFICA
        //==================================================================================

        vtkOtsuSphereSource *otsu = new vtkOtsuSphereSource();
        thresholdOTSU = otsu->calculateOptimalThreshold(patImageData);


        //==================================================================================
        // MARCHING CUBES
        //==================================================================================

        patMCLevel = 128.0;
        cubesFilter = vtkMarchingCubes::New();
        cubesFilter->SetInput(otsu->getImageForSegmentation());
        cubesFilter->SetValue(0,patMCLevel); 
        cubesFilter->SetNumberOfContours(1);

        //==================================================================================
        // POLYDATA
        //==================================================================================

    polyUnico = vtkPolyDataConnectivityFilter::New();
        polyUnico->SetInput(cubesFilter->GetOutput());
        polyUnico->SetExtractionModeToLargestRegion();

        stripper = vtkStripper::New();
        stripper->SetInput(polyUnico->GetOutput());


        //==================================================================================
        // DATA MAPPING
        //==================================================================================
        
        dataMapper = vtkPolyDataMapper::New();
        dataMapper->SetInput(stripper->GetOutput());
        dataMapper->ScalarVisibilityOff();

        //==================================================================================
        // ACTORES
        //==================================================================================

        
        dataActor = vtkActor::New();
        dataActor->SetMapper(dataMapper);
        dataActor->SetPosition(lowX*spacing[0],lowY*spacing[1],lowZ*spacing[2]);
        dataActor->GetProperty()->SetOpacity(patOpacityLevel);
        dataActor->GetProperty()->SetColor( 0.20, 0.20, 1.0 );
        dataActor->GetProperty()->SetRepresentationToWireframe();

        renderer->AddActor(dataActor);
        
        
        
        outlineFilter = vtkOutlineFilter::New();
    outlineFilter->SetInput(voiFilter->GetOutput());

        vtkPolyDataMapper* outlineMapper = vtkPolyDataMapper::New();
        outlineMapper->SetInput(outlineFilter->GetOutput());
                
        if(outlineActor){
                renderer->RemoveActor(outlineActor);
        }

        if(patActor){
                renderer->RemoveActor(patActor);
        }

        if(p3Actor){
                renderer->RemoveActor(p3Actor);
        }

        
        outlineActor = vtkActor::New();
        outlineActor->SetMapper(outlineMapper);
        outlineActor->GetProperty()->SetOpacity(1);


        /*p3SphereSource = vtkSphereSource::New();
    p3Mapper = vtkPolyDataMapper::New();
    p3Actor = vtkActor::New();

        p3Mapper->SetInput(p3SphereSource->GetOutput());
        p3Actor->SetMapper(p3Mapper);
        p3Actor->PickableOff( );

        p3SphereSource->SetCenter(p2Center[0]-(lowX*spacing[0]), p2Center[1]-(lowY*spacing[1]) , p2Center[2]-(lowZ*spacing[2]) );
        p3SphereSource->SetRadius(0.5);
        p3Actor->GetProperty()->SetColor( 0.9, 0.4, 0.2);*/

        
        //renderer->AddActor(outlineActor);
        //renderer->AddActor(p3Actor);
        
        
        Refresh();
}


void wxPathologyWidget_01::OnChangePatOpacity(wxScrollEvent& event){
        patOpacityLevel = ((double)patSliderOpacity->GetValue())/100.0;
    dataActor->GetProperty( )->SetOpacity(patOpacityLevel);
        Refresh();
}

void wxPathologyWidget_01::OnChangePatMarchingCubes(wxScrollEvent& event){
        patMCLevel = ((double)patSliderMarchingCubes->GetValue());
        cubesFilter->SetValue(0,patMCLevel);
        cubesFilter->Update();
        Refresh();
}



//Metodo auxiliar que devuelve la distancia minima de un punto a un polydata.
//Se evalua la distancia del punto a todos los centros de las celdas del polydata y se devuelve la menor
double GetDistanciaAPolydata(double* punto, vtkPoints* centers){
        double distancia = 0.0;
        double minima = 1e100;
        double vertice[3];
        int N = centers->GetNumberOfPoints();

        for(int i = 0; i < N; i++){
                 centers->GetPoint(i,vertice);
                 distancia =  sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));

                 if (distancia < minima){
                         minima = distancia;
                 }
        }
        return minima;
}


double GetProductoPunto(double* vectorA, double* vectorB){
        return (vectorA[0]*vectorB[0] +  vectorA[1]*vectorB[1] + vectorA[2]*vectorB[2]);
}


double GetNormaVector(double* vector){
        return sqrt(GetProductoPunto(vector, vector));
}

//construye el vector como la diferencia de dos puntos (vectores en el origen)
void GetVector(double* p0, double* p1, double* resultado){ //p1:destino - p0:origen
        
        resultado[0] = p1[0] - p0[0];
        resultado[1] = p1[1] - p0[1];
        resultado[2] = p1[2] - p0[2];

}


void VectorUnitario(double *v){
        double norma = GetNormaVector(v);
        v[0] = v[0] / norma;
        v[1] = v[1] / norma;
        v[2] = v[2] / norma;
}

//calculo de las normales
void GetNormalCelda(double* p0, double* p1, double* p2, double* normal){
        
        double U[3],V[3];
        GetVector(p1,p2,U); //p2-p1
        GetVector(p1,p0,V); //p0-p1
        

        normal[0] = (U[1]*V[2])-(U[2]*V[1]);
        normal[1] = -((U[0]*V[2])-(U[2]*V[0]));
        normal[2] = (U[0]*V[1])-(U[1]*V[0]);
        VectorUnitario(normal);
}

void GetNormal(vtkPoints* puntosCelda, double* normal){
        double p0[3], p1[3], p2[3];

        puntosCelda->GetPoint(0, p0);
        puntosCelda->GetPoint(1, p1);
        puntosCelda->GetPoint(2, p2);

        GetNormalCelda(p0,p1,p2, normal);
}


// Obtiene el ID de la celda mas cercana
int GetClosestCellOnPolyData(double* punto, vtkPoints* centers){
        double distancia = 0.0;
        double minima = 1e100;
        int ID = -1;
        double vertice[3];
        int N = centers->GetNumberOfPoints();

        for(int i = 0; i < N; i++){
                 centers->GetPoint(i,vertice);
                 distancia =  sqrt(((punto[0]-vertice[0])*(punto[0]-vertice[0]))+((punto[1]-vertice[1])*(punto[1]-vertice[1]))+((punto[2]-vertice[2])*(punto[2]-vertice[2])));

                 if (distancia < minima){
                         minima = distancia;
                         ID = i;
                 }
        }
        
        //printf("distancia del punto p: (%.2f,%.2f,%.2f) al polydata es %.2f (ID %i)\n", punto[0], punto[1], punto [2], minima, ID);
        return ID;
}

double GetProjectionOnNormal(double* vector, double* normal){
        double pp = GetProductoPunto(vector, normal);
        double vv = GetNormaVector(normal);
        double pr = pp / vv;
        return pr;
}





void wxPathologyWidget_01::OnBtnExtractAxis(wxCommandEvent& event){



        //================================================================================
        // RENDERER
        //================================================================================
        vtkRenderer* renderer = _maracasSurfaceWidget->GetVtk3DSurfaceWidget()->GetRenderer();

        //===================================================================================
        // IMAGEN ORIGINAL.
        //===================================================================================
        int extentOriginal[6];
        this->patImageData->GetExtent(extentOriginal);  
        

        //===================================================================================
        // REMUESTREO
        //===================================================================================

        int factor = 1;
        double resampling = 1.0/(double)factor; 

        

        // REMUESTREO DE LA IMAGEN ORIGINAL
        //-----------------------------------------------------------------------------------
        vtkImageResample* resamplingFilter = vtkImageResample::New();
        resamplingFilter->SetAxisOutputSpacing(0,resampling);
        resamplingFilter->SetAxisOutputSpacing(1,resampling);
        resamplingFilter->SetAxisOutputSpacing(2,resampling);
        resamplingFilter->ReleaseDataFlagOff();
        resamplingFilter->SetInput(patImageData);
        resamplingFilter->SetInterpolationModeToCubic();
        resamplingFilter->Update();


        //IMAGEN ORIGINAL REMUESTREADA
        //-----------------------------------------------------------------------------------
        vtkImageData* imagenRemuestreada = resamplingFilter->GetOutput();
        imagenRemuestreada->ReleaseDataFlagOff();
        


        //===================================================================================
        // IMAGEN BINARIA = CRECIMIENTO DE REGIONES
        //===================================================================================

        this->isoBinaria = vtkImageData::New();
        this->isoBinaria->CopyStructure(resamplingFilter->GetOutput());
        this->isoBinaria->SetScalarTypeToDouble();
        
        
        
        //CALCULO DE LAS DIMENSIONES Y EL EXTENT
        //-----------------------------------------------------------------------------------
        int extentBinaria[6];
        this->isoBinaria->GetExtent(extentBinaria);
        int *dim = isoBinaria->GetDimensions();

        
        
        // IMAGEN EN NEGRO.
        //-----------------------------------------------------------------------------------
        int eX, eY, eZ, idpoint; // extent
        double *ptr = NULL;

        FILE* logger = freopen("tests.txt","w",stdout);



        for(eX = extentBinaria[0]; eX <= extentBinaria[1]; eX++){
                for(eY = extentBinaria[2]; eY <= extentBinaria[3]; eY++){
                        for (eZ = extentBinaria[4]; eZ <= extentBinaria[5]; eZ++){
                                idpoint = isoBinaria->FindPoint(eX, eY, eZ);
                                ptr = (double *)isoBinaria->GetScalarPointer(eX,eY,eZ); // El apuntador se obtiene sobre el extent.
                                *ptr = 0;
                        }
                }
        }


    //===================================================================================
        // POLYDATA
        //===================================================================================
     vtkPolyData *polyData= stripper->GetOutput();  //OJO SE TOMA SOLO LA SUPERFICIE MAS GRANDE (connectivity + stripper)
         //===================================================================================

         vtkPoints *puntosPolyData = polyData->GetPoints();
         vtkPointData *pointData = polyData->GetPointData();
        
         vtkPolyDataWriter* writePoly = vtkPolyDataWriter::New();
         writePoly->SetFileName("poly.vtk");
         writePoly->SetInput(polyData);
         writePoly->Write();
         
         
         
        // PROCESAR LOS VERTICES
        //----------------------------------------------------------------------------------

         int numVertices = puntosPolyData->GetNumberOfPoints();
    
         // mover el sistema de referencias del polydata para poder determinar lo que se encuentra
         //por dentro segun las semillas.
         int j;
         double vertice[3];
         /*for(j=0;j<numVertices;j++){
                 puntosPolyData->GetPoint(j,vertice);
                 vertice[0] = vertice[0] + extentOriginal[0];
                 vertice[1] = vertice[1] + extentOriginal[2];
                 vertice[2] = vertice[2] + extentOriginal[4];
                 puntosPolyData->SetPoint(j,vertice);
         }*/


        
        // PROCESAR LAS CELDAS
        //----------------------------------------------------------------------------------
        

         //1. Obtener normales de la superficie (celdas)
         vtkPolyDataNormals* polyNormalsFilter = vtkPolyDataNormals::New(); //estas son las normales de los vertices NO de las celdas
         polyNormalsFilter->SetInput(polyData);
         polyNormalsFilter->SplittingOn();
         polyNormalsFilter->ConsistencyOn();
         polyNormalsFilter->NonManifoldTraversalOn();
         polyNormalsFilter->ComputeCellNormalsOn();
         polyNormalsFilter->Update();
        
         
         vtkPolyData* polyDataConNormales = polyNormalsFilter->GetOutput();
     

         //2. Obtener los centros de las celdas
         // Necesitamos el centro de las celdas para trazar el vector interno de un voxel 
         //a una cara de la superficie.
         vtkCellCenters* cellCentersFilter = vtkCellCenters::New();
         cellCentersFilter->SetInput(polyDataConNormales);
         cellCentersFilter->VertexCellsOn(); 
         cellCentersFilter->Update();
         
         vtkPolyData* polydataCenters = cellCentersFilter->GetOutput();
         vtkPoints*   polyCenterPoints  = polydataCenters->GetPoints();

         vtkPoints*   puntosCelda = NULL;
         vtkIdList*   lista = NULL;
         vtkCell*     celda = NULL;
     double normal[3]; 
         int numCeldas = polyCenterPoints->GetNumberOfPoints();


         vtkPoints *normales = vtkPoints::New();//almacena los componentes cartesianos de la normal de cada celda
         
     int i;     
     for(i = 0; i < numCeldas ; i++){
                
                polyCenterPoints->GetPoint(i,vertice);
        
            celda = polyDataConNormales->GetCell(i);
                puntosCelda = celda->GetPoints();
        GetNormal(puntosCelda, normal);
                normales->InsertPoint(i, normal);

                                
                lista = celda->GetPointIds();
                for(j = 0; j < lista->GetNumberOfIds(); j++){
                        int idPunto = lista->GetId(j);
                        puntosCelda->GetPoint(j, vertice);
                        
                }

         }

        
         
         //==========================================================================================
         //IMPLEMENTAR EL CRECIMIENTO DE REGIONES
         //==========================================================================================
         
        double spacing[3];
        patImageData->GetSpacing(spacing);
        
        double p1crem[3], p2crem[3];

        p1crem[0] = p1Center[0] / spacing[0];
        p1crem[1] = p1Center[1] / spacing[1];
        p1crem[2] = p1Center[2] / spacing[2];

        p2crem[0] = p2Center[0] / spacing[0];
        p2crem[1] = p2Center[1] / spacing[1];
        p2crem[2] = p2Center[2] / spacing[2];

        

        

    //COORDENADAS DE LA IMAGEN 

        double radio = GetDistanciaEuclideana(p1crem,p2crem) / 2;
        double centro[3];
        
        centro[0] = (p1crem[0] + p2crem[0]) / 2;
        centro[1] = (p1crem[1] + p2crem[1]) / 2;
        centro[2] = (p1crem[2] + p2crem[2]) / 2;

        

        //Obtener los ejes sobre la esfera para delimitar el area.

        double lowX             = centro[0] - (2*radio);
        double highX    = centro[0] + (2*radio);
        double lowY             = centro[1] - (2*radio);
        double highY    = centro[1] + (2*radio);
        double lowZ             = centro[2] - (2*radio);
        double highZ    = centro[2] + (2*radio);

        //================================================================================
        //0. Inicializacion
        //================================================================================
        
        int numeroPuntos = dim[0] * dim[1] * dim[2];


        int sem[3];


        sem[0] = (int)(p2Center[0]-(lowX*spacing[0]));
        sem[1] = (int)(p2Center[1]-(lowY*spacing[1]));
        sem[2] = (int)(p2Center[2]-(lowZ*spacing[2]));

                
        int idSemilla = isoBinaria->FindPoint(sem[0],sem[1],sem[2]);

        if (idSemilla == -1) return;

        int graphSize = (extentBinaria[1]-extentBinaria[0]) * (extentBinaria[3]-extentBinaria[2])*(extentBinaria[5]-extentBinaria[4]);

                
        vtkPriorityQueue *colaEvaluacion = vtkPriorityQueue::New();
                
        colaEvaluacion->Allocate(graphSize);
                
        colaEvaluacion->Insert(0, idSemilla);

        vtkFloatingPointType prioridad;
                
        int longitudCola = colaEvaluacion->GetNumberOfItems();

                
        int vecinoID = -1;

        int puntoID = -1;

        double* coordsVecino;
        
        double* coordsPunto;

        int cont=0;
                
        vtkPoints *resultado = vtkPoints::New();

        int premiere = 1;


        double centroCelda[3], vectorInterno[3], vectorNormal[3];
        
        while(longitudCola > 0){

                        //================================================================================
                        //1. Procesar el punto
                        //================================================================================
                                                
                    puntoID = colaEvaluacion->Pop(0,prioridad);

                        coordsPunto = isoBinaria->GetPoint(puntoID);
                        
                        
                    ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsPunto[0]*factor), (int)(coordsPunto[1]*factor), (int)(coordsPunto[2]*factor));
                        *ptr = 255.0;

                        if(premiere == 1){
                                vtkImageWriter* writeImage = vtkImageWriter::New();
                                writeImage->SetInput(isoBinaria);
                                writeImage->SetFileName("sem.vtk");
                                writeImage->Write();
                                premiere = 0;
                        }

                        resultado->InsertPoint(cont, coordsPunto);
                        cont++;
                
                        //printf("iniciando en el punto (%.1f,%.1f,%.1f) \n", coordsPunto[0], coordsPunto[1], coordsPunto [2]);

                        //================================================================================
                        //2. Evaluar los vecinos del punto
                        //================================================================================
                        int vk,vj,vi;
                        for(vk = -1; vk<2; vk++) {
                                for(vj = -1; vj<2; vj++) {
                                        for(vi = -1; vi<2; vi++) {         
    
                                                vecinoID = puntoID + (vk * dim[1]*dim[0]) + (vj * dim[0]) + vi;
                                                
                                                if( vecinoID >= 0 && vecinoID < numeroPuntos && vecinoID != 0 && vecinoID != puntoID ) {
                                                
                                                                coordsVecino = isoBinaria->GetPoint(vecinoID);

                                                                //ptr = (unsigned char *)isoBinaria->GetScalarPointer(coordsVecino[0],coordsVecino[1],coordsVecino[2]);
                                                                ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));
                                                                
                                                                //================================================================================
                                                                // Para cada vecino no evaluado evaluar la distancia al polydata
                                                                //================================================================================
                                                                if (*ptr == 0){
                                                                        //printf("evaluando punto (%.1f,%.1f,%.1f) \n", coordsVecino[0], coordsVecino[1], coordsVecino [2]);
                                                                        
                                                                        int normalID = GetClosestCellOnPolyData(coordsVecino, polyCenterPoints);
                                                                        
                                                                        
                                                                        polyCenterPoints->GetPoint(normalID,centroCelda);
                                                                        normales->GetPoint(normalID, vectorNormal);

                                                                        GetVector(coordsVecino, centroCelda, vectorInterno);
                                                                        VectorUnitario(vectorInterno);
                                                                        
                                                                        double projection = GetProjectionOnNormal(vectorInterno, vectorNormal);
                                                                        ptr = (double *)isoBinaria->GetScalarPointer((int)(coordsVecino[0]*factor),(int)(coordsVecino[1]*factor),(int)(coordsVecino[2]*factor));

                                                                        if (projection > 0){

                                                                        //================================================================================
                                                                        //2.1 Si la proyeccion es positiva
                                                                        //colocar el vecino en la cola y poner el punto correspondiente en la imagen a 255.
                                                                        // (no ha pasado la frontera, el marcado a 255 se hace arriba en el while, ver paso 1)
                                                                        //================================================================================
                                                                                colaEvaluacion->Insert(prioridad,vecinoID);

                                                                        } else {
                                                                        
                                                                        //================================================================================      
                                                                        //2.2 Si la proyeccion es negativa con la normal mas cercana del polydata,
                                                                        //descartar el vecino == no poner en la cola == alcanzamos la superficie
                                                                        //================================================================================
                                                                        
                                                                                *ptr = 0;
                                                                        }
                                                                }
                                                        
                                                }
                                        }
                                }
                        }

                        longitudCola = colaEvaluacion->GetNumberOfItems();
        }
        colaEvaluacion->Delete();
                
        
        //================================================================================
        //  MAPA DE DISTANCIAS
        //================================================================================

        this->isoDistance = vtkImageEuclideanDistance::New();
        this->isoDistance->SetInput(this->isoBinaria);
        this->isoDistance->InitializeOn();
        this->isoDistance->Update();
        vtkImageData *imDistancias = isoDistance->GetOutput();

        
        //================================================================================
        //  CAMINO MINIMO SOBRE EL MAPA DE DISTANCIAS
        //================================================================================
        this->dijkstraFilter = vtkDijkstraImageData::New();
        this->dijkstraFilter->SetInput(isoDistance->GetOutput());

        int idA = (int)isoBinaria->FindPoint(p1Center[0],p1Center[1], p1Center[2]);
        int idB = (int)isoBinaria->FindPoint(p2Center[0],p2Center[1], p2Center[2]);

        dijkstraFilter->SourceID = idA;
        dijkstraFilter->SinkID = idB;
        dijkstraFilter->Update();
        
        caminoMapper = vtkPolyDataMapper::New();
        caminoMapper->SetInput(dijkstraFilter->GetOutput());
        caminoMapper->ScalarVisibilityOff();

        renderer->RemoveActor(dataActor);

        //=======================================================================================================       
        //RECORRIDO DE LOS PUNTOS DEL EJE DETECTADO
        //=======================================================================================================       
        vtkIdList *camino = dijkstraFilter->GetShortestPathIdList();

        int size = camino->GetNumberOfIds();
        double coords[3];
 

        for(i=0;i < size;i++){
                int pointId = camino->GetId(i);
                isoBinaria->GetPoint(pointId, coords);
                double scalar = imagenRemuestreada->GetPointData()->GetScalars()->GetTuple1(pointId);
                        
                //----------------------------------------------------------------------------------------
                //Visualizacion del nodo del eje
                //----------------------------------------------------------------------------------------
                vtkSphereSource   *_nodo   = vtkSphereSource::New();
                _nodo->SetCenter(coords[0],coords[1],coords[2]);
                _nodo->SetRadius(0.2);

                vtkPolyDataMapper *_mapper = vtkPolyDataMapper::New();
                _mapper->SetInput(_nodo->GetOutput());


        vtkActor          *_actor  = vtkActor::New();
                _actor->GetProperty()->SetColor( 0.8, 0.2, 0.3);
            _actor->SetMapper(_mapper);
                _actor->PickableOff( );
                
                renderer->AddActor(_actor);
                
        
        }


}
// ------------------------------------------------------------------------
// END PATHOLOGY FUNCTIONS - DHC
// ------------------------------------------------------------------------